text:
"""
Code for creating and manipulating chunklists (i.e., list of views into a spectrum).
For example, creating a list of views into the orders of a spectrum to be analyzed or
creating a list of views into chunks around entries in a line list.

Author: Eric Ford
Created: August 2020
Contact: https://github.com/eford/
"""

function make_chunk_list_from_loc_df(spectra::AS, inst::AbstractInstrument2D, df_λ_good::DataFrame ) where {AS<:AbstractSpectra }
    chunk_locs_df = make_order_pixel_list_for_chunks_in_order(spectra, inst, df_λ_good)
    ChunkList( map(r-> ChunkOfSpectrum(spectra,(pixels=r.pixels,order=r.order)),  eachrow(chunk_locs_df) ),
            map(r->r.order,eachrow(chunk_locs_df) ) )

end

function make_chunk_list_timeseries_from_λ_ranges( spectra::AS, df_λ_good::DataFrame ) where {ST<:AbstractSpectra, AS<:AbstractArray{ST,1} }
    times = map(s->s.metadata[:bjd],spectra)
    inst = get_inst(spectra)
    metadata = make_vec_metadata_from_spectral_timeseries(spectra)
    chunk_lists = map( spec->RvSpectMLBase.make_chunk_list_from_loc_df(spec,get_inst(spectra), df_λ_good), spectra)

    chunk_list_timeseries = ChunkListTimeseries(times, chunk_lists, inst=first(spectra).inst, metadata=metadata )
end

function make_chunk_list_from_orders_pixels_df(spectra::AS, inst::AbstractInstrument2D, df_orders_pixels::DataFrame ) where {AS<:AbstractSpectra }
    ChunkList( map(r-> ChunkOfSpectrum(spectra,(pixels=r.pixels,order=r.order)),  eachrow(df_orders_pixels) ),
            map(r->r.order,eachrow(df_orders_pixels) ) )

end

function make_chunk_list_timeseries_from_orders_pixels_df(spectra::AS, inst::AbstractInstrument2D, df_orders_pixels::DataFrame ) where {ST<:AbstractSpectra, AS<:AbstractArray{ST,1} }
    times = map(s->s.metadata[:bjd],spectra)
    inst = get_inst(spectra)
    metadata = make_vec_metadata_from_spectral_timeseries(spectra)
    chunk_lists = map( spec->RvSpectMLBase.make_chunk_list_from_orders_pixels_df(spec,get_inst(spectra), df_orders_pixels), spectra)

    chunk_list_timeseries = ChunkListTimeseries(times, chunk_lists, inst=first(spectra).inst, metadata=metadata )
end

function make_chunk_list_timeseries_telluric_free(spectra::AS, df_λ_good::DataFrame; min_pixels_in_chunk::Integer = 250, verbose::Bool = false ) where {ST<:AbstractSpectra, AS<:AbstractArray{ST,1} }
    max_Δv_match_chunk = 2*max_bc
    inst = first(spectra).inst
    times = map(s->s.metadata[:bjd],spectra)
    num_obs = length(spectra)
    metadata = make_vec_metadata_from_spectral_timeseries(spectra)
    time_series_of_chunk_lists = map(spec->RvSpectMLBase.make_chunk_list_from_loc_df(spec, inst, df_λ_good), spectra)
    @assert all(map(cl->length(cl.data),time_series_of_chunk_lists) .>= 1)
    (min_num_chunks, obs_idx_min_num_chunks) = findmin(map(obs_idx->length(time_series_of_chunk_lists[obs_idx].data),1:length(time_series_of_chunk_lists) ))
    #time_series_of_chunk_lists_output = Array{Any,1}(undef,min_num_chunks)
    #=
    # Try to average over obs_idx with the minimum number of chunks
    obs_idx_min_num_chunks = findall(map(obs_idx->length(time_series_of_chunk_lists[obs_idx].data)==min_num_chunks,1:length(time_series_of_chunk_lists) ))

    λ_min_targets = map(ch_idx->median(map(obs_idx->minimum(time_series_of_chunk_lists[obs_idx][ch_idx].λ),obs_idx_min_num_chunks)),1:min_num_chunks)
    λ_max_targets = map(ch_idx->median(map(obs_idx->maximum(time_series_of_chunk_lists[obs_idx][ch_idx].λ),obs_idx_min_num_chunks)),1:min_num_chunks)
    order_idx_targets = map(ch_idx->round(Int,median(map(obs_idx->time_series_of_chunk_lists[obs_idx][ch_idx].λ.indices[2],obs_idx_min_num_chunks))),
                                1:min_num_chunks)
    =#
    # Try just using first obs_idx with the minimum number of chunks
    λ_min_targets = map(ch_idx->minimum(time_series_of_chunk_lists[obs_idx_min_num_chunks][ch_idx].λ),1:min_num_chunks)
    λ_max_targets = map(ch_idx->maximum(time_series_of_chunk_lists[obs_idx_min_num_chunks][ch_idx].λ),1:min_num_chunks)
    order_idx_targets = map(ch_idx->time_series_of_chunk_lists[obs_idx_min_num_chunks][ch_idx].λ.indices[2],1:min_num_chunks)

    similarity_of_chunks = fill(Inf, min_num_chunks, num_obs)
    similarity_chunk_idx = zeros(Int64,min_num_chunks, num_obs)
    chunk_cant_be_this_big = 99999999999
    min_pixels_in_chunk_in_any_obs = fill(chunk_cant_be_this_big, min_num_chunks)
    for obs_idx in 1:length(time_series_of_chunk_lists)
        #=
        if length(time_series_of_chunk_lists[obs_idx].data) == min_num_chunks
            continue  # TODO: Technically should still check that they all match
        end
        =#
        #num_chunks_this_obs = length(time_series_of_chunk_lists[obs_idx].data)
        for ch_idx in 1:min_num_chunks
            if verbose    println("# Trying to match ch_idx= ",ch_idx, " lmin= ", λ_min_targets[ch_idx], " lmax= ", λ_max_targets[ch_idx])   end
            idx_order_match = findall(map(ch->time_series_of_chunk_lists[obs_idx][ch].λ.indices[2], 1:length(time_series_of_chunk_lists[obs_idx].data)).==order_idx_targets[ch_idx])
            if verbose  &&  ch_idx >= 110
                println("# idx_order_match= ", idx_order_match)
            end

            if !(length(idx_order_match) >= 1)   continue    end
            # Find if any ch_idx2 is good match to ch_idx's wavelength range
            #similarity_of_chunks .= Inf
            #similarity_chunk_idx .= 0
            for (i,ch_idx2) in enumerate(idx_order_match)
                delta_λmax = maximum(time_series_of_chunk_lists[obs_idx][ch_idx2].λ)-λ_max_targets[ch_idx]
                delta_λmin = minimum(time_series_of_chunk_lists[obs_idx][ch_idx2].λ)-λ_min_targets[ch_idx]
                if verbose    && ch_idx >= 110
                    println("# ch_idx2 = ", ch_idx2, " dlmin= ", delta_λmin/λ_max_targets[ch_idx]*speed_of_light_mps, " dlmax= ", delta_λmax/λ_max_targets[ch_idx]*speed_of_light_mps)
                end
                similarity_score = ((abs(delta_λmax)/λ_max_targets[ch_idx]*speed_of_light_mps) / max_Δv_match_chunk)^2 +
                                          ((abs(delta_λmin)/λ_max_targets[ch_idx]*speed_of_light_mps) / max_Δv_match_chunk)^2
                if similarity_score < similarity_of_chunks[ch_idx,obs_idx]
                    similarity_of_chunks[ch_idx,obs_idx] = similarity_score
                    similarity_chunk_idx[ch_idx,obs_idx] = ch_idx2
                end
                #= if (abs(delta_λmax)/λ_max_targets[ch_idx]*speed_of_light_mps < max_Δv_match_chunk) &&
                    (abs(delta_λmax)/λ_max_targets[ch_idx]*speed_of_light_mps < max_Δv_match_chunk)
                    ch_perm_idx[ch_idx] = ch_idx2
                end =#
            end # for over ch_idx2 to match target

            #ch_idx2 = similarity_chunk_idx[findmin(similarity_of_chunks)[2]]
            ch_idx2 = similarity_chunk_idx[ch_idx,obs_idx]
            delta_λmax = maximum(time_series_of_chunk_lists[obs_idx][ch_idx2].λ)-λ_max_targets[ch_idx]
            delta_λmin = minimum(time_series_of_chunk_lists[obs_idx][ch_idx2].λ)-λ_min_targets[ch_idx]
            if (abs(delta_λmax)/λ_max_targets[ch_idx]*speed_of_light_mps < max_Δv_match_chunk) &&
               (abs(delta_λmin)/λ_max_targets[ch_idx]*speed_of_light_mps < max_Δv_match_chunk)
                  pixels_in_chunk = length(time_series_of_chunk_lists[obs_idx][ch_idx2].λ)
                  if pixels_in_chunk < min_pixels_in_chunk_in_any_obs[ch_idx]
                     min_pixels_in_chunk_in_any_obs[ch_idx] = pixels_in_chunk
                  end
                  if pixels_in_chunk < min_pixels_in_chunk
                     similarity_chunk_idx[ch_idx,obs_idx] = 0
                  end
            else
                similarity_chunk_idx[ch_idx,obs_idx] = 0
            end

        end # for over ch_idx that is target
        #=
        ch_perm_idx = zeros(Int64,min_num_chunks)
        if verbose    println("# obs_idx = ", obs_idx)   end
        if verbose    println("# ch_perm_idx = ",ch_perm_idx)   end
        ch_perm_idx_non_zero = findall(x->x!=0,ch_perm_idx)
        println("# obs_idx = ", obs_idx, "  length(ch_perm_idx_non_zero) = ", length(ch_perm_idx_non_zero))
        @assert length(ch_perm_idx_non_zero) >= 1
        time_series_of_chunk_lists_output[obs_idx] = ChunkList( map(ch_idx-> time_series_of_chunk_lists[obs_idx].data[ch_idx],  ch_perm_idx[ch_perm_idx_non_zero]),
                                                             map(ch_idx-> time_series_of_chunk_lists[obs_idx].order[ch_idx], ch_perm_idx[ch_perm_idx_non_zero]) )
        =#

    end # for over obs_idx
    #num_chunks_to_keep = sum(map(ch_idx->all(similarity_chunk_idx[ch_idx,:] .!= 0), 1:min_num_chunks))
    for ch_idx in 1:min_num_chunks
        if any(similarity_chunk_idx[ch_idx,:] .== 0)
            similarity_chunk_idx[ch_idx,:] .= 0
        end
    end
    time_series_of_chunk_lists_to_keep = map(obs_idx -> ChunkList( map(keep_idx-> time_series_of_chunk_lists[obs_idx].data[similarity_chunk_idx[keep_idx,obs_idx]],  findall(x->x!=0,similarity_chunk_idx[:,obs_idx]) ),
                                                                   map(keep_idx-> time_series_of_chunk_lists[obs_idx].order[similarity_chunk_idx[keep_idx,obs_idx]], findall(x->x!=0,similarity_chunk_idx[:,obs_idx]) ) ), 1:length(time_series_of_chunk_lists) )

    chunk_list_timeseries = ChunkListTimeseries(times, time_series_of_chunk_lists_to_keep, inst=inst, metadata=metadata )
end

function make_order_pixel_list_for_chunks_in_order(spectra::AS, inst::AbstractInstrument2D, df_λ_good::DataFrame;
            orders_to_use = orders_to_use_default(spectra.inst), verbose::Bool = false
            ) where {AS<:AbstractSpectra }
    @assert eltype(orders_to_use) <: Integer
    @assert all( min_order(spectra.inst) .<= orders_to_use .<= max_order(spectra.inst) )
    @assert hasproperty(df_λ_good, :lambda_lo)
    @assert hasproperty(df_λ_good, :lambda_hi)
    min_pixels_per_chunk = 5  # arbitrary for now

    df_out = DataFrame(:pixels=>UnitRange[], :order=>Int[])
    for order_idx in orders_to_use
        pixels_to_use = min_col_default(spectra.inst,order_idx):max_col_default(spectra.inst,order_idx)
        #min_col = max(first(get_pixel_range(spectra.inst,order_idx)), 1825)
        #max_col = min(last(get_pixel_range(spectra.inst,order_idx)), 7000)
        #pixels_to_use = min_col:max_col
        if length(pixels_to_use) < min_pixels_per_chunk
            if verbose   println("# Skipping order_idx= ", order_idx, " due to pixel length")   end
            continue
        end
        if verbose    println("# Order idx= ",order_idx)   end
        # TODO OPT: could speed up by scanning through wavelengths and list together) or avoiding this allocation all together
        idx_good = is_in_wavelength_range_list.(spectra.λ[pixels_to_use,order_idx], list=df_λ_good)

        idx_start = findfirst(idx_good)
        idx_stop = findlast(idx_good)
        if isnothing(idx_start) || isnothing(idx_stop)
            if verbose   println("# Skipping order_idx= ", order_idx, " due to finding no good pixels")   end
            continue
        end
        idx_start = idx_start
        idx_stop = idx_stop
        if verbose   println("# Searching pixels:", pixels_to_use[idx_start]:pixels_to_use[idx_stop], " or λ= ",spectra.λ[pixels_to_use[idx_start],order_idx], " - ", spectra.λ[pixels_to_use[idx_stop],order_idx])   end
        idx_lo = idx_start
        idx_hi = idx_start
        while idx_hi < idx_stop
            idx_lo = findfirst( view(idx_good,idx_hi:idx_stop) )
            if isnothing(idx_lo)
                if verbose   println("# Searching pixels:", pixels_to_use[idx_hi]:pixels_to_use[idx_stop], " or λ= ",spectra.λ[pixels_to_use[idx_start],order_idx], " - ", spectra.λ[pixels_to_use[idx_stop],order_idx])   end
                break
            end
            idx_lo += idx_hi-1
            idx_hi = findfirst( view(.!idx_good,idx_lo:idx_stop) )
            if isnothing(idx_hi)
                idx_hi = idx_stop
            else
                idx_hi += idx_lo-2
            end
            if idx_hi-idx_lo > min_pixels_per_chunk # TODO Make min_pixels_per_chunk
                if verbose   println("# Found pixels:", pixels_to_use[idx_lo]:pixels_to_use[idx_hi], " or λ= ",spectra.λ[pixels_to_use[idx_lo],order_idx], " - ", spectra.λ[pixels_to_use[idx_hi],order_idx])   end
                push!(df_out,Dict(:pixels=>pixels_to_use[idx_lo]:pixels_to_use[idx_hi], :order=>order_idx))
                idx_hi += 1
            else
                idx_hi += 1
            end
        end # while still in order
    end  # for over order_to_use
    return df_out
end

""" `make_orders_into_chunks`
Return a ChunkList with a region of spectrum from each order in orders_to_use.
# Arguments
- spectra<:AbstractSpectra
- inst:  Instrument trait that provides default values
# Optional arguments
- orders_to_use: Range or Array (orders_to_use(inst))
- pixels_to_use: Array of Ranges (each from min_col to max_col)
or
- min_col: (min_col_default(inst,order)) and
- max_col: (max_col_default(inst,order))
"""
function make_orders_into_chunks
end

function make_orders_into_chunks(spectra::AS, inst::AbstractInstrument2D;
        orders_to_use = orders_to_use_default(spectra.inst) # 1:size(spectra.flux,2),
        ) where {AS<:AbstractSpectra }
    @assert eltype(orders_to_use) <: Integer
    @assert all( min_order(spectra.inst) .<= orders_to_use .<= max_order(spectra.inst) )
    pixels_to_use = map(ord->min_col_default(spectra.inst,ord):max_col_default(spectra.inst,ord),orders_to_use)
    make_orders_into_chunks(spectra,orders_to_use=orders_to_use, pixels_to_use=pixels_to_use)
end

function make_orders_into_chunks(spectra::AS;
        orders_to_use::Union{AR,AA1}, pixels_to_use::AAR ) where {
         AS<:AbstractSpectra2D, AR<:AbstractRange{Int64}, AA1<:AbstractArray{Int64,1}, AAR<:AbstractArray{AR,1} }
    @assert eltype(orders_to_use) <: Integer
    #@assert all( min_order(inst) .<= orders_to_use .<= max_order(inst) )
    #@assert minimum(pixels_to_use) >= min_pixel_in_order(inst)
    #@assert maximum(pixels_to_use) <= max_pixel_in_order(inst)
    ChunkList( map(order_idx->
                    ChunkOfSpectrum(spectra,(pixels=pixels_to_use[order_idx],order=orders_to_use[order_idx])),
                    1:length(orders_to_use) ), orders_to_use )
end

#=
# Doesn't really make sense for 1D spectra
function make_orders_into_chunks(spectra::AS;
        orders_to_use::Integer, pixels_to_use::AR ) where {
         AS<:AbstractSpectra, AR<:AbstractRange{Int64} }
    @assert eltype(orders_to_use) <: Integer
    #@assert all( min_order(inst) .<= orders_to_use .<= max_order(inst) )
    #@assert minimum(pixels_to_use) >= min_pixel_in_order(inst)
    #@assert maximum(pixels_to_use) <= max_pixel_in_order(inst)
    ChunkList( [ChunkOfSpectrum(spectra,(pixels=pixels_to_use,order=orders_to_use))], [orders_to_use] )
end
=#

""" `make_grid_for_chunk`
Create a range with equal spacing between points with end points set based on union of all chunks in timeseries.
# Arguments:
- timeseries: ChunkListTimeseries
- chunk index:
- oversample_factor: (1)
"""
function make_grid_for_chunk(timeseries::ACLT, c::Integer; oversample_factor::Real = 1.0, spacing::Symbol = :Linear, remove_rv_est::Bool = false ) where { ACLT<:AbstractChunkListTimeseries }
    num_obs = length(timeseries.chunk_list)
    @assert num_obs >= 1
    @assert 1<= c <= length(first(timeseries.chunk_list).data)
    @assert allequal(map(chunk->length(chunk.data),timeseries.chunk_list))
    @assert spacing == :Log || spacing == :Linear
    if spacing == :Log
        @warn "There's some issues with end points exceeding the bounds.  Round off error?  May cause bounds errors."
    end
    # Create grid, so that chunks at all times include the grid's minimum and maximum wavelength.
    if remove_rv_est   @assert haskey(first(timeseries.metadata),:rv_est)   end
    boost_factor = [ remove_rv_est ? calc_doppler_factor(timeseries.metadata[t][:rv_est]) : 1 for t in 1:num_obs ]
    λ_min = maximum(minimum(timeseries.chunk_list[t].data[c].λ)/boost_factor[t] for t in 1:num_obs)
    λ_max = minimum(maximum(timeseries.chunk_list[t].data[c].λ)/boost_factor[t] for t in 1:num_obs)
    Δλ_grid_obs = mean(log(timeseries.chunk_list[t].data[c].λ[end]/
                           timeseries.chunk_list[t].data[c].λ[1]   )/
                         (length(timeseries.chunk_list[t].data[c].λ)-1) for t in 1:num_obs)
    num_pixels_obs = log(λ_max/λ_min)/Δλ_grid_obs
    num_pixels_gen = (num_pixels_obs-1) * oversample_factor + 1
    if spacing == :Log
        Δlnλ_grid_obs = mean(log(timeseries.chunk_list[t].data[c].λ[end]/
                                 timeseries.chunk_list[t].data[c].λ[1]   )/
                                 (length(timeseries.chunk_list[t].data[c].λ)-1) for t in 1:num_obs)
        num_pixels_obs = log(λ_max/λ_min)/Δlnλ_grid_obs
        num_pixels_gen = (num_pixels_obs-1) * oversample_factor + 1
            Δlnλ_grid_gen = log(λ_max/λ_min)/ (num_pixels_gen-1)
        return exp.(range(log(λ_min),stop=log(λ_max),step=Δlnλ_grid_gen))
    elseif spacing == :Linear
        Δλ_grid_obs = mean((timeseries.chunk_list[t].data[c].λ[end]-timeseries.chunk_list[t].data[c].λ[1] )/
                            (length(timeseries.chunk_list[t].data[c].λ)-1) for t in 1:num_obs)
        num_pixels_obs = (λ_max-λ_min)/Δλ_grid_obs
        num_pixels_gen = (num_pixels_obs-1) * oversample_factor + 1
        Δλ_grid_gen = (λ_max-λ_min)/ (num_pixels_gen-1)
        return range(λ_min,stop=λ_max,step=Δλ_grid_gen)
    end
end

""" `make_chunk_list_around_lines( spectra, line_list)`
Return a ChunkList of best regions of spectrum with lines in line_line.
    line_list is a DataFrame containing :lambda_lo and :lambda_hi.
    Pads edges by Δ.
"""
function make_chunk_list_around_lines(spectra::AS, line_list::DataFrame; rv_shift::Real = 0, Δ::Real=Δλoλ_edge_pad_default) where { AS<:AbstractSpectra }
    @assert hasproperty(line_list,:lambda_lo)
    @assert hasproperty(line_list,:lambda_hi)
    boost_factor = calc_doppler_factor(rv_shift)
    if rv_shift != 0
        @warn("I haven't tested this yet, especially the sign.")  # TODO
    end
    line_locs = map(row->find_line_best(row.lambda_lo*boost_factor,row.lambda_hi*boost_factor,spectra,Δ=Δ), eachrow(line_list) )
    line_locs = line_locs[.!ismissing.(line_locs)]
    # TODO: Update to add order to ChunkList
    cl = ChunkList(map(loc->ChunkOfSpectrum(spectra,loc),line_locs), map(loc->loc[2], line_locs) )
end

function make_chunk_list_timeseries_around_lines(spectra::AS,line_list_df::DataFrame; rv_shift::Real = 0) where {ST<:AbstractSpectra, AS<:AbstractArray{ST,1} }
    times = map(s->s.metadata[:bjd],spectra)
    metadata = make_vec_metadata_from_spectral_timeseries(spectra)
    time_series_of_chunk_lists = map(spec->RvSpectMLBase.make_chunk_list_around_lines(spec,line_list_df, rv_shift=rv_shift), spectra)
    chunk_list_timeseries = ChunkListTimeseries(times, time_series_of_chunk_lists, inst=first(spectra).inst, metadata=metadata )
end

function extract_chunk_list_timeseries_for_order(clt::AbstractChunkListTimeseries, order::Integer) # where {ST<:AbstractSpectra, AS<:AbstractArray{ST,1} }
    # TODO: Decide/document if using order or order_idx for selecting chunks
    @assert min_order(clt.inst) <= order <= max_order(clt.inst)
    num_chunks_to_search = num_chunks(clt)
    chunk_order = map(ch->first(clt.chunk_list).data[ch].λ.indices[2], 1:num_chunks_to_search )
    chunks_in_order = findall(chunk_order .== order)
    num_chunks_in_order = sum(chunks_in_order)
    if !(num_chunks_in_order >= 1)
        return nothing
    end
    new_chunk_list_timeseries = Vector{ChunkList}(undef, length(clt.times))
    for i in 1:length(clt.times)
        chunk_list = ChunkList(clt.chunk_list[i].data[chunks_in_order], fill(order,num_chunks_in_order) )
        new_chunk_list_timeseries[i] = chunk_list
    end
    output = ChunkListTimeseries(clt.times, new_chunk_list_timeseries, inst=clt.inst, metadata=clt.metadata )
    return output
end

#=
function make_order_list_timeseries(spectra::AS) #= , order_list::AOL ) =# where {ST<:AbstractSpectra, AS<:AbstractArray{ST,1} #=, CLT<:AbstractChunkList, AOL::AbstractArray{CLT,1} =# }
    times = map(s->s.metadata[:bjd],spectra)
    inst = first(spectra).inst
    metadata = make_vec_metadata_from_spectral_timeseries(spectra)
    order_list = map( spec->RvSpectMLBase.make_orders_into_chunks(spec,inst), spectra)
    chunk_list_timeseries = ChunkListTimeseries(times, order_list, inst=first(spectra).inst, metadata=metadata )
end
=#

function make_order_list_timeseries(spectra::AS; orders_to_use = orders_to_use_default(first(spectra).inst) ) #= , order_list::AOL ) =# where {ST<:AbstractSpectra, AS<:AbstractArray{ST,1} #=, CLT<:AbstractChunkList, AOL::AbstractArray{CLT,1} =# }
    # TODO: Decide/document if using order or order_idx for selecting chunks
    times = map(s->s.metadata[:bjd],spectra)
    inst = first(spectra).inst
    metadata = make_vec_metadata_from_spectral_timeseries(spectra)
    order_list = map( spec->RvSpectMLBase.make_orders_into_chunks(spec,inst, orders_to_use=orders_to_use), spectra)
    chunk_list_timeseries = ChunkListTimeseries(times, order_list, inst=first(spectra).inst, metadata=metadata )
end

function make_chunk_list_expr(spectra::AS, range_list::DataFrame ) where { AS<:AbstractSpectra }
    # TODO: Decide/document if using order or order_idx for selecting chunks
    @assert hasproperty(range_list,:lambda_lo)
    @assert hasproperty(range_list,:lambda_hi)
    cl = typeof(ChunkOfSpectrum(spectra,1,1:size(spectra.flux,1)))[]
    ord_list = Int64[]
    for row in eachrow(range_list)
        orders = find_orders_in_range(row.lambda_lo,row.lambda_hi, spectra.λ)
        if length(orders) == 1
            order = orders[1]
            pixels = find_cols_to_fit(spectra.λ[:,order],row.lambda_lo,row.lambda_hi)
            if 1 <= first(pixels) <= size(spectra.λ,1) && 1 <= last(pixels) <= size(spectra.λ,1) &&  # valid location
                calc_snr(spectra.flux[pixels,order],spectra.var[pixels,order]) > 0     # not just NaN's
                println("# Found one order (", order, ") for range ", row.lambda_lo, " - ", row.lambda_hi )
                push!(cl, ChunkOfSpectrum(spectra,order,pixels) )
                push!(ord_list, order)
            else
                println("# Found one order (", order, ") for range ", row.lambda_lo, " - ", row.lambda_hi, " but not valid or all NaNs.")
                continue
            end
        elseif length(orders) > 1
            pixels = map(ord->find_cols_to_fit(spectra.λ[:,ord],row.lambda_lo,row.lambda_hi), orders)
            scores = map( i->calc_snr(spectra.flux[pixels[i],orders[i]],spectra.var[pixels[i],orders[i]]), 1:length(orders) )
            idx_best = findmax(scores)[2]
            println("# Found ", length(orders), " orders for range ", row.lambda_lo, " - ", row.lambda_hi, " picking ", orders[idx_best])
            push!(cl, ChunkOfSpectrum(spectra,orders[idx_best],pixels[idx_best]) )
            push!(ord_list, orders[idx_best])
            idx_all = findall(scores)

        else # length(orders) == 0
            #println("# Found no orders for range ", row.lambda_lo, " - ", row.lambda_hi )

        end
    end
    return ChunkList(cl,ord_list)
end

""" Return (chunk_timeseries) that have been trimmed of any chunks that are bad based on any spectra in the chunk_timeseries.
For now just checks for NaNs.  Instruments can provide their own checks.
Inputs:
- `chunk_timeseries`: ChunkListTimeseries
Optional arguemnts:
- `verbose`: print debugging info (false)
Returns:
- ChunkListTimeseries that has removed problematic chunks.
"""
function filter_bad_chunks(chunk_list_timeseries::ACLT; verbose::Bool = false) where { ACLT<:AbstractChunkListTimeseries }
    @assert(length(chunk_list_timeseries)>=1)
    idx_keep = trues(num_chunks(chunk_list_timeseries))
    for t in 1:length(chunk_list_timeseries)
        idx_bad_λ = findall(c->any(isnan.(chunk_list_timeseries.chunk_list[t].data[c].λ)),1:num_chunks(chunk_list_timeseries))
        idx_bad_flux = findall(c->any(isnan.(chunk_list_timeseries.chunk_list[t].data[c].flux)),1:num_chunks(chunk_list_timeseries))
        idx_bad_var = findall(c->any(isnan.(chunk_list_timeseries.chunk_list[t].data[c].var)),1:num_chunks(chunk_list_timeseries))
        idx_keep[idx_bad_λ] .= false
        idx_keep[idx_bad_flux] .= false
        idx_keep[idx_bad_var] .= false
        if verbose && (length(idx_bad_λ)+length(idx_bad_flux)+length(idx_bad_var) > 0)
               println("# Removing chunks", vcat(idx_bad_λ,idx_bad_flux,idx_bad_var), " at time ", t, " due to NaNs (",length(idx_bad_λ),",",length(idx_bad_flux),",",length(idx_bad_var),").")
        end
    end
    chunks_to_remove = findall(.!idx_keep)
    if length(chunks_to_remove) == 0
        if verbose   println("# Nothing to remove.")   end
        return chunk_list_timeseries
    else
        if verbose
            println("# Removing ", length(chunks_to_remove), " chunks due to NaNs.")
            map(c->println("# ",c),chunks_to_remove)
        end
        new_chunk_list_timeseries = [ChunkList(chunk_list_timeseries.chunk_list[t].data[idx_keep], chunk_list_timeseries.chunk_list[t].order[idx_keep]) for t in 1:length(chunk_list_timeseries) ]
        #return ChunkListTimeseries(chunk_list_timeseries.times[idx_keep],new_chunk_list_timeseries, inst=chunk_list_timeseries.inst, metadata=chunk_list_timeseries.metadata[idx_keep])
        return ChunkListTimeseries(chunk_list_timeseries.times,new_chunk_list_timeseries, inst=chunk_list_timeseries.inst, metadata=chunk_list_timeseries.metadata)
    end
end

""" Create/update metadata entry rv_est for each observation in a ChunkListTimeseries """
function set_rv_est!(clt::ACLT, rvs::AbstractVector{T})  where { ACLT<:AbstractChunkListTimeseries, T<:Real }
  @assert num_times(clt) == length(rvs)
  for i in 1:num_times(clt)
    clt.metadata[i][:rv_est] = rvs[i]
  end
  return clt
end

#= This version has some EXPRES-specific code.  Could incorporate into to EchelleInstruments.EXPRES module.
function filter_bad_chunks(chunk_list_timeseries::ACLT, line_list::DataFrame; verbose::Union{Int,Bool} = false) where { ACLT<:AbstractChunkListTimeseries }
    @assert(length(chunk_list_timeseries)>=1)
    @assert(hasproperty(line_list,:lambda_lo))
    @assert(hasproperty(line_list,:lambda_hi))
    inst = chunk_list_timeseries.inst
    idx_keep = trues(num_chunks(chunk_list_timeseries))
    for t in 1:length(chunk_list_timeseries)
        idx_bad_λ = findall(c->any(isnan.(chunk_list_timeseries.chunk_list[t].data[c].λ)),1:num_chunks(chunk_list_timeseries))
        idx_bad_flux = findall(c->any(isnan.(chunk_list_timeseries.chunk_list[t].data[c].flux)),1:num_chunks(chunk_list_timeseries))
        idx_bad_var = findall(c->any(isnan.(chunk_list_timeseries.chunk_list[t].data[c].var)),1:num_chunks(chunk_list_timeseries))
        idx_not_sorted = findall(c->!issorted(chunk_list_timeseries.chunk_list[t].data[c].λ),1:num_chunks(chunk_list_timeseries))
        idx_keep[idx_bad_λ] .= false
        idx_keep[idx_bad_flux] .= false
        idx_keep[idx_bad_var] .= false
        idx_keep[idx_not_sorted] .= false
        if verbose && (length(idx_bad_λ)+length(idx_bad_flux)+length(idx_bad_var)+length(idx_not_sorted))
               flush(stdout)
               println("# Removing chunks", vcat(idx_bad_λ,idx_bad_flux,idx_bad_var), " at time ", t, " due to NaNs (",
                            length(idx_bad_λ),",",length(idx_bad_flux),",",length(idx_bad_var),") and ",
                            length(idx_not_sorted), " for not being sorted.")
        end
        #=
        # TODO: Move these kinds of checks to a traits/plan-based system
        if hasproperty(chunk_list_timeseries.metadata[t],:pixel_mask)
            idx_keep[.!chunk_list_timeseries.metadata[t][:pixel_mask]] .= false
        end
        if hasproperty(chunk_list_timeseries.metadata[t],:excalibur_mask)
            idx_keep[.!chunk_list_timeseries.metadata[t][:excalibur_mask]] .= false
        end
        =#
    end
    chunks_to_remove = findall(.!idx_keep)
    if length(chunks_to_remove) == 0
        println("# No lines to remove.")
        return (chunk_timeseries=chunk_list_timeseries, line_list=line_list)
    else
        println("# Removing ", length(chunks_to_remove), " chunks.")
        map(c->println("# ",c,": ",line_list.lambda_lo[c]," - ",line_list.lambda_hi[c]),chunks_to_remove)
        new_line_list = line_list[findall(idx_keep),:]
        new_chunk_list_timeseries = [ChunkList(chunk_list_timeseries.chunk_list[t].data[idx_keep]) for t in 1:length(chunk_list_timeseries) ]
        return (chunk_timeseries=ChunkListTimeseries(chunk_list_timeseries.times,new_chunk_list_timeseries, inst=inst, metadata=chunk_list_timeseries.metadata), line_list=new_line_list)
    end
end
=#

""" Calc normalization of chunk based on average flux in a ChunkOfSpectrum. """
function calc_normalization(chunk::AC) where { AC<:AbstractChunkOfSpectrum}
    total_flux = NaNMath.sum(Float64.(chunk.flux))
    num_pixels = length(chunk.flux)
    scale_fac = num_pixels / total_flux
end

""" Calc normalization of chunk based on average flux in a ChunkOfSpectrum using inverse variance weighting. """
function calc_normalization_var_weighted(chunk::AC) where { AC<:AbstractChunkOfSpectrum}
    sum_weighted_flux = NaNMath.sum(Float64.(chunk.flux) ./ chunk.var )
    sum_weights = NaNMath.sum(1.0 ./ chunk.var)
    scale_fac = sum_weights / sum_weighted_flux
end

""" Calc normalization of spectra based on average flux in a ChunkList. """
function calc_normalization(chunk_list::ACL) where { ACL<:AbstractChunkList}
    #total_flux = NaNMath.sum(NaNMath.sum(Float64.(chunk_list.data[c].flux))
    #                       for c in 1:length(chunk_list) )
    total_flux = sum(sum(Float64.(chunk_list.data[c].flux))
                        for c in 1:length(chunk_list) )
    num_pixels = sum( length(chunk_list.data[c].flux) for c in 1:length(chunk_list) )
    scale_fac = num_pixels / total_flux
end

""" Calc normalization of spectra based on average flux in a ChunkList using inverse variance weighting. """
function calc_normalization_var_weighted(chunk_list::ACL) where { ACL<:AbstractChunkList}
    #total_flux = NaNMath.sum(NaNMath.sum(Float64.(chunk_list.data[c].flux))
    #                    for c in 1:length(chunk_list) )
    total_flux = sum(sum(Float64.(chunk_list.data[c].flux))
                    for c in 1:length(chunk_list) )
    num_pixels = sum( length(chunk_list.data[c].flux) for c in 1:length(chunk_list) )
    scale_fac = num_pixels / total_flux
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
# Only run coverage from linux build on travis.
get(ENV, "TRAVIS_OS_NAME", "")       == "linux"   || exit()
using Pkg
Pkg.add("Coverage")
using Coverage

cd(joinpath(dirname(@__FILE__), "..")) do
    Codecov.submit(Codecov.process_folder())
    Coveralls.submit(Coveralls.process_folder())
end

ext: jl


text:
# daq io types
abstract type AbstractIO end
abstract type AnalogIn   <: AbstractIO end
abstract type AnalogOut  <: AbstractIO end
abstract type DigitalIn  <: AbstractIO end
abstract type DigitalOut <: AbstractIO end
abstract type CounterIn  <: AbstractIO end
abstract type CounterOut <: AbstractIO end

# node for channels, perhaps unneeded...
abstract type DAQChannel end

# type aliasing for convenience
const AI = AnalogIn
const AO = AnalogOut
const DI = DigitalIn
const DO = DigitalOut
const CI = CounterIn
const CO = CounterOut

const ChannelIndex = Union{Integer,UnitRange{<:Integer}}
const AnalogIO = Union{AnalogIn, AnalogOut}

# holds non-null-terminated strings (buffers)
struct DAQStringBuffer
     str::Vector{UInt8}
    size::UInt32

    function DAQStringBuffer(size = 256)
        str  = Vector{UInt8}(undef, size)
        size = UInt32(size)
        return new(str,size)
    end
end

mutable struct DAQDevice
        name::String
    channels::Union{LittleDict,Nothing}
end

struct PhysicalChannel{T<:AbstractIO} <: DAQChannel
      name::String
    parent::DAQDevice
    ranges::Union{Vector{Tuple{Float64,Float64}},Nothing}
end

# This hasn't been implemented yet, but will act as a
# wrapper for a virtual channel that holds the underlying
# physical channel information, named alias of the virtual
# channel, and an in-memory dictionary of its attributes.
mutable struct TaskChannel{T<:AbstractIO} <: DAQChannel
     alias::String
      phys::PhysicalChannel
      attr::OrderedDict
    function TaskChannel{<:AbstractIO}(alias::String,
                                       pchan::PhysicalChannel)
        return new(alias, pchan, OrderedDict())
    end
end

mutable struct DAQTask{T<:AbstractIO}
        name::String
      handle::TaskHandle
      device::Union{DAQDevice,Nothing}
    channels::Union{Vector{PhysicalChannel{T}},Nothing}
    conditions::Vector{Base.AsyncCondition}
end


text:
module CreditApproval

export check_background, create_account, notify_downstream, open_account

# Background check.
# In practice, we would call a remote service for this.
# For this example, we just return true.
function check_background(first_name, last_name)
    println("Doing background check for $first_name $last_name")
    return true
end

# Create an account.
# In practice, we would actually create a record in database.
# For this example, we return an account number of 1.
function create_account(first_name, last_name, email)
    println("Creating an account for $first_name $last_name")
    return 1
end

# Notify downstream system by sending a message.
# For this example, we just print to console and returns nothing.
function notify_downstream(account_number)
    println("Notifying downstream system about new account $account_number")
    return nothing
end

# Open a new account.
# Returns `:success` if account is created successfully.
# Returns `:failure` if background check fails.
function open_account(first_name, last_name, email)
    check_background(first_name, last_name) || return :failure
    account_number = create_account(first_name, last_name, email)
    notify_downstream(account_number)
    return :success
end

end # module


text:
export zerossiz, onessiz, randsiz, randnsiz
export shzeros, shones, shrand, shrandn
export shzerossiz, shonessiz, shrandsiz, shrandnsiz
export zeroel, oneel
export +, *, repeat, nop, id,  istrue, isfalse, not, or, and, @dict
export plus, minus, times, divby, square, power
export any, anyequal, all, allequal, unequal
export inside
export minimum, maximum

#######################################
##  zerossiz, onessiz, randsiz, randnsiz

zerossiz(a::Array, typ::Type = Float64) = zeros(typ, a...)
onessiz(a::Array, typ::Type = Float64) = ones(typ, a...)
randsiz(a::Array, typ::Type = Float64) = rand(typ, a...)
randnsiz(a::Array) = randn(a...)
randnsiz(a::Array, typ::Type) = (r = Array(typ, a...); fillrandn(r))

shzerossiz(a::Array, typ::Type = Float64) = shzeros(typ, a...)
shonessiz(a::Array, typ::Type = Float64) = shones(typ, a...)
shrandsiz(a::Array, typ::Type = Float64) = shrand(typ, a...)
shrandnsiz(a::Array, typ::Type = Float64) = shrandn(a...)

fillrand{T}(a::SharedArray{T}) = (for i = 1:length(a) a[i] = rand(T) end; a)
fillrandn{T}(a::SharedArray{T}) = (for i = 1:length(a) a[i] = randn() end; a)
shzeros(a...) = SharedArray(Float64, a...)
shzeros(typ::Type, a...) = SharedArray(typ, a...)
shones(a...) = fill!(SharedArray(Float64, a...), one(typ))
shones(typ::Type, a...) = fill!(SharedArray(typ, a...), one(typ))
shrand(a...) = fillrand(SharedArray(Float64, a...))
shrand(typ::Type, a...) = fillrand(SharedArray(typ, a...))
shrandn(a...) = fillrandn(SharedArray(Float64, a...))
shrandn(typ::Type, a...) = fillrandn(SharedArray(typ, a...))

#######################################
## zeroel, oneel

zeroel(a) = zero(eltype(a))
oneel(a) = one(eltype(a))

macro dict(a...)
    esc(Expr(:call, :Dict, [Expr(:(=>), QuoteNode(x), x) for x in a]...))
end

import Base.*
*(a::Union{Char,AbstractString}...) = string(a...)

import Base.repeat
repeat(a::Char,n::Int) = repeat(string(a), n)
repeat(a, n) = flatten(map(unstack(1:n), x->a))

nop(a...) = return
id(a) = a
id(a...) = a

istrue(a::Bool) = a
istrue(f::Function) = istrue(f())
isfalse(a) = !(istrue(a))

#######################################
##  not, or, and

not = !
or(a,b) =  a || b
and(a,b) = a && b

#######################################
###  any, anyequal, all, allequal

import Base.any
function any(a::AbstractArray, f::Function)
    v = trytoview(a,1)
    for i = 1:len(a)
        if f(trytoview(a,i,v))
            return true
        end
    end
    return false
end

anyequal(a::AbstractArray, x) = return any(a, y -> isequal(y, x))

import Base.all
function all(a::AbstractArray, f::Function)
    v = trytoview(a,1)
    for i = 1:len(a)
        if !f(trytoview(a,i,v))
            return false
        end
    end
    return true
end

allequal(a::AbstractArray, x) = all(a, y -> isequal(y, x))

#######################################
##  internal data management

newarraysize(a::Number,n::Int) = (n,)
function newarraysize(a,n::Int)
    s = siz(a)
    if s[end] == 1
        s = s[1:end-1]
    end
    tuple(s...,n)
end

arraylike{T<:AbstractString}(a::T, n::Int, array = nothing) = Array(T, n)
function arraylike(a, n::Int, array = nothing)
    if a == nothing || (array != nothing && !(eltype(array) <: Number))
        return Array{Any}(n)
    end
    r = Array(eltype(a), newarraysize(a,n)...)
end

# function arraylike(a::Array, n::Int, array = nothing)
#     if isempty(a)
#         return Array{Any}(n)
#     end
#     r = Array(eltype(a), newarraysize(a,n)...)
# end

# arraylike(a, n::Int, array = nothing) = return Array{Any}(n)

sharraylike(a, n::Int) = SharedArray(eltype(a), newarraysize(a,n)...)

function copy(from, to)
    if length(from)!=length(to)
        error("in copy(from,to::SubArray): lengths do not match. length(from): $(length(from))  length(to): $(length(to))")
    end
    for i = 1:length(from)
        to[i] = from[i]
    end
end

plus(a,b) = a.+b
minus(a,b) = a.-b
times(a,b) = a.*b
divby(a,b) = a./b
unequal(a,b) = !(isequal(a,b))
square(a) = a.*a
power(a,n) = a.^n

inside = in

if VERSION >= v"0.4-"
    import Base: minimum, maximum
    minimum{T<:AbstractFloat}(::Type{T} = Float64) = -realmax(T)
    maximum{T<:AbstractFloat}(::Type{T} = Float64) =  realmax(T)
    minimum{T<:Number}(::Type{T}) = typemin(T)
    maximum{T<:Number}(::Type{T}) = typemax(T)
end


text:
#!/usr/bin/julia

# Julia Lang

# Written by Katrina Siegfried
# 2018-09-15
# CSCI3412_Algorithms

# This file reads in a txt file contianing an array that has the following format
# Line 1	Title of txt file
# Line 2	<number of elements>
# Line 3	n1
# Line 4	n2
# Line 5	n3
# 		...
# Line n	n

using DelimitedFiles

function read_in_file()
	# open the input file, save the first two lines as title and number of indices
	stream = open("/home/katrina/Documents/CSCI3412_Algorithms/Test.txt", "r")
	# stream = open("/home/katrina/Documents/CSCI3412_Algorithms/one-million-randoms.txt", "r")
	# stream = open("/home/katrina/Documents/CSCI3412_Algorithms/nearly-sorted.txt", "r")
	title = readline(stream, keep=false)
	numIndicesString = readline(stream, keep=false)
	numIndices = parse(Int64, numIndicesString)
	# create an empty array of length passed in the file as numIndices
	array = Array{Int64}(undef, numIndices, 1)
	index = 1
	# while not at the end of the file, read in all the remaining values as Int64
	while eof(stream) != true
		tmp = readline(stream, keep=true)
		array[index] = parse(Float64, tmp)
		index += 1
	end
	# return the array from the file for sorting
	return array
end

# CALL FUNCTION TO READ IN FILE
listToSort  = read_in_file()

# UNCOMMENT CODE BELOW TO PRINT UNSORTED LIST
for i = 1:length(listToSort)
	# println(listToSort[i])
end

# WRITE RESULTS TO NEW TXT FILE
writedlm("sortResults.txt", listToSort)

max_stars_repo_name: siegfrkn/JuliaDemoScripts

max_issues_repo_name: siegfrkn/JuliaDemoScripts

max_forks_repo_name: siegfrkn/JuliaDemoScripts


text:
module BeachLine

using ..Geometry
using ..EventQueue
using ..Diagram

@enum SIDE LEFT RIGHT

mutable struct Arc
   region::Diagram.Region
   disappearsAt::Union{EventQueue.CircleEvent, Nothing}
   parent#::Union{Breakpoint, Nothing}
   prev::Union{Arc, Nothing}
   next::Union{Arc, Nothing}
end

mutable struct Breakpoint
   leftFocus::Tuple{Real, Real}
   rightFocus::Tuple{Real, Real}
   parent::Union{Breakpoint, Nothing}
   leftChild::Union{Arc, Breakpoint}
   rightChild::Union{Arc, Breakpoint}
   halfEdge::Union{Diagram.HalfEdge, Nothing}
end

mutable struct BST
   root::Union{Arc, Breakpoint, Nothing}
end

function BST()
   return BST(nothing)
end

function fixBreakpointsUpwards(subTree)
   newBreakpoint = nothing

   leftExtreme = rightExtreme = subTree
   while isa(leftExtreme, Breakpoint)
      leftExtreme = leftExtreme.leftChild
   end
   while isa(rightExtreme, Breakpoint)
      rightExtreme = rightExtreme.rightChild
   end

   parent = subTree.parent
   while parent != nothing
      if leftExtreme == rightExtreme == nothing
         break # nothing else upwards can change
      end

      if parent.leftChild == subTree
         if rightExtreme != nothing
            if parent.leftFocus != rightExtreme.region.generator
               parent.leftFocus = rightExtreme.region.generator
               newBreakpoint = parent
            end
            rightExtreme = nothing # nothing means no change
         end
      elseif parent.rightChild == subTree
         if leftExtreme != nothing
            if parent.rightFocus != leftExtreme.region.generator
               parent.rightFocus = leftExtreme.region.generator
               newBreakpoint = parent
            end
            leftExtreme = nothing # nothing means no change
         end
      end

      subTree = subTree.parent
      parent = subTree.parent
   end

   return newBreakpoint
end

function insert(T::BST, region::Diagram.Region, ly::Real)
   arc = Arc(region, nothing, nothing, nothing, nothing)

   sideOnSpecialCase = nothing

   if T.root == nothing
      T.root = arc
      node = nothing
   else # look for the parabola immediately over arc.region.generator
      parent = nothing
      node = T.root
      side = LEFT # indicate which pointer from parent points to node

      while !isa(node, Arc)
         parent = node

         bp = Geometry.parabolaIntersection(node.leftFocus, node.rightFocus, ly)

         if arc.region.generator[1] <= bp[1]
             node = node.leftChild
             side = LEFT
         else
             node = node.rightChild
             side = RIGHT
         end
      end

      if (node.region.generator[2] == ly)
          if node.region.generator[1] <= arc.region.generator[1]
              sideOnSpecialCase = RIGHT
              newSubTree = Breakpoint(node.region.generator, arc.region.generator, parent, node, arc, nothing)
              arc.next = node.next
              if arc.next != nothing
                  arc.next.prev = arc
              end
              node.next = arc
              arc.prev = node
          else
              sideOnSpecialCase = LEFT
              newSubTree = Breakpoint(arc.region.generator, node.region.generator, parent, arc, node, nothing)
              arc.next = node
              arc.prev = node.prev
              if arc.prev != nothing
                  arc.prev.next = arc
              end
              node.prev = arc
          end

          node.parent = arc.parent = newSubTree
          node = nothing
      else
          # arrived at a leaf 'node'; switch it for subtree with tree leaves, 'arc' being the middle one, 'node' being on both the others
          newNode = Arc(node.region, nothing, nothing, arc, node.next)
          arc.prev = node
          arc.next = newNode
          if node.next != nothing
              node.next.prev = newNode
          end
          node.next = arc

          childTree = Breakpoint(arc.region.generator, node.region.generator, nothing, arc, newNode, nothing)
          newSubTree = Breakpoint(node.region.generator, arc.region.generator, parent, node, childTree, nothing)

          node.parent = newSubTree
          childTree.parent = newSubTree
          arc.parent = childTree
          newNode.parent = childTree
      end

      if parent == nothing
          T.root = newSubTree
      else
          if side == LEFT
              parent.leftChild = newSubTree
          else # RIGHT
              parent.rightChild = newSubTree
          end
      end

      fixBreakpointsUpwards(newSubTree)

      # TODO: BALANCE TREE
   end

   return arc, node, sideOnSpecialCase
end

function remove(T::BST, arc::Arc)
   if arc == T.root
      T.root = nothing

      return
   end

   parent = arc.parent

   if parent.leftChild == arc
      other = parent.rightChild
   else
      other = parent.leftChild
   end

   other.parent = parent.parent

   if parent.parent.leftChild == parent
      parent.parent.leftChild = other
   else
      parent.parent.rightChild = other
   end

   newBreakpoint = fixBreakpointsUpwards(other)

   # TODO: BALANCE TREE

   return newBreakpoint
end

# traverses tree returning only the leaves and the breakpoint x coordinates
# between them; updates positions of the extremes of the Voronoi edges
# lazy version. if we wanted speed, should be iterative
# this is also O(n), it's only used so we can draw the beachLine
# production version should not use this (see Fortune.compute())
function beachLine(T::BST, ly)
   function beachLine(node::Arc)
      return [node.region.generator]
   end

   function beachLine(node::Breakpoint)
      bp = Geometry.parabolaIntersection(node.leftFocus, node.rightFocus, ly)
      node.halfEdge.origin = bp

      return vcat(beachLine(node.leftChild), [bp], beachLine(node.rightChild))
   end

   if T.root == nothing
      return []
   else
      return beachLine(T.root)
   end
end

function printNode(node::Nothing, depth) # for debugging
   return
end

function printNode(node::Arc, depth) # for debugging
   for i in 1:depth
      print(".")
   end
   if node.parent == nothing
      print(node.region.generator, "; orphan")
   else
      print(node.region.generator, "; parent = ", node.parent.leftFocus, ",", node.parent.rightFocus)
   end
   println()
end

function printNode(node::Breakpoint, depth) # for debugging
   for i in 1:depth
      print(".")
   end
   if node.parent == nothing
      print(node.leftFocus, node.rightFocus, "; orphan")
   else
      print(node.leftFocus, node.rightFocus, "; parent = ", node.parent.leftFocus, ",", node.parent.rightFocus)
   end
   println()
   printNode(node.leftChild, depth+1)
   printNode(node.rightChild, depth+1)
end

function printTree(T::BST) # for debugging
   if T.root != nothing
      printNode(T.root, 0)
   end
end

end # module

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
module TestOptions

using ProfileVega: _options, set_options, set_default_size, Options, width, height
using Test

@testset "general options setup" begin
    set_options(width = 1000)
    @test width(_options) == 1000

    set_options(height = 600)
    @test height(_options) == 600

    set_options(negate = true)
    @test _options.negate

    set_options(simplify = false)
    @test !_options.simplify
end

@testset "size options" begin
    set_default_size(100, 200)
    @test _options.width == 100
    @test _options.height == 200
end

end # module

max_stars_repo_name: davidanthoff/ProfileVega

max_issues_repo_name: davidanthoff/ProfileVega

max_forks_repo_name: davidanthoff/ProfileVega


text:
lattice = load_lattice(joinpath(@__DIR__, "Al40.json"))
atoms   = load_atoms(  joinpath(@__DIR__, "Al40.json"))
atoms   = attach_pseudos(atoms, Al="hgh/pbe/al-q3.hgh")

cases = standard_cases()
cases["adaptive_05"] = (damping=DFTK.AdaptiveDamping(;α_trial_min=0.5, α_min=0.05), )

function load_Al40_potential(basis)
    guessfile = joinpath(@__DIR__, "Al40_initV.jld2")
    if !isfile(guessfile)
        ρ0 = guess_density(basis)
        _, ham = DFTK.energy_hamiltonian(basis, nothing, nothing; ρ=ρ0)
        V = DFTK.total_local_potential(ham)
        V .+= 3.0rand(Float64, size(V)...)

        jldopen(guessfile, "w") do file
            file["V"] = V
        end
    end

    jldopen(guessfile, "r") do file
        file["V"]
    end
end

config[@__DIR__] = (
    directory        =  @__DIR__,
    lattice          = lattice,
    atoms            = atoms,
    magnetic_moments = [],
    #
    kwargs_model     = (smearing=Smearing.Gaussian(), temperature=0.001),
    Ecut             = 20,
    kgrid            = [3, 3, 1],
    supersampling    = 1.5,
    potential        = load_Al40_potential,
    kwargs_scf       = (; ),
    mixing           = KerkerMixing(),
    cases            = cases,
)

max_stars_repo_name: mfherbst/supporting-adaptive-damping

max_issues_repo_name: mfherbst/supporting-adaptive-damping

max_forks_repo_name: mfherbst/supporting-adaptive-damping


text:
int_rules_4_3_11 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.3.11*(e*x)^m*(a+b*tan(c+d*x^n))^p =#
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (1 / ~n - 1) * (~a + ~b * tan(~c + ~d * ~x)) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && (IGtQ(1 / ~n, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~(a') + ~(b') * cot(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (1 / ~n - 1) * (~a + ~b * cot(~c + ~d * ~x)) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && (IGtQ(1 / ~n, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Unintegrable((~a + ~b * tan(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * cot(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Unintegrable((~a + ~b * cot(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(c') + ~(d') * (~u) ^ ~n)) ^ ~(p'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~a + ~b * tan(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (LinearQ(~u, ~x) && NeQ(~u, ~x))
    @apply_utils Antiderivative((~(a') + ~(b') * cot(~(c') + ~(d') * (~u) ^ ~n)) ^ ~(p'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~a + ~b * cot(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (LinearQ(~u, ~x) && NeQ(~u, ~x))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~u)) ^ ~(p'), ~x) => Antiderivative((~a + ~b * tan(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~(a') + ~(b') * cot(~u)) ^ ~(p'), ~x) => Antiderivative((~a + ~b * cot(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * tan(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * (~a + ~b * tan(~c + ~d * ~x)) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (IGtQ(Simplify((~m + 1) / ~n), 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * cot(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * (~a + ~b * cot(~c + ~d * ~x)) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (IGtQ(Simplify((~m + 1) / ~n), 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~x) ^ ~(m') * tan(~(c') + ~(d') * (~x) ^ ~n) ^ 2, ~x) => (((~x) ^ ((~m - ~n) + 1) * tan(~c + ~d * (~x) ^ ~n)) / (~d * ~n) - Antiderivative((~x) ^ ~m, ~x)) - (((~m - ~n) + 1) / (~d * ~n)) * Antiderivative((~x) ^ (~m - ~n) * tan(~c + ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~m, ~n], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * cot(~(c') + ~(d') * (~x) ^ ~n) ^ 2, ~x) => ((-((~x) ^ ((~m - ~n) + 1)) * cot(~c + ~d * (~x) ^ ~n)) / (~d * ~n) - Antiderivative((~x) ^ ~m, ~x)) + (((~m - ~n) + 1) / (~d * ~n)) * Antiderivative((~x) ^ (~m - ~n) * cot(~c + ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~m, ~n], ~x)
    #= Antiderivative((~x)^(~m')*tan((~a')+(~b')*(~x)^(~n))^(~p),~x) :=  x^(m-n+1)*tan(a+b*x^n)^(p-1)/(b*n*(p-1)) - (m-n+1)/(b*n*(p-1))*Antiderivative(x^(m-n)*tan(a+b*x^n)^(p-1),x) - Antiderivative(x^m*tan(a+b*x^n)^(p-2),x) <-- FreeQ([a,b],x) && LtQ(0,n,m+1) &&  GtQ(p,1) =#
    #= Antiderivative((~x)^(~m')*cot((~a')+(~b')*(~x)^(~n))^(~p),~x) :=  -x^(m-n+1)*cot(a+b*x^n)^(p-1)/(b*n*(p-1)) + (m-n+1)/(b*n*(p-1))*Antiderivative(x^(m-n)*cot(a+b*x^n)^(p-1),x) - Antiderivative(x^m*cot(a+b*x^n)^(p-2),x) <-- FreeQ([a,b],x) && LtQ(0,n,m+1) &&  GtQ(p,1) =#
    #= Antiderivative((~x)^(~m')*tan((~a')+(~b')*(~x)^(~n))^(~p),~x) :=  x^(m-n+1)*tan(a+b*x^n)^(p+1)/(b*n*(p+1)) - (m-n+1)/(b*n*(p+1))*Antiderivative(x^(m-n)*tan(a+b*x^n)^(p+1),x) - Antiderivative(x^m*tan(a+b*x^n)^(p+2),x) <-- FreeQ([a,b],x) && LtQ(0,n,m+1) &&  LtQ(p,-1) =#
    #= Antiderivative((~x)^(~m')*cot((~a')+(~b')*(~x)^(~n))^(~p),~x) :=  -x^(m-n+1)*cot(a+b*x^n)^(p+1)/(b*n*(p+1)) + (m-n+1)/(b*n*(p+1))*Antiderivative(x^(m-n)*cot(a+b*x^n)^(p+1),x) - Antiderivative(x^m*cot(a+b*x^n)^(p+2),x) <-- FreeQ([a,b],x) && LtQ(0,n,m+1) &&  LtQ(p,-1) =#
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * tan(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Unintegrable((~x) ^ ~m * (~a + ~b * tan(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * cot(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Unintegrable((~x) ^ ~m * (~a + ~b * cot(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~e * ~x) ^ ~(m') * (~(a') + ~(b') * tan(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * tan(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~e * ~x) ^ ~(m') * (~(a') + ~(b') * cot(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * cot(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~e * ~x) ^ ~(m') * (~(a') + ~(b') * tan(~u)) ^ ~(p'), ~x) => Antiderivative((~e * ~x) ^ ~m * (~a + ~b * tan(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~e, ~m, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~e * ~x) ^ ~(m') * (~(a') + ~(b') * cot(~u)) ^ ~(p'), ~x) => Antiderivative((~e * ~x) ^ ~m * (~a + ~b * cot(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~e, ~m, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~x) ^ ~(m') * sec(~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(p') * tan(~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(q'), ~x) => ((~x) ^ ((~m - ~n) + 1) * sec(~a + ~b * (~x) ^ ~n) ^ ~p) / (~b * ~n * ~p) - (((~m - ~n) + 1) / (~b * ~n * ~p)) * Antiderivative((~x) ^ (~m - ~n) * sec(~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (IntegerQ(~n) && (GeQ(~m, ~n) && EqQ(~q, 1)))
    @apply_utils Antiderivative((~x) ^ ~(m') * csc(~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(p') * cot(~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(q'), ~x) => (-((~x) ^ ((~m - ~n) + 1)) * csc(~a + ~b * (~x) ^ ~n) ^ ~p) / (~b * ~n * ~p) + (((~m - ~n) + 1) / (~b * ~n * ~p)) * Antiderivative((~x) ^ (~m - ~n) * csc(~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (IntegerQ(~n) && (GeQ(~m, ~n) && EqQ(~q, 1)))
end

max_stars_repo_path: src/rules/4 Trig functions/4.3 Tangent/4.3.11 (e x)^m (a+b tan(c+d x^n))^p.jl


text:
using Distances
using Optim
using DataFrames
using RCall
using Statistics
using Distributions
using Optim, NLSolversBase, Random
using LinearAlgebra: diag

function myoptimfunc(dnds, A, Δ, rλ; t = 10.0, λ = 1.0, Amin = 0.01, ρ = 100.0)
    euclidean(cumulativednds(Amin, A, ρ; t = t, r = rλ, λ = λ, Δ = Δ), dnds)
end
mutable struct dndsfit
    Δ
    rλ
    DF
    plot
end

using Bridge
#myexpint(x) = map(val -> -Bridge.expint(-val), x)
function myexpint(x)
    if length(x) == 1
        if x > 0.0
            x = -1.0
        end
        x = -Bridge.expint(-x)
        return x
    else
        ei = zeros(Float64, length(x))
        for i in 1:length(x)
            if x[i] > 0.0
                x[i] = -1.0
            end
            ei[i] = -Bridge.expint(-x[i])
        end
        return ei
    end
end

expintdiff(N, Amax, Amin, ρ) = (myexpint(-(Amax .* ρ) ./ N) .- myexpint(-(Amin .* ρ) / N))
EulerMaclaurin(N, Amax, Amin, ρ) = ((exp.(-(Amax .* ρ) ./ N) ./ (Amax .* ρ))
                .+ (exp.(-Amin .* ρ / N) ./ (Amin .* ρ))) ./ 2

function Cnintervaltheory(Amax, Amin, ρ; t = 10.0, r = 0.5, λ = 1.0, Δ = 0.0, μ = 0.001, n0 = 1)
    if Δ > 0.0
        Nave = n0 / (1 - ((r * λ * t) / (1 + r * λ * t))^n0)
        Nt = ((1 + Δ) * exp.(2 * r * λ * t * Δ) - (1 - Δ)) / (2 * Δ)
    else
        Nave = n0 / (1 - ((r * λ * t) / (1 + r * λ * t))^n0)
        Nt = 1 + r * λ * t
    end
    corr = EulerMaclaurin(Nt, Amax, Amin, ρ)
    Cnint = μ * Nave .* (1 / (1 + Δ)) .* ((myexpint(-(Amax .* ρ) ./ Nt) .- myexpint(-(Amin .* ρ) / Nt)) .+ corr)
    return Cnint
end

function cumulativednds(Amin, A, ρ; t = 10.0, r = 0.5, λ = 1.0, Δ = 0.00001)
    if Δ == 0.0
        Δ = 0.00001
    end
    Nt1 = ((1 + Δ) * exp.(2 * r * λ * t * Δ) - (1 - Δ)) / (2 * Δ)
    Nt2 = 1 + r * λ * t
    cdnds = (1 / (1 + Δ)) .* (expintdiff(Nt1, A, Amin, ρ) .+ EulerMaclaurin(Nt1, A, Amin, ρ) ) ./
    (expintdiff(Nt2, A, Amin, ρ) .+ EulerMaclaurin(Nt2, A, Amin, ρ))
    return cdnds
end

function optimizationresults(dnds, A; t = 10.0, λ = 1.0, Amin = 0.01, ρ = 100.0)
    res = optimize(x -> myoptimfunc(dnds, A, x[1], x[2]; t = t, λ = λ, Amin = Amin, ρ = ρ), [-1.0, 0.0001], [1.0, 100.0], [0.1, 0.1])
    deltares, lambdares = Optim.minimizer(res)
    DF = DataFrame(dnds = dnds, A = A,
    dndsfit = cumulativednds(Amin, A, ρ; t = t, r = lambdares, λ = λ, Δ = deltares),
    deltafit = deltares,
    lambdarfit = lambdares)

    residuals = DF[:dnds] .- DF[:dndsfit]
    SSres = sum(residuals.^2)
    SStot = sum( (DF[:dnds] .- mean(DF[:dnds])).^2)
    rsq = 1 - (SSres / SStot)
    DF[:rsq] = rsq
    @rput DF
    @rput deltares
    @rput lambdares
    x = R"""
    library(ggplot2)
    library(cowplot)
    g <- ggplot(DF, aes(x = A, y = dnds)) +
    geom_point(alpha = 0.5, col = "plum4") +
    geom_line(aes(y = dndsfit), alpha = 0.7, col = "firebrick") +
    xlab("Clone area") +
    ylab("Interval dN/dS") +
    ggtitle(paste0("Δ = ", round(deltares, 2), ", λr = ", round(lambdares, 2)))
    """
    return dndsfit(deltares, lambdares, DF, x)
end

function LLfunc(dnds, A, Δ, rλ, sigma; t = 10.0, λ = 1.0, Amin = 0.01, ρ = 100.0)
    #print(c(A, s))
    R = dnds - cumulativednds(Amin, A, ρ; t = t, r = rλ, λ = λ, Δ = Δ)
    R = sum(logpdf.(Normal(0, sigma), R))
    -sum(R)
end

function LLoptimizationresults(dnds, A; t = 10.0, λ = 1.0, Amin = 0.01, ρ = 100.0)
    function myLLfunc(x)
        LLfunc(dnds, A, x[1], x[2], x[3]; t = t, λ = λ, Amin = Amin, ρ = ρ)
    end
    res = optimize(myLLfunc, [-1.0, 0.01, 0.001], [1.0, 100.0, 1.0], [0.1, 0.1, 0.1])
    #deltares, lambdares = Optim.minimizer(res)
    params = Optim.minimizer(res)
    func = TwiceDifferentiable(myLLfunc, ones(3))
    numerical_hessian = hessian!(func, params)
    if (sum(numerical_hessian.==0.0) > 5 )
        CI = zeros(Float64, 2)
    else
        var_cov_matrix = inv(numerical_hessian)
        CI = 1.96 .* sqrt.(abs.(diag(var_cov_matrix)))
    end

    deltares = params[1]
    lambdares = params[2]

    dndsmodel = hcat(cumulativednds(Amin, A, ρ; t = t, r = lambdares - CI[2], λ = λ, Δ = deltares - CI[1]),
    cumulativednds(Amin, A, ρ; t = t, r = lambdares + CI[2], λ = λ, Δ = deltares - CI[1]),
    cumulativednds(Amin, A, ρ; t = t, r = lambdares - CI[2], λ = λ, Δ = deltares + CI[1]),
    cumulativednds(Amin, A, ρ; t = t, r = lambdares + CI[2], λ = λ, Δ = deltares + CI[1]))

    mindnds = minimum(dndsmodel, dims = 2)[:]
    maxdnds = maximum(dndsmodel, dims = 2)[:]

    #sometimes the CI calculation fails due to inability to calculate the hessian
    # this is often due to a poor fit, if this is the case we'll just set the CI to the dnds

    if sum(mindnds .< 10^-5) > 2
        @warn "Cannot compute CI intervals, returning MLE for CIs"
        mindnds = cumulativednds(Amin, A, ρ; t = t, r = lambdares, λ = λ, Δ = deltares)
    end

    if sum(maxdnds .> 10^6) > 2
        maxdnds = cumulativednds(Amin, A, ρ; t = t, r = lambdares, λ = λ, Δ = deltares)
    end

    DF = DataFrame(dnds = dnds, A = A,
    dndsfit = cumulativednds(Amin, A, ρ; t = t, r = lambdares, λ = λ, Δ = deltares),
    dndsfitlq = mindnds,
    dndsfituq = maxdnds,
    deltafit = deltares,
    lambdarfit = lambdares,
    deltafitlq = deltares - CI[1],
    lambdarfitlq = lambdares - CI[2],
    deltafituq = deltares + CI[1],
    lambdarfituq = lambdares + CI[2],
    sedelta = CI[1],
    selambda = CI[2])

    residuals = DF[:dnds] .- DF[:dndsfit]
    SSres = sum(residuals.^2)
    SStot = sum( (DF[:dnds] .- mean(DF[:dnds])).^2)
    rsq = 1 - (SSres / SStot)
    DF[:rsq] = rsq
    @rput DF
    @rput deltares
    @rput lambdares
    x = R"""
    library(ggplot2)
    library(cowplot)
    g <- ggplot(DF, aes(x = A, y = dnds)) +
    geom_point(alpha = 0.5, col = "plum4") +
    geom_line(aes(y = dndsfit), alpha = 0.8, col = "firebrick") +
    geom_ribbon(aes(ymin = dndsfitlq, ymax = dndsfituq),alpha = 0.3, fill = "firebrick") +
    xlab("Clone area") +
    ylab("Interval dN/dS") +
    ggtitle(paste0("Δ = ", round(deltares, 2), ", λr = ", round(lambdares, 2)))
    """
    return dndsfit(deltares, lambdares, DF, x)
end


text:
#=
This file implements the "agent-space interaction API" for `nothing`, ie
no space type
=#

function kill_agent!(agent::A, model::ABM{Nothing,A}) where {A<:AbstractAgent}
    delete!(model.agents, agent.id)
end

function add_agent!(agent::A, model::ABM{Nothing,A}) where {A<:AbstractAgent}
    add_agent_pos!(agent, model)
end

function add_agent_pos!(agent::A, model::ABM{Nothing,A}) where {A<:AbstractAgent}
    model[agent.id] = agent
    return model[agent.id]
end

function add_agent!(
        model::ABM{Nothing, A},
        properties...;
        kwargs...,
    ) where {A<:AbstractAgent}
    add_agent_pos!(A(nextid(model), properties...; kwargs...), model)
end


text:
"""
Cross-entropy method variants for optimization.
    - Cross-entropy method (standard)
    - Cross-entropy surrogate method
    - Cross-entropy mixture method

See paper for full description:
    http://web.stanford.edu/~mossr/pdf/cem_variants.pdf
"""
module CrossEntropyVariants

export ce_surrogate,
       ce_mixture,
       cross_entropy_method,
       sierra,
       branin,
       ackley,
       paraboloid,
       PlotSettings

using Distributions
using DataStructures
using GaussianMixtures
using GaussianProcesses
using Parameters
using Optim
using StatsBase
using Random
using Colors
using ColorSchemes
using PGFPlots
using PyPlot
using LinearAlgebra
using Suppressor

const ℝ = Float64
const ℤ = Int64

include("utils.jl")
include("sierra.jl")
include("test_objective_functions.jl")
include("surrogate_models.jl")
include("evaluation_schedules.jl")
include("plotting.jl")

dist_manhattan(𝐯, 𝐯′) = norm(𝐯 - 𝐯′, 1)
dist_euclidean(𝐯, 𝐯′) = norm(𝐯 - 𝐯′, 2)
dist_supremum(𝐯, 𝐯′) = norm(𝐯 - 𝐯′, Inf)

function nearest_neighbors(x′, 𝒟, dist)
    𝒟[argmin([dist(x, x′) for (x,y) ∈ 𝒟])][2]
end

import Distributions.MvNormal
function MvNormal(mu::Vector{Int}, Sigma::Matrix{Int})
    # Fix PDMat errors when Sigma is an integer matrix
    # Temporarty fix: https://github.com/JuliaStats/PDMats.jl/pull/117
    return MvNormal(mu, convert(Matrix{Float64}, Sigma))
end

"""
    The cross-entropy surrogate method.
"""
function ce_surrogate(f, 𝐌; # objective function `f` and initial distribution `𝐌`
                      m=10, # Number of samples per iteration
                      m_elite=5, # Number of elite samples
                      α=1, # Covariance scale for model-elites
                      k_max=10, # Maximum number of iterations
                      ϵ=1e-9, # Change in objective stopping condition (unused)
                      plot=false, # Plotting indicator
                      redraw=false, # Redraw plot every time
                      p_geo=NaN, # Evaluation scheduling parameter for geometric distribution
                      plot_settings=PlotSettings(), # Plotting keywork arguments (see `plotting`)
                      surrogate_kwargs...) # Surrogate model parameters (see `surrogate_model`)
    plot ? initial_plot(f) : nothing
    𝒟 = Queue{Vector}() # population data

    bᵥ = Inf
    bₓ = missing
    beᵥ = Inf # Best estimate value (from surrogate model).
    beₓ = missing # Best estimate x-value (from surrogate model).
    Σ = deepcopy(isa(𝐌, MixtureModel) ? 𝐌.components[1].Σ : 𝐌.Σ)

    if isnan(p_geo)
        use_evaluation_schedule = false
    else
        use_evaluation_schedule = true
        Pₛ = truncated(Geometric(p_geo), 0, k_max) # for evaluation schedule.
        Int(Random.GLOBAL_RNG.seed[1]) == 1 ? @info("Evaluation schedule: $(map(k->evaluation_schedule(Pₛ, k, k_max, m, m_elite)[1], 1:k_max))") : nothing
    end

    for k in 1:k_max
        if use_evaluation_schedule
            mₑ, m_elite = evaluation_schedule(Pₛ, k, k_max, m, m_elite)
        else
            mₑ = m
        end

        samples = rand(𝐌, mₑ)
        Y = [f(samples[:,i]) for i in 1:mₑ]
        order = sortperm(Y)
        elite = samples[:,order[1:m_elite]]

        𝐄, 𝐌, 𝒮, model_elite = model_elite_set!(𝒟, samples, elite, Y, 𝐌, Σ, m, m_elite, α, k; surrogate_kwargs...)

        plot ? plotting(k, missing, 𝐌, 𝒮, missing, m_elite, samples, model_elite, elite, true; settings=plot_settings) : nothing

        𝐌 = fit(𝐌, 𝐄)

        # Monitoring. Non-algorithmic.
        if length(order) > 1 # For the cases where evaluation_schedule returned 0
            yₜ = Y[order[1]] # Top elite y value.
            if yₜ < bᵥ
                bᵥ = yₜ # Better value.
                bₓ = samples[:,order[1]]
            end
        end

    end

    return 𝐌, bₓ, bᵥ
end

"""
Model sub-component elite set.
    ℓ: short-term memory length
"""
function model_elite_set!(𝒟, samples, elite, Y, 𝐌, Σ, m, m_elite, α, k; ℓ=3, surrogate_kwargs...)
    # Surrogate model: Gaussian process.
    𝒟ₚ = map(tuple, columns(samples), Y) # fit to entire population

    # Use short-term memory
    enqueue!(𝒟, 𝒟ₚ)
    while length(𝒟) > ℓ
        dequeue!(𝒟)
    end
    𝒟 = reduce(vcat, 𝒟) # Re-cast, abuse of variable notation
    𝒮hat = surrogate_model(𝒟; surrogate_kwargs...)

    m_model = 10m
    m_model_elite = 10m_elite
    model_samples = rand(𝐌, m_model)
    model_Y = [𝒮hat(model_samples[:,i]) for i in 1:(m_model)]
    model_order = sortperm(model_Y)
    model_elite = model_samples[:,model_order[1:(m_model_elite)]]

    # Elites form their own distributions based on the surrogate model.
    cov = Σ/α
    𝐦 = map(e->MvNormal(e, cov), columns(elite))

    𝐄 = Matrix{ℝ}(undef, 2, 0) # Elite set.

    # For each elite in 𝐦, run CE-method to completion to choose new "sub-elites" (use μ = eₓ and Σ = M.Σ)
    sub_elite = Matrix{ℝ}(undef, 2, 0) # Sub-elite set.
    for i in 1:length(𝐦)
        # Note. If this k_max is too large, then we could converge/overfit to the surrogate model
        # (especially when the outer `k` is low, the surrogate model has not matured yet)
        𝐦[i], cem_bₓ, cem_bᵥ = cross_entropy_method(𝒮hat, 𝐦[i]; m=100,
                                                                 m_elite=10,
                                                                 k_max=2)
        # Add top elite, not the mean.
        sub_elite = hcat(sub_elite, cem_bₓ)
    end

    if isa(𝐌, MixtureModel)
        # Mix 𝐦 into 𝐌 (mixture)
        𝐌 = MixtureModel(𝐦)
    end

    𝐄 = hcat(𝐄, elite, model_elite, sub_elite)
    return 𝐄, 𝐌, 𝒮hat, model_elite
end

"""
Cross-entropy mixture method. Same as `ce_surrogate` but using mixture models.
"""
function ce_mixture(f, 𝐌; kwargs...)
    𝐌 = MixtureModel([𝐌])
    return ce_surrogate(f, 𝐌; kwargs...)
end

import Distributions.fit
fit(𝐌::MixtureModel, 𝐄; kwargs...) = fit(GMM(𝐌), 𝐄; kwargs...) # Cast input 𝐌 to GMM
"""
Gaussian mixture model fitting using the Expectation Maximization algorithm.
"""
function fit(𝐌::GMM, 𝐄)
    try
        # Expectation maximization to fit mixture model 𝐌 to elite set (data) 𝐄
        @suppress em!(𝐌, permutedims(𝐄))
    catch err
        @warn err
    end

    return MixtureModel(𝐌) # re-cast
end

"""
Fit the multivariate Normal distribution 𝐌 to data set 𝐄.
"""
function fit(𝐌::MvNormal, 𝐄)
    try
        𝐌 = fit(typeof(𝐌), 𝐄)
    catch err
        @warn err
    end

    return 𝐌
end

# Kochenderfer and Wheeler, Algorithms for Optimization, Algorithm 8.7, Page 135
"""
    The cross-entropy method, which takes an objective
    function `f` to be minimized, a proposal distribution
    `P`, an iteration count `k_max`, a sample size `m`, and
    the number of samples to use when refitting the
    distribution `m_elite`. It returns the updated distribution
    over where the global minimum is likely to exist.
"""
function cross_entropy_method(f, P; m=10, m_elite=5, k_max=10, ϵ=1e-9, plot=false, redraw=false, p_geo=nothing)
    bᵥ = Inf
    bₓ = missing
    plot ? initial_plot(f) : nothing
    μ = mean(P)

    for k in 1:k_max
        samples = rand(P, m)
        Y = [f(samples[:,i]) for i in 1:m]
        order = sortperm(Y)
        elite = samples[:,order[1:m_elite]]

        # Top elite.
        yₜ = Y[order[1]] # Top elite y value
        if yₜ < bᵥ
            bᵥ = yₜ # Found better. (Monitoring only)
            bₓ = samples[:,order[1]]
        end

        plot ? plot_cem(P, samples, order, m_elite) : nothing

        try
            P = fit(typeof(P), elite)
        catch err
            @warn err
            @info k
        end

        μ = mean(P)
    end

    return (P, bₓ, bᵥ)
end

end # module CrossEntropyVariants


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 92ce9460-f62b-11ea-1a8c-179776b5a0b4
using Distributions, Parameters, POMDPStressTesting, Latexify, PlutoUI

# ╔═╡ 2978b840-f62d-11ea-2ea0-19d7857208b1
md"""
# Black-Box Stress Testing
"""

# ╔═╡ 1e00230e-f630-11ea-3e40-bf8c852f78b8
# begin
#   using Pkg
#   pkg"registry add https://github.com/JuliaPOMDP/Registry"
#   pkg"add https://github.com/sisl/RLInterface.jl"
#   pkg"add https://github.com/sisl/POMDPStressTesting.jl"
# end
md"*Unhide for installation (waiting on Julia registry).*"

# ╔═╡ 40d3b1e0-f630-11ea-2160-01338d9f2209
md"""
To find failures in a black-box autonomous system, we can use the `POMDPStressTesting` package which is part of the POMDPs.jl ecosystem.

Various solvers—which adhere to the POMDPs.jl interface—can be used:
- `MCTSPWSolver` (MCTS with action progressive widening)
- `TRPOSolver` and `PPOSolver` (deep reinforcement learning policy optimization)
- `CEMSolver` (cross-entropy method)
- `RandomSearchSolver`
"""

# ╔═╡ 86f13f60-f62d-11ea-3241-f3f1ffe37d7a
md"""
## Simple Problem: One-Dimensional Walk
We define a simple problem for adaptive stress testing (AST) to find failures. This problem, called Walk1D, samples random walking distances from a standard normal distribution $\mathcal{N}(0,1)$ and defines failures as walking past a certain threshold (which is set to $\pm 10$ in this example). AST will either select the seed which deterministically controls the sampled value from the distribution (i.e. from the transition model) or will directly sample the provided environmental distributions. These action modes are determined by the seed-action or sample-action options. AST will guide the simulation to failure events using a notion of distance to failure, while simultaneously trying to find the set of actions that maximizes the log-likelihood of the samples.
"""

# ╔═╡ d3411dd0-f62e-11ea-27d7-1b2ed8edc415
md"""
## Gray-Box Simulator and Environment
The simulator and environment are treated as gray-box because we need access to the state-transition distributions and their associated likelihoods.
"""

# ╔═╡ e37d7542-f62e-11ea-0b61-513a4b44fc3c
md"""
##### Parameters
First, we define the parameters of our simulation.
"""

# ╔═╡ fd7fc880-f62e-11ea-15ac-f5407aeff2a6
@with_kw mutable struct Walk1DParams
    startx::Float64 = 0   # Starting x-position
    threshx::Float64 = 10 # +- boundary threshold
    endtime::Int64 = 30   # Simulate end time
end;

# ╔═╡ 012c2eb0-f62f-11ea-1637-c113ad01b144
md"""
##### Simulation
Next, we define a `GrayBox.Simulation` structure.
"""

# ╔═╡ 0d7049de-f62f-11ea-3552-214fc4e7ec98
@with_kw mutable struct Walk1DSim <: GrayBox.Simulation
    params::Walk1DParams = Walk1DParams() # Parameters
    x::Float64 = 0 # Current x-position
    t::Int64 = 0 # Current time ±
    distribution::Distribution = Normal(0, 1) # Transition distribution
end;

# ╔═╡ 11e445d0-f62f-11ea-305c-495272981112
md"""
### GrayBox.environment
Then, we define our `GrayBox.Environment` distributions. When using the `ASTSampleAction`, as opposed to `ASTSeedAction`, we need to provide access to the sampleable environment.
"""

# ╔═╡ 43c8cb70-f62f-11ea-1b0d-bb04a4176730
GrayBox.environment(sim::Walk1DSim) = GrayBox.Environment(:x => sim.distribution)

# ╔═╡ 48a5e970-f62f-11ea-111d-35694f3994b4
md"""
### GrayBox.transition!
We override the `transition!` function from the `GrayBox` interface, which takes an environment sample as input. We apply the sample in our simulator, and return the log-likelihood.
"""

# ╔═╡ 5d0313c0-f62f-11ea-3d33-9ded1fb804e7
function GrayBox.transition!(sim::Walk1DSim, sample::GrayBox.EnvironmentSample)
    sim.t += 1 # Keep track of time
    sim.x += sample[:x].value # Move agent using sampled value from input
    return logpdf(sample)::Real # Summation handled by `logpdf()`
end

# ╔═╡ 6e111310-f62f-11ea-33cf-b5e943b2f088
md"""
## Black-Box System
The system under test, in this case a simple single-dimensional moving agent, is always treated as black-box. The following interface functions are overridden to minimally interact with the system, and use outputs from the system to determine failure event indications and distance metrics.
"""

# ╔═╡ 7c84df7e-f62f-11ea-3b5f-8b090654df19
md"""
### BlackBox.initialize!
Now we override the `BlackBox` interface, starting with the function that initializes the simulation object. Interface functions ending in `!` may modify the `sim` object in place.
"""

# ╔═╡ 9b736bf2-f62f-11ea-0330-69ffafe9f200
function BlackBox.initialize!(sim::Walk1DSim)
    sim.t = 0
    sim.x = sim.params.startx
end

# ╔═╡ a380e250-f62f-11ea-363d-2bf2b59d5eed
md"""
### BlackBox.distance
We define how close we are to a failure event using a non-negative distance metric.
"""

# ╔═╡ be39db60-f62f-11ea-3a5c-bd57114455ff
BlackBox.distance(sim::Walk1DSim) = max(sim.params.threshx - abs(sim.x), 0)

# ╔═╡ bf8917b0-f62f-11ea-0e77-b58065b0da3e
md"""
### BlackBox.isevent
We define an indication that a failure event occurred.
"""

# ╔═╡ c5f03110-f62f-11ea-1119-81f5c9ec9283
BlackBox.isevent(sim::Walk1DSim) = abs(sim.x) ≥ sim.params.threshx

# ╔═╡ c378ef80-f62f-11ea-176d-e96e1be7736e
md"""
### BlackBox.isterminal
Similarly, we define an indication that the simulation is in a terminal state.
"""

# ╔═╡ cb5f7cf0-f62f-11ea-34ca-5f0656eddcd4
function BlackBox.isterminal(sim::Walk1DSim)
    return BlackBox.isevent(sim) || sim.t ≥ sim.params.endtime
end

# ╔═╡ e2f34130-f62f-11ea-220b-c7fc7de2c7e7
md"""
### BlackBox.evaluate!
Lastly, we use our defined interface to evaluate the system under test. Using the input sample, we return the log-likelihood, distance to an event, and event indication.
"""

# ╔═╡ f6213a50-f62f-11ea-07c7-2dcc383c8042
function BlackBox.evaluate!(sim::Walk1DSim, sample::GrayBox.EnvironmentSample)
    logprob::Real = GrayBox.transition!(sim, sample) # Step simulation
    d::Real       = BlackBox.distance(sim) # Calculate miss distance
    event::Bool   = BlackBox.isevent(sim) # Check event indication
    return (logprob::Real, d::Real, event::Bool)
end

# ╔═╡ 01da7aa0-f630-11ea-1262-f50453455766
md"""
## AST Setup and Running
Setting up our simulation, we instantiate our simulation object and pass that to the Markov decision proccess (MDP) object of the adaptive stress testing formulation. We use Monte Carlo tree search (MCTS) with progressive widening on the action space as our solver. Hyperparameters are passed to `MCTSPWSolver`, which is a simple wrapper around the POMDPs.jl implementation of MCTS. Lastly, we solve the MDP to produce a planner. Note we are using the `ASTSampleAction`.
"""

# ╔═╡ fdf55130-f62f-11ea-33a4-a783b4d216dc
function setup_ast(seed=0)
    # Create gray-box simulation object
    sim::GrayBox.Simulation = Walk1DSim()

    # AST MDP formulation object
    mdp::ASTMDP = ASTMDP{ASTSampleAction}(sim)
    mdp.params.debug = true # record metrics
    mdp.params.top_k = 10   # record top k best trajectories
    mdp.params.seed = seed  # set RNG seed for determinism

    # Hyperparameters for MCTS-PW as the solver
    solver = MCTSPWSolver(n_iterations=1000,        # number of algorithm iterations
                          exploration_constant=1.0, # UCT exploration
                          k_action=1.0,             # action widening
                          alpha_action=0.5,         # action widening
                          depth=sim.params.endtime) # tree depth

    # Get online planner (no work done, yet)
    planner = solve(solver, mdp)

    return planner
end;

# ╔═╡ 09c928f0-f631-11ea-3ef7-512a6bececcc
md"""
### Searching for Failures
After setup, we search for failures using the planner and output the best action trace.
"""

# ╔═╡ 17d0ed20-f631-11ea-2e28-3bb9ca9a445f
planner = setup_ast();

# ╔═╡ 1c47f652-f631-11ea-15f6-1b9b59700f36
action_trace = search!(planner)

# ╔═╡ 21530220-f631-11ea-3994-319c862d51f9
md"""
### Playback
We can also playback specific trajectories and print intermediate $x$-values.
"""

# ╔═╡ 3b282ae0-f631-11ea-309d-639bf4411bb3
playback_trace = playback(planner, action_trace, sim->sim.x, return_trace=true)

# ╔═╡ 7473adb0-f631-11ea-1c87-0f76b18a9ab6
failure_rate = print_metrics(planner)

# ╔═╡ b6244db0-f63a-11ea-3b48-89d427664f5e
md"""
## Other Solvers: Cross-Entropy Method
We can easily take our `ASTMDP` object (`planner.mdp`) and re-solve the MDP using a different solver—in this case the `CEMSolver`.
"""

# ╔═╡ e15cc1b0-f63a-11ea-2401-5321d48118c3
mdp = planner.mdp; # reused from above

# ╔═╡ f5a15af0-f63a-11ea-1dd7-593d7cb01ee4
cem_solver = CEMSolver(n_iterations=1000, episode_length=mdp.sim.params.endtime)

# ╔═╡ fb3fa610-f63a-11ea-2663-17224dc8aade
cem_planner = solve(cem_solver, mdp);

# ╔═╡ 09c9e0b0-f63b-11ea-2d50-4154e3432fa0
cem_action_trace = search!(cem_planner);

# ╔═╡ 46b40e10-f63b-11ea-2375-1976bb637d51
md"Notice the failure rate is about 10x of `MCTSPWSolver`."

# ╔═╡ 32fd5cf2-f63b-11ea-263a-39f013ef6d68
cem_failure_rate = print_metrics(cem_planner)

# ╔═╡ 801f8080-f631-11ea-0728-f15dddc3ef5d
md"""
## AST Reward Function
The AST reward function gives a reward of $0$ if an event is found, a reward of negative distance if no event is found at termination, and the log-likelihood during the simulation.
"""

# ╔═╡ 8f06f650-f631-11ea-1c52-697060322173
@latexify function R(p,e,d,τ)
    if τ && e
        return 0
    elseif τ && !e
        return -d
    else
        return log(p)
    end
end

# ╔═╡ 9463f6e2-f62a-11ea-1cef-c3fa7d4f19ad
md"""
## References
1. Robert J. Moss, Ritchie Lee, Nicholas Visser, Joachim Hochwarth, James G. Lopez, and Mykel J. Kochenderfer, "Adaptive Stress Testing of Trajectory Predictions in Flight Management Systems", *Digital Avionics Systems Conference, 2020.*
"""

# ╔═╡ 05be9a80-0877-11eb-3d88-efbd306754a2
PlutoUI.TableOfContents("POMDPStressTesting.jl")

# ╔═╡ Cell order:
# ╟─2978b840-f62d-11ea-2ea0-19d7857208b1
# ╟─1e00230e-f630-11ea-3e40-bf8c852f78b8
# ╠═92ce9460-f62b-11ea-1a8c-179776b5a0b4
# ╟─40d3b1e0-f630-11ea-2160-01338d9f2209
# ╟─86f13f60-f62d-11ea-3241-f3f1ffe37d7a
# ╟─d3411dd0-f62e-11ea-27d7-1b2ed8edc415
# ╟─e37d7542-f62e-11ea-0b61-513a4b44fc3c
# ╠═fd7fc880-f62e-11ea-15ac-f5407aeff2a6
# ╟─012c2eb0-f62f-11ea-1637-c113ad01b144
# ╠═0d7049de-f62f-11ea-3552-214fc4e7ec98
# ╟─11e445d0-f62f-11ea-305c-495272981112
# ╠═43c8cb70-f62f-11ea-1b0d-bb04a4176730
# ╟─48a5e970-f62f-11ea-111d-35694f3994b4
# ╠═5d0313c0-f62f-11ea-3d33-9ded1fb804e7
# ╟─6e111310-f62f-11ea-33cf-b5e943b2f088
# ╟─7c84df7e-f62f-11ea-3b5f-8b090654df19
# ╠═9b736bf2-f62f-11ea-0330-69ffafe9f200
# ╟─a380e250-f62f-11ea-363d-2bf2b59d5eed
# ╠═be39db60-f62f-11ea-3a5c-bd57114455ff
# ╟─bf8917b0-f62f-11ea-0e77-b58065b0da3e
# ╠═c5f03110-f62f-11ea-1119-81f5c9ec9283
# ╟─c378ef80-f62f-11ea-176d-e96e1be7736e
# ╠═cb5f7cf0-f62f-11ea-34ca-5f0656eddcd4
# ╟─e2f34130-f62f-11ea-220b-c7fc7de2c7e7
# ╠═f6213a50-f62f-11ea-07c7-2dcc383c8042
# ╟─01da7aa0-f630-11ea-1262-f50453455766
# ╠═fdf55130-f62f-11ea-33a4-a783b4d216dc
# ╟─09c928f0-f631-11ea-3ef7-512a6bececcc
# ╠═17d0ed20-f631-11ea-2e28-3bb9ca9a445f
# ╟─1c47f652-f631-11ea-15f6-1b9b59700f36
# ╟─21530220-f631-11ea-3994-319c862d51f9
# ╠═3b282ae0-f631-11ea-309d-639bf4411bb3
# ╠═7473adb0-f631-11ea-1c87-0f76b18a9ab6
# ╟─b6244db0-f63a-11ea-3b48-89d427664f5e
# ╠═e15cc1b0-f63a-11ea-2401-5321d48118c3
# ╠═f5a15af0-f63a-11ea-1dd7-593d7cb01ee4
# ╠═fb3fa610-f63a-11ea-2663-17224dc8aade
# ╠═09c9e0b0-f63b-11ea-2d50-4154e3432fa0
# ╟─46b40e10-f63b-11ea-2375-1976bb637d51
# ╠═32fd5cf2-f63b-11ea-263a-39f013ef6d68
# ╟─801f8080-f631-11ea-0728-f15dddc3ef5d
# ╠═8f06f650-f631-11ea-1c52-697060322173
# ╟─9463f6e2-f62a-11ea-1cef-c3fa7d4f19ad
# ╠═05be9a80-0877-11eb-3d88-efbd306754a2


text:
module DenseMapPosterior

include("RGBD_Fusion.jl")
include("DMP_RGBD.jl")

end # module


text:
# GlaThiDa ice thickness loader
#
# This is a bit of a mess as it is hard to find the right glacier.
# Strategy is to just use Matthias' compiled data.

"""
Returns thickness data if there is any.  Uses Matthias's compiled data instead of directly GlaThiDa.
"""
function load!(ds::DataSet{ThicknessData, GlaThiDaLoader}, dem,
               glaciermask, target_proj, outline::Outline{RGIGlacier}, cropbox,
               pp::Phys)
    if isfile(ds.files[1])
        # get radar lines
        out, head = readdlm(ds.files[1], F, header=true)
        x = out[:,2]
        y = out[:,3]
        z = out[:,4]
        z[z.==-9999] = FILL
        bed = out[:,5]
        bed[bed.==-9999] = FILL
        thick = out[:,6]
        thick[thick.==-9999] = FILL
        if all(z.==FILL) && all(thick.!=FILL) && all(bed.!=FILL)
            z .= thick.+bed
        elseif all(thick.==FILL)
            thick .= z.-bed
        end
        keep = Int[]
        for (i,t) in enumerate(thick)
            t!=FILL && push!(keep,i)
        end

        # coordinate transform
        for (i,(xx,yy,zz)) in enumerate(zip(x,y,z))
            x[i],y[i],z[i] = transform_proj([xx,yy,zz], "+proj=longlat +datum=WGS84" , target_proj)
        end
        tr = Traj(x=x[keep],y=y[keep],v=thick[keep],err=z[keep])

        # some radar data has all zero, if so discard
        if all(abs.(tr.v).<1)
            println("All thickness data less than 1m, not using GlaThiDa measurements.")
            return SomeData(nothing)
        end

        # Sort it.  Do this a few times, makes things better
        if length(tr)>3
            for i = 1:3
                if !VAWTools.issorted_traj(tr, 1:min(length(tr), 15)) && length(tr)<2*10^4
                    tr = VAWTools.sort_traj(tr)
                end

                # purge points closer than some distance apart
                len = sqrt(outline.gid.attrs.area)
                window = len/100
                tr = VAWTools.downsample(tr, window)
            end
        end

        SomeData(vals=tr, sigma=ds.opts[:sigma])
    elseif isfile(ds.files[2])
        return SomeData(nothing)

        # # Get a mean thickness in AVG_*.dat
        # out, head = readdlm(ds.files[2], F, header=true)
        # rginr = out[:,3]
        # @assert length(unique(rginr))==length(rginr)
        # i = findfirst(rginr.==getrgi_nr(outline.gid))
        # if i==0
        #     return SomeData(nothing)
        # else
        #     meanh = out[i,end]
        #     return SomeData(vals=meanh, sigma=ds.opts[:sigma_mean_h_rel]*meanh)
        # end
    else
        return SomeData(nothing)
    end
end

# ##
# ttt = CSV.read(GZip.open("/home/mauro/julia/mypkgs/BITEModel/scripts/data/GlaThiDa_2016/TTT_2016_corr.csv.gz"), datarow=5,  rows_for_type_detect=10^6, header=3)
# pu = by(ttt, :GlaThiDa_ID, d->DataFrame(c=first(d[:POLITICAL_UNIT]),
#                                    lat=mean(d[:POINT_LAT]),
#                                    lon=mean(d[:POINT_LON])
#                                    ))
# pu[pu[:c].=="NO",:] # two in Svalbard

# t = CSV.read("/home/mauro/julia/mypkgs/BITEModel/scripts/data/GlaThiDa_2016/T_2016_corr.csv", datarow=3,  rows_for_type_detect=10^6, header=1)
# no = t[t[:POLITICAL_UNIT].=="NO",:]
# no = no[no[:LAT].>75,:]


text:
# test_suite.jl

using SignalExtensions
using Base.Test
using OffsetArrays

include("test_util.jl")
include("test_sequences.jl")
include("test_extensions.jl")
include("test_lazy.jl")

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

function run_tests()
    delimit("Utility functions")
    test_util()

    delimit("Sequences")
    test_sequences()

    delimit("Extensions")
    test_extensions()

    delimit("Lazy sequences")
    test_lazy()
end

run_tests()


text:
# filehandling.jl (v0.1)

import JSON

export readfile
export writefile

function readfile(fname)::Array{Any, 1}
	open(fname, "r") do f
		return JSON.parse(read(f, String));
	end
end

function writefile(fname, content)
	cp(fname, fname * ".backup"; force = true);

	try
		rm(fname);
		open(fname, "w") do f
			write(f, JSON.json(content))
		end
	catch
		println("Removal and rewriting of " * fname * " was unsuccessful.");
		println("Restoring " * fname * ".backup");
	end

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to the API function to invert rotations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/inv_rotations.jl
# ============================

# Functions: inv_rotation
# -----------------------

@testset "Invert rotations" begin
    for T in (Float32, Float64)
        # DCM
        # ======================================================================

        # Create a random DCM.
        D = rand(DCM{T})

        Di = inv_rotation(D)
        @test eltype(Di) === T

        Die = inv(D)
        @test Di ≈ Die

        # Quaternion
        # ======================================================================

        # Create a random quaternion.
        q = rand(Quaternion{T})

        qi = inv_rotation(q)
        @test eltype(qi) === T

        qie = inv(q)
        @test qi ≈ qie

        # Euler angle and axis
        # ======================================================================

        # Create a random Euler angle and axis.
        av = rand(EulerAngleAxis{T})

        avi = inv_rotation(av)
        @test eltype(avi) === T

        avie = inv(av)
        @test avi ≈ avie

        # Euler angles
        # ======================================================================

        # Create random Euler angles.
        ea = rand(EulerAngles{T})

        eai = inv_rotation(ea)
        @test eltype(eai) === T

        eaie = inv(ea)
        @test eai ≈ eaie
    end
end


text:
# This file was generated, do not modify it. # hide
using ScientificTypes
data = coerce(boston, :Tax=>Continuous, :Rad=>Continuous);

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
using IncrementalInference
using Test

##

@testset "Only do UPWARD_COMMON message likelihoods if hasPriors" begin

##

fg = generateGraph_LineStep(4;
                            poseEvery=1,
                            landmarkEvery=5,
                            posePriorsAt=[],
                            landmarkPriorsAt=[],
                            sightDistance=5,
                            solverParams=SolverParams(algorithms=[:default, :parametric]))

deleteFactor!.(fg, [Symbol("x$(i)lm0f1") for i=1:3])

#force wrong init
prpo = Prior(Normal(5, 0.01))
addFactor!(fg, [:x0], prpo)
initAll!(fg)
deleteFactor!(fg, :x0f1)

# now the correct prior
prpo = Prior(Normal(0, 0.01))
addFactor!(fg, [:x0], prpo)

fg.solverParams.useMsgLikelihoods = true

smtasks = Task[]
tree = solveTree!(fg; smtasks=smtasks, verbose=true, timeout=30);

@warn("hasPriors test needs multiple solves")
tree = solveTree!(fg);
tree = solveTree!(fg);
# tree = solveTree!(fg; smtasks, verbose=true, timeout=20, recordcliqs=ls(fg));

for i = 0:4
  ppe = getPPE(getVariable(fg, Symbol("x$i"))).suggested[1]
  @show i ppe
  @test isapprox(ppe, i; atol=0.7)
end

##

end

#


text:
"""
Represents an input to a program implemented by a Koza tree.
"""
type KozaInput
  """
  The label, or name of this input.
  """
  label::AbstractString

  """
  The (name of the) type used by this argument.
  """
  ty::AbstractString

  KozaInput(def::AbstractString) = new(def[1:Base.search(def, ':')-1],
    def[Base.last(Base.search(def, "::"))+1:end])
end

"""
Compiles an ordered list of koza inputs into a partial function header
signature.
"""
compile(il::Vector{KozaInput}) =
  join(["$(i.label)::$(i.ty)" for i in il], ",")


text:
#using FactCheck
facts("Top-level interface") do
  context("run a simple optimization with mostly defaults") do
    res = bboptimize(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 2,
      MaxTime = 0.5, TraceMode = :silent)
    @fact bestfitness(res) < 0.01 => true
    xbest = best(res)
    @fact typeof(xbest) => Vector{Float64}
    @fact length(xbest) => 2
  end

  context("continue running an optimization after it finished") do
    optctrl = bbsetup(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 2,
      MaxTime = 0.5, TraceMode = :silent)
    res1 = bboptimize(optctrl)
    @fact numruns(optctrl) => 1

    res2 = bboptimize(optctrl; MaxTime = 1.0)
    @fact numruns(optctrl) => 2

    @fact best_fitness(res2) <= best_fitness(res1) => true

    # parameters should be the same except for MaxTime
    for p in keys(parameters(res1))
      if p != :MaxTime
        @fact parameters(res1)[p] == parameters(res2)[p] => true
      end
    end
    @fact parameters(res1)[:MaxTime] => 0.5
    @fact parameters(res2)[:MaxTime] => 1.0
  end

  context("continue running an optimization after serializing to disc") do
    optctrl = bbsetup(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 2,
      MaxTime = 0.5, TraceMode = :silent)
    res1 = bboptimize(optctrl)

    tempfilename = "./temp" * string(rand(1:int(1e8))) * ".tmp"
    open(tempfilename, "w") do fh
      serialize(fh, optctrl)
    end

    # Try to make sure its not in mem:
    opctrl = nothing; gc()

    open(tempfilename, "w") do fh
      ocloaded = deserialize(fh)
    end

    @fact numruns(ocloaded) => 1
    res2 = bboptimize(ocloaded; MaxTime = 1.0)
    @fact numruns(ocloaded) => 2

    @fact bestfitness(res2) <= bestfitness(res1) => true
  end

  context("custom user termination decision") do
    # You can implement your own termination condition by inserting a callback function that returns
    # an OptimizationAction (such as, for example, StopOptimization) or just nothing for no action.
    terminate_func = (optstate) -> numsteps(optstate) >= 100 ? BlackBoxOptim.StopOptimization : nothing
    res = bboptimize(rosenbrock; MaxSteps = 500, SearchRange = (-5.0, 5.0), NumDimensions = 2,
      TraceMode = :silent, StepCallback => terminate_func)
    @fact numsteps(res) => 100
  end

  context("accessing the trace of the optimization") do
    res = bboptimize(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 2,
      MaxTime = 0.5, TraceMode = :silent)
    tr = trace(res)
    numimprovs = length(tr)
    # Last fitness value in trace is the best fitness achieved
    @fact bestfitness(res) == tr[numimprovs] => true
    if numimprovs > 1
      @fact tr[numimprovs] < tr[numimprovs-1] => true
    end
  end
end


text:
using Dates: parse_components, default_format

@testset "parse" begin
    @test isequal(
        parse(ZonedDateTime, "2017-11-14 11:03:53 +0100", dateformat"yyyy-mm-dd HH:MM:SS zzzzz"),
        ZonedDateTime(2017, 11, 14, 11, 3, 53, tz"UTC+01"),
    )
    @test isequal(
        parse(ZonedDateTime, "2016-04-11 08:00 UTC", dateformat"yyyy-mm-dd HH:MM ZZZ"),
        ZonedDateTime(2016, 4, 11, 8, tz"UTC"),
    )
    # two-digit time zone
    @test isequal(
        parse(ZonedDateTime, "2000+00", dateformat"yyyyz"),
        ZonedDateTime(2000, tz"UTC"),
    )
    @test_throws ArgumentError parse(ZonedDateTime, "2016-04-11 08:00 EST", dateformat"yyyy-mm-dd HH:MM zzz")
    # test AbstractString
    @test isequal(
        parse(ZonedDateTime, Test.GenericString("2018-01-01 00:00 UTC"), dateformat"yyyy-mm-dd HH:MM ZZZ"),
        ZonedDateTime(2018, 1, 1, 0, tz"UTC"),
    )

end

@testset "tryparse" begin
    @test isequal(
        tryparse(ZonedDateTime, "2013-03-20 11:00:00+04:00", dateformat"y-m-d H:M:SSz"),
        ZonedDateTime(2013, 3, 20, 11, tz"UTC+04"),
    )
    @test isequal(
        tryparse(ZonedDateTime, "2016-04-11 08:00 EST", dateformat"yyyy-mm-dd HH:MM zzz"),
        nothing,
    )
end

@testset "parse components" begin
    local test = ("2017-11-14 11:03:53 +0100", dateformat"yyyy-mm-dd HH:MM:SS zzzzz")
    local expected = [
        Dates.Year(2017),
        Dates.Month(11),
        Dates.Day(14),
        Dates.Hour(11),
        Dates.Minute(3),
        Dates.Second(53),
        tz"UTC+01",
    ]
    @test parse_components(test...) == expected
end

@testset "default format" begin
    @test default_format(ZonedDateTime) === TimeZones.ISOZonedDateTimeFormat
end

@testset "parse constructor" begin
    @test isequal(
        ZonedDateTime("2000-01-02T03:04:05.006+0700"),
        ZonedDateTime(2000, 1, 2, 3, 4, 5, 6, tz"UTC+07")
    )
    @test isequal(
        ZonedDateTime("2018-11-01-0600", dateformat"yyyy-mm-ddzzzz"),
        ZonedDateTime(2018, 11, 1, tz"UTC-06"),
    )

    # Validate that error message contains the original string and the format used
    str = "2018-11-01"
    try
        ZonedDateTime(str)
    catch e
        @test e isa ArgumentError
        @test occursin(str, e.msg)
        @test occursin(string(TimeZones.ISOZonedDateTimeFormat), e.msg)
    end
end


text:
using Documenter
using RandomTree

makedocs(
    sitename="RandomTree.jl",
    modules = [RandomTree],
    pages = [
        "index.md",
        "Installation" => "install.md",
        "Usage" => "usage.md",
        "API reference" => "api.md",
        "Future Plan" => "plan.md",
    ]
)

deploydocs(
    repo = "github.com/newptcai/RandomTree.jl.git",
)


text:
"""
    TypeL(var_1::TypeM)

Lorem ipsum dolor sit amet.
"""
struct TypeL

  var_1::TypeM

end