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

using Base.Pkg.Types
using Base.Pkg.Query
using Base.Pkg.Resolve
using Base.Pkg.Resolve.VersionWeights
import Base.Pkg.PkgError

# Check that VersionWeight keeps the same ordering as VersionNumber

vlst = [
    v"0.0.0",
    v"0.0.1",
    v"0.1.0",
    v"0.1.1",
    v"1.0.0",
    v"1.0.1",
    v"1.1.0",
    v"1.1.1",
    v"1.0.0-pre",
    v"1.0.0-pre1",
    v"1.0.1-pre",
    v"1.0.0-0.pre.2",
    v"1.0.0-0.pre.3",
    v"1.0.0-0.pre1.tst",
    v"1.0.0-pre.1+0.1",
    v"1.0.0-pre.1+0.1plus",
    v"1.0.0-pre.1-+0.1plus",
    v"1.0.0-pre.1-+0.1Plus",
    v"1.0.0-pre.1-+0.1pLUs",
    v"1.0.0-pre.1-+0.1pluS",
    v"1.0.0+0.1plus",
    v"1.0.0+0.1plus-",
    v"1.0.0+-",
    v"1.0.0-",
    v"1.0.0+",
    v"1.0.0--",
    v"1.0.0---",
    v"1.0.0--+-",
    v"1.0.0+--",
    v"1.0.0+-.-",
    v"1.0.0+0.-",
    v"1.0.0+-.0",
    v"1.0.0-a+--",
    v"1.0.0-a+-.-",
    v"1.0.0-a+0.-",
    v"1.0.0-a+-.0"
    ]

for v1 in vlst, v2 in vlst
    vw1 = VersionWeight(v1)
    vw2 = VersionWeight(v2)
    clt = v1 < v2
    @test clt == (vw1 < vw2)
    ceq = v1 == v2
    @test ceq == (vw1 == vw2)
end

# auxiliary functions

function deps_from_data(deps_data)
    deps = Dict{String,Dict{VersionNumber,Available}}()
    for d in deps_data
        p = d[1]; vn = d[2]; r = d[3:end]
        if !haskey(deps, p)
            deps[p] = Dict{VersionNumber,Available}()
        end
        if !haskey(deps[p], vn)
            deps[p][vn] = Available("$(p)_$(vn)_sha1", Dict{String,VersionSet}())
        end
        isempty(r) && continue
        rp = r[1]
        if length(r) > 1
            rvs = VersionSet(VersionNumber[r[2:end]...])
        else
            rvs = VersionSet()
        end
        deps[p][vn].requires[rp] = rvs
    end
    deps
end
function reqs_from_data(reqs_data)
    reqs = Dict{String,VersionSet}()
    for r in reqs_data
        p = r[1]
        reqs[p] = VersionSet(VersionNumber[r[2:end]...])
    end
    reqs
end
function sanity_tst(deps_data, expected_result; pkgs=[])
    deps = deps_from_data(deps_data)
    #println("deps=$deps")
    #println()
    result = sanity_check(deps, Set(String[pkgs...]))
    length(result) == length(expected_result) || return false
    for (p, vn, pp) in result
        in((p, vn), expected_result) || return  false
    end
    return true
end
sanity_tst(deps_data; kw...) = sanity_tst(deps_data, []; kw...)

function resolve_tst(deps_data, reqs_data)
    deps = deps_from_data(deps_data)
    reqs = reqs_from_data(reqs_data)

    #println()
    #println("deps=$deps")
    #println("reqs=$reqs")
    deps = Query.prune_dependencies(reqs, deps)
    want = resolve(reqs, deps)
end

## DEPENDENCY SCHEME 1: TWO PACKAGES, DAG
deps_data = Any[
    ["A", v"1", "B", v"1"],
    ["A", v"2", "B", v"2"],
    ["B", v"1"],
    ["B", v"2"]
]

@test sanity_tst(deps_data)
@test sanity_tst(deps_data, pkgs=["A", "B"])
@test sanity_tst(deps_data, pkgs=["B"])
@test sanity_tst(deps_data, pkgs=["A"])

# require just B
reqs_data = Any[
    ["B"]
]

want = resolve_tst(deps_data, reqs_data)
@test want == Dict("B"=>v"2")

# require just A: must bring in B
reqs_data = Any[
    ["A"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"2", "B"=>v"2")

## DEPENDENCY SCHEME 2: TWO PACKAGES, CYCLIC
deps_data = Any[
    ["A", v"1", "B", v"2"],
    ["A", v"2", "B", v"1"],
    ["B", v"1", "A", v"2"],
    ["B", v"2", "A", v"1"]
]

@test sanity_tst(deps_data)

# require just A
reqs_data = Any[
    ["A"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"2", "B"=>v"2")

# require just B, force lower version
reqs_data = Any[
    ["B", v"1", v"2"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"2", "B"=>v"1")

# require just A, force lower version
reqs_data = Any[
    ["A", v"1", v"2"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"1", "B"=>v"2")

## DEPENDENCY SCHEME 3: THREE PACKAGES, CYCLIC, TWO MUTUALLY EXCLUSIVE SOLUTIONS
deps_data = Any[
    ["A", v"1", "B", v"2"],
    ["A", v"2", "B", v"1", v"2"],
    ["B", v"1", "C", v"2"],
    ["B", v"2", "C", v"1", v"2"],
    ["C", v"1", "A", v"1", v"2"],
    ["C", v"2", "A", v"2"]
]

@test sanity_tst(deps_data)

# require just A (must choose solution which has the highest version for A)
reqs_data = Any[
    ["A"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"2", "B"=>v"1", "C"=>v"2")

# require just B (must choose solution which has the highest version for B)
reqs_data = Any[
    ["B"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"1", "B"=>v"2", "C"=>v"1")

# require just A, force lower version
reqs_data = Any[
    ["A", v"1", v"2"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"1", "B"=>v"2", "C"=>v"1")

# require A and C, incompatible versions
reqs_data = Any[
    ["A", v"1", v"2"],
    ["C", v"2"]
]
@test_throws PkgError resolve_tst(deps_data, reqs_data)

## DEPENDENCY SCHEME 4: TWO PACKAGES, DAG, WITH TRIVIAL INCONSISTENCY
deps_data = Any[
    ["A", v"1", "B", v"2"],
    ["B", v"1"]
]

@test sanity_tst(deps_data, [("A", v"1")])
@test sanity_tst(deps_data, [("A", v"1")], pkgs=["B"])

# require B (must not give errors)
reqs_data = Any[
    ["B"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("B"=>v"1")

## DEPENDENCY SCHEME 5: THREE PACKAGES, DAG, WITH IMPLICIT INCONSISTENCY
deps_data = Any[
    ["A", v"1", "B", v"2"],
    ["A", v"1", "C", v"2"],
    ["A", v"2", "B", v"1", v"2"],
    ["A", v"2", "C", v"1", v"2"],
    ["B", v"1", "C", v"2"],
    ["B", v"2", "C", v"2"],
    ["C", v"1"],
    ["C", v"2"]
]

@test sanity_tst(deps_data, [("A", v"2")])
@test sanity_tst(deps_data, [("A", v"2")], pkgs=["B"])
@test sanity_tst(deps_data, [("A", v"2")], pkgs=["C"])

# require A, any version (must use the highest non-inconsistent)
reqs_data = Any[
    ["A"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"1", "B"=>v"2", "C"=>v"2")

# require A, force highest version (impossible)
reqs_data = Any[
    ["A", v"2"]
]
@test_throws PkgError resolve_tst(deps_data, reqs_data)

## DEPENDENCY SCHEME 6: TWO PACKAGES, CYCLIC, TOTALLY INCONSISTENT
deps_data = Any[
    ["A", v"1", "B", v"2"],
    ["A", v"2", "B", v"1", v"2"],
    ["B", v"1", "A", v"1", v"2"],
    ["B", v"2", "A", v"2"]
]

@test sanity_tst(deps_data, [("A", v"1"), ("A", v"2"),
                             ("B", v"1"), ("B", v"2")])

# require A (impossible)
reqs_data = Any[
    ["A"]
]
@test_throws PkgError resolve_tst(deps_data, reqs_data)

# require B (impossible)
reqs_data = Any[
    ["B"]
]
@test_throws PkgError resolve_tst(deps_data, reqs_data)

## DEPENDENCY SCHEME 7: THREE PACKAGES, CYCLIC, WITH INCONSISTENCY
deps_data = Any[
    ["A", v"1", "B", v"1", v"2"],
    ["A", v"2", "B", v"2"],
    ["B", v"1", "C", v"1", v"2"],
    ["B", v"2", "C", v"2"],
    ["C", v"1", "A", v"2"],
    ["C", v"2", "A", v"2"],
]

@test sanity_tst(deps_data, [("A", v"1"), ("B", v"1"),
                             ("C", v"1")])

# require A
reqs_data = Any[
    ["A"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"2", "B"=>v"2", "C"=>v"2")

# require C
reqs_data = Any[
    ["C"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"2", "B"=>v"2", "C"=>v"2")

# require C, lowest version (impossible)
reqs_data = Any[
    ["C", v"1", v"2"]
]
@test_throws PkgError resolve_tst(deps_data, reqs_data)

## DEPENDENCY SCHEME 8: THREE PACKAGES, CYCLIC, TOTALLY INCONSISTENT
deps_data = Any[
    ["A", v"1", "B", v"1", v"2"],
    ["A", v"2", "B", v"2"],
    ["B", v"1", "C", v"1", v"2"],
    ["B", v"2", "C", v"2"],
    ["C", v"1", "A", v"2"],
    ["C", v"2", "A", v"1", v"2"],
]

@test sanity_tst(deps_data, [("A", v"1"), ("A", v"2"),
                             ("B", v"1"), ("B", v"2"),
                             ("C", v"1"), ("C", v"2")])

# require A (impossible)
reqs_data = Any[
    ["A"]
]
@test_throws PkgError resolve_tst(deps_data, reqs_data)

# require B (impossible)
reqs_data = Any[
    ["B"]
]
@test_throws PkgError resolve_tst(deps_data, reqs_data)

# require C (impossible)
reqs_data = Any[
    ["C"]
]
@test_throws PkgError resolve_tst(deps_data, reqs_data)

## DEPENDENCY SCHEME 9: SIX PACKAGES, DAG
deps_data = Any[
    ["A", v"1"],
    ["A", v"2"],
    ["A", v"3"],
    ["B", v"1", "A", v"1", v"2"],
    ["B", v"2", "A"],
    ["C", v"1", "A", v"2", v"3"],
    ["C", v"2", "A", v"2"],
    ["D", v"1", "B", v"1"],
    ["D", v"2", "B", v"2"],
    ["E", v"1", "D"],
    ["F", v"1", "A", v"1", v"3"],
    ["F", v"1", "E"],
    ["F", v"2", "C", v"2"],
    ["F", v"2", "E"],
]

@test sanity_tst(deps_data)

# require just F
reqs_data = Any[
    ["F"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"3", "B"=>v"2", "C"=>v"2",
                   "D"=>v"2", "E"=>v"1", "F"=>v"2")

# require just F, lower version
reqs_data = Any[
    ["F", v"1", v"2"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"2", "B"=>v"2", "D"=>v"2",
                   "E"=>v"1", "F"=>v"1")

# require F and B; force lower B version -> must bring down F, A, and D versions too
reqs_data = Any[
    ["F"],
    ["B", v"1", v"2"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"1", "B"=>v"1", "D"=>v"1",
                   "E"=>v"1", "F"=>v"1")

# require F and D; force lower D version -> must not bring down F version
reqs_data = Any[
    ["F"],
    ["D", v"1", v"2"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"3", "B"=>v"2", "C"=>v"2",
                   "D"=>v"1", "E"=>v"1", "F"=>v"2")

# require F and C; force lower C version -> must bring down F and A versions
reqs_data = Any[
    ["F"],
    ["C", v"1", v"2"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"2", "B"=>v"2", "C"=>v"1",
                   "D"=>v"2", "E"=>v"1", "F"=>v"1")

## DEPENDENCY SCHEME 10: SIX PACKAGES, DAG, WITH PRERELEASE/BUILD VERSIONS
deps_data = Any[
    ["A", v"1"],
    ["A", v"2-rc.1"],
    ["A", v"2-rc.1+bld"],
    ["A", v"2"],
    ["A", v"2.1.0"],
    ["B", v"1", "A", v"1", v"2-"],
    ["B", v"1.0.1-beta", "A", v"2-rc"],
    ["B", v"1.0.1", "A"],
    ["C", v"1", "A", v"2-", v"2.1"],
    ["C", v"1+BLD", "A", v"2-rc.1", v"2.1"],
    ["C", v"2", "A", v"2-rc.1"],
    ["D", v"1", "B", v"1"],
    ["D", v"2", "B", v"1.0.1-"],
    ["E", v"1-plztst", "D"],
    ["E", v"1", "D"],
    ["F", v"1.1", "A", v"1", v"2.1"],
    ["F", v"1.1", "E", v"1"],
    ["F", v"2-rc.1", "A", v"2-", v"2.1"],
    ["F", v"2-rc.1", "C", v"1"],
    ["F", v"2-rc.1", "E"],
    ["F", v"2-rc.2", "A", v"2-", v"2.1"],
    ["F", v"2-rc.2", "C", v"2"],
    ["F", v"2-rc.2", "E"],
    ["F", v"2", "C", v"2"],
    ["F", v"2", "E"],
]

@test sanity_tst(deps_data)

# require just F
reqs_data = Any[
    ["F"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"2.1", "B"=>v"1.0.1", "C"=>v"2",
                   "D"=>v"2", "E"=>v"1", "F"=>v"2")

# require just F, lower version
reqs_data = Any[
    ["F", v"1", v"2-"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"2", "B"=>v"1.0.1", "D"=>v"2",
                   "E"=>v"1", "F"=>v"1.1")

# require F and B; force lower B version -> must bring down F, A, and D versions too
reqs_data = Any[
    ["F"],
    ["B", v"1", v"1.0.1-"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"1", "B"=>v"1", "D"=>v"1",
                   "E"=>v"1", "F"=>v"1.1")

# require F and D; force lower D version -> must not bring down F version
reqs_data = Any[
    ["F"],
    ["D", v"1", v"2"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"2.1", "B"=>v"1.0.1", "C"=>v"2",
                   "D"=>v"1", "E"=>v"1", "F"=>v"2")

# require F and C; force lower C version -> must bring down F and A versions
reqs_data = Any[
    ["F"],
    ["C", v"1", v"2"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"2", "B"=>v"1.0.1", "C"=>v"1+BLD",
                   "D"=>v"2", "E"=>v"1", "F"=>v"2-rc.1")

## DEPENDENCY SCHEME 11: FIVE PACKAGES, SAME AS SCHEMES 5 + 1, UNCONNECTED
deps_data = Any[
    ["A", v"1", "B", v"2"],
    ["A", v"1", "C", v"2"],
    ["A", v"2", "B", v"1", v"2"],
    ["A", v"2", "C", v"1", v"2"],
    ["B", v"1", "C", v"2"],
    ["B", v"2", "C", v"2"],
    ["C", v"1"],
    ["C", v"2"],
    ["D", v"1", "E", v"1"],
    ["D", v"2", "E", v"2"],
    ["E", v"1"],
    ["E", v"2"]
]

@test sanity_tst(deps_data, [("A", v"2")])
@test sanity_tst(deps_data, [("A", v"2")], pkgs=["B"])
@test sanity_tst(deps_data, pkgs=["D"])
@test sanity_tst(deps_data, pkgs=["E"])
@test sanity_tst(deps_data, [("A", v"2")], pkgs=["B", "D"])

# require A, any version (must use the highest non-inconsistent)
reqs_data = Any[
    ["A"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"1", "B"=>v"2", "C"=>v"2")

# require just D: must bring in E
reqs_data = Any[
    ["D"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("D"=>v"2", "E"=>v"2")

# require A and D, must be the merge of the previous two cases
reqs_data = Any[
    ["A"],
    ["D"]
]
want = resolve_tst(deps_data, reqs_data)
@test want == Dict("A"=>v"1", "B"=>v"2", "C"=>v"2", "D"=>v"2", "E"=>v"2")


text:
Base.convert{T}(::Type{Matrix{Float64}}, img::Image{Gray{T}}) =
    convert(Array{Float64, 2}, convert(Array, img))

# Base.convert(::Type{Array{Float64, 3}}, img::Image{RGB}) =
#     convert(Array{Float64}, separate(img).data)

rgb2arr(img) = convert(Array{Float64}, separate(img).data)

function walkdir(dir; pred=(_ -> true))
    paths = [joinpath(dir, filename) for filename in readdir(dir)]
    result = Array(AbstractString, 0)
    for path in paths
        if isdir(path)
            append!(result, walkdir(path, pred=pred))
        elseif pred(path)
            push!(result, path)
        end
    end
    return result
end

ext: jl


text:
using RecipesBase

#=
The plotting functions in this file are adapted from those in LazySets.jl.
See the relevant docstrings in LazySets for the description of the available options.
=#

using LazySets: plot_recipe,
                PLOT_PRECISION,
                PLOT_POLAR_DIRECTIONS,
                DEFAULT_LABEL,
                DEFAULT_GRID,
                DEFAULT_ASPECT_RATIO,
                DEFAULT_ALPHA,
                DEFAULT_COLOR,
                _extract_limits,
                _extract_extrema,
                _set_auto_limits_to_extrema!,
                _bounding_hyperrectangle,
                _update_plot_limits!,
                isapproxzero,
                _plot_singleton_list

import LazySets: _plot_singleton_list_1D,
                 _plot_singleton_list_2D

# heuristics for projecting a reach-set either concretely or lazily to the space
# spanned by vars; the returned set is concrete when the exact projection can be done
# efficiently; in all other cases a lazy set is returned
function _project_reachset(R::AbstractLazyReachSet, vars)
    ST = setrep(R)
    return _project_reachset(ST, R, vars)
end

# zonotopes are projected concretely unless they have too many generators
function _project_reachset(::Type{<:Zonotope}, R, vars)
    # TODO review/lift the condition on max number of generators after LazySets#2288
    if (ngens(set(R)) <= 15) || (0 ∈ vars)
        # concrete projection is efficient
        πR = project(R, vars) |> set

        # zonotopes usually contain lots of redundant generators
        πR = remove_redundant_generators(πR)
    else
        # avoid expensive vertex enumeration
        πR = Projection(R, vars) |> set
    end
    return πR
end

function _project_reachset(::Union{Type{<:AbstractZonotope}, VPOLY}, R, vars)
    # concrete projection is efficient
    πR = project(R, vars)
    return set(πR)
end

function _project_reachset(::Type{<:AbstractPolyhedron}, R, vars)
    # if the set is polyhedral and two-dimensional
    #   - if time is not required -> don't project
    #   - otherwise, make the projection (lazily)
    if (dim(R) == 2) &&  (0 ∉ vars)
        πR = set(R) # no-op
    else
        πR = Projection(R, vars) |> set # lazy projection
    end
    return πR
end

function _project_reachset(::Type{<:LazySet}, R, vars)
    πR = Projection(R, vars) # lazy projection
    X = set(πR)
end

# ---------------------------------------
# Specialization for template reach-sets
# ---------------------------------------
using LazySets: _check_constrained_dimensions

function _is_compatible_projection(R::TemplateReachSet, vars)
    X = set(R)
    for ci in constraints_list(X)
        status = _check_constrained_dimensions(ci, vars)
        status == 0 && return false
    end
    return true
end

# concrete projection for template reach-sets
# TODO specialize for vars = (0, v) and when the direction v is known from a SEV template
function _project_reachset(R::TemplateReachSet, vars)
    # if the exact projection is cheap, we compute it;
    # otherwise we take the lazy projection
    # in both cases 2D bounded sets are eps-close approximated
    if _is_compatible_projection(R, vars)
        πR = project(R, vars)
    else
       πR = Projection(R, vars)
    end
    X = set(πR)
    if dim(X) == 2 && isbounded(X)
        X = overapproximate(X, HPolygon, 1e-3)
    end
    return X
end

# --------------------------------------------
# Specialization for Taylor model reach-sets
# --------------------------------------------

function _project_reachset(T::TaylorModelReachSet, vars)
    R = overapproximate(T, Zonotope)
    _project_reachset(R, vars)
end

function _check_vars(vars)
    if vars == nothing
        throw(ArgumentError("default plotting variables not implemented yet; you need " *
              "to pass the `vars=(...)` option, e.g. `vars=(0, 1)` to plot variable with " *
              "index 1 vs. time, or `vars=(1, 2)` to plot variable with index 2 vs. variable with index 1`"))
    end
    D = length(vars)
    @assert (D == 1) || (D == 2) "can only plot in one or two dimensions, " *
                                 "but received $D variable indices where `vars = ` $vars"
end

# ========================
# Reach-set plot recipes
# ========================

@recipe function plot_reachset(R::AbstractReachSet{N};
                               vars=nothing,
                               ε=N(PLOT_PRECISION)) where {N}
    _check_vars(vars)
    X = _project_reachset(R, vars)

    dim(X) == 1 && return plot_recipe(X, ε)

    label --> DEFAULT_LABEL
    grid --> DEFAULT_GRID
    if DEFAULT_ASPECT_RATIO != :none
        aspect_ratio --> DEFAULT_ASPECT_RATIO
    end
    seriesalpha --> DEFAULT_ALPHA
    seriescolor --> DEFAULT_COLOR

    # extract limits and extrema of already plotted sets
    p = plotattributes[:plot_object]
    lims = _extract_limits(p, plotattributes)
    extr = _extract_extrema(p)

    if !isbounded(X)
        _set_auto_limits_to_extrema!(lims, extr)
        X = intersection(X, _bounding_hyperrectangle(lims, eltype(X)))

    elseif length(p) > 0
        # if there is already a plotted set and the limits are fixed,
        # automatically adjust the axis limits (e.g. after plotting a unbounded set)
        _update_plot_limits!(lims, X)
    end

    xlims --> lims[:x]
    ylims --> lims[:y]

    x, y = plot_recipe(X, ε)
    if !isempty(x)
        m = length(x)
        if m == 1
            seriestype := :scatter
        elseif m == 2 && isapproxzero((x[1] - x[2])^2 + (y[1] - y[2])^2)
            seriestype := :scatter
        else
            seriestype := :shape
        end
    end
    return (x, y)
end

# TODO use LazySets._plot_singleton_list_1D if the reach-sets are stored as a struct
# array. Otherwise, we could use pass the list through x -> set(x)
function _plot_singleton_list_1D(list::Union{Flowpipe{N}, AbstractVector{RN}}) where {N, RN<:AbstractReachSet{N}}
    m = length(list)

    x = Vector{N}(undef, m)
    y = zeros(N, m)

    @inbounds for (i, Xi) in enumerate(list)
        p = element(set(Xi))
        x[i] = p[1]
    end
    x, y
end

function _plot_singleton_list_2D(list::Union{Flowpipe{N}, AbstractVector{RN}}) where {N, RN<:AbstractReachSet{N}}
    m = length(list)
    x = Vector{N}(undef, m)
    y = Vector{N}(undef, m)

    @inbounds for (i, Xi) in enumerate(list)
        p = element(set(Xi))
        x[i] = p[1]
        y[i] = p[2]
    end
    x, y
end

function _plot_reachset_list(list, N, vars, ε, Nφ)

    first = true
    x = Vector{N}()
    y = Vector{N}()
    for Ri in list
        Xi = _project_reachset(Ri, vars)

        if Xi isa Intersection
            xcoords, ycoords = plot_recipe(Xi, ε, Nφ)

        else
            Pi = isoperation(Xi) ? overapproximate(Xi, ε) : Xi
            vlist = convex_hull(vertices_list(Pi))
            m = length(vlist)
            if m == 0
                @warn "overapproximation during plotting was empty"
                continue
            end
            xcoords = Vector{N}(undef, m)
            ycoords = Vector{N}(undef, m)
            @inbounds for (i, v) in enumerate(vlist)
                xcoords[i] = v[1]
                ycoords[i] = v[2]
            end
            if m > 1
                # add first vertex to "close" the polygon
                push!(xcoords, xcoords[1])
                push!(ycoords, ycoords[1])
            end
        end
        isempty(xcoords) && continue

        x_new = xcoords
        y_new = ycoords

        if first
            first = false
        else
            push!(x, N(NaN))
            push!(y, N(NaN))
        end
        append!(x, x_new)
        append!(y, y_new)
    end
    x, y
end

# ========================
# Flowpipe plot recipes
# ========================

@recipe function plot_list(list::Union{Flowpipe{N}, AbstractVector{RN}};
                           vars=nothing,
                           ε=N(PLOT_PRECISION),
                           Nφ=PLOT_POLAR_DIRECTIONS) where {N, RN<:AbstractReachSet{N}}
    _check_vars(vars)

    label --> DEFAULT_LABEL
    grid --> DEFAULT_GRID
    if DEFAULT_ASPECT_RATIO != :none
        aspect_ratio --> DEFAULT_ASPECT_RATIO
    end
    seriesalpha --> DEFAULT_ALPHA
    seriescolor --> DEFAULT_COLOR

    if !(0 ∈ vars) && (setrep(list) <: AbstractSingleton)
        seriestype --> :scatter
        _plot_singleton_list(list)
    else
        seriestype --> :shape
        _plot_reachset_list(list, N, vars, ε, Nφ)
    end
end

# composite flowpipes
@recipe function plot_list(fp::Union{HF, MF};
                           vars=nothing,
                           ε=Float64(PLOT_PRECISION),
                           Nφ=PLOT_POLAR_DIRECTIONS
                          ) where {N, HF<:HybridFlowpipe{N}, MF<:MixedFlowpipe{N}}

    _check_vars(vars)

    label --> DEFAULT_LABEL
    grid --> DEFAULT_GRID
    if DEFAULT_ASPECT_RATIO != :none
        aspect_ratio --> DEFAULT_ASPECT_RATIO
    end
    seriesalpha --> DEFAULT_ALPHA
    seriescolor --> DEFAULT_COLOR
    seriestype --> :shape

    x = Vector{N}()
    y = Vector{N}()

    for F in fp
        x_new, y_new = _plot_reachset_list(F, N, vars, ε, Nφ)
        append!(x, x_new)
        append!(y, y_new)
        push!(x, N(NaN))
        push!(y, N(NaN))
    end
    x, y
end

# ========================
# Solution plot recipes
# ========================

@recipe function plot_list(sol::ReachSolution{FT};
                           vars=nothing,
                           ε=Float64(PLOT_PRECISION),
                           Nφ=PLOT_POLAR_DIRECTIONS
                          ) where {N, FT<:Flowpipe{N}}

    _check_vars(vars)

    label --> DEFAULT_LABEL
    grid --> DEFAULT_GRID
    if DEFAULT_ASPECT_RATIO != :none
        aspect_ratio --> DEFAULT_ASPECT_RATIO
    end
    seriesalpha --> DEFAULT_ALPHA
    seriescolor --> DEFAULT_COLOR

    fp = flowpipe(sol)
    if !(0 ∈ vars) && (setrep(fp) <: AbstractSingleton)
        seriestype --> :scatter
        _plot_singleton_list(fp)
    else
        seriestype --> :shape
        _plot_reachset_list(fp, N, vars, ε, Nφ)
    end
end

# compound solution flowpipes
@recipe function plot_list(sol::Union{SMF, SHF};
                           vars=nothing,
                           ε=Float64(PLOT_PRECISION),
                           Nφ=PLOT_POLAR_DIRECTIONS
                          ) where {N, MF<:MixedFlowpipe{N}, HF<:HybridFlowpipe{N},
                                   SMF<:ReachSolution{MF}, SHF<:ReachSolution{HF}}
    _check_vars(vars)

    label --> DEFAULT_LABEL
    grid --> DEFAULT_GRID
    if DEFAULT_ASPECT_RATIO != :none
        aspect_ratio --> DEFAULT_ASPECT_RATIO
    end
    seriesalpha --> DEFAULT_ALPHA
    seriescolor --> DEFAULT_COLOR
    seriestype --> :shape

    x = Vector{N}()
    y = Vector{N}()

    fp = flowpipe(sol)
    for F in fp
        x_new, y_new = _plot_reachset_list(F, N, vars, ε, Nφ)
        append!(x, x_new)
        append!(y, y_new)
        push!(x, N(NaN))
        push!(y, N(NaN))
    end
    x, y
end

# TODO new plot recipe to dispatch on ShiftedFlowpipe
# πRi = project(F, i, vars) # project the reach-set

# TODO: refactor with Flowpipe
# TODO extend projection of shifted flowpipes to use lazy projection if needed
@recipe function plot_list(fp::ShiftedFlowpipe{N};
                           vars=nothing,
                           ε=Float64(PLOT_PRECISION),
                           Nφ=PLOT_POLAR_DIRECTIONS
                          ) where {N}
    _check_vars(vars)

    label --> DEFAULT_LABEL
    grid --> DEFAULT_GRID
    if DEFAULT_ASPECT_RATIO != :none
        aspect_ratio --> DEFAULT_ASPECT_RATIO
    end
    seriesalpha --> DEFAULT_ALPHA
    seriescolor --> DEFAULT_COLOR
    seriestype --> :shape

    first = true
    x = Vector{N}()
    y = Vector{N}()
    for (i, Ri) in enumerate(fp)
        πRi = project(fp, i, vars) # project the reach-set
        Xi = set(πRi) # extract the set representation

        if Xi isa Intersection
            res = plot_recipe(Xi, ε, Nφ)
        else
            # hard-code overapproximation here to avoid individual
            # compilations for mixed sets
            Pi = isa(Xi, AbstractPolygon) ? Xi : overapproximate(Xi, ε)
            vlist = transpose(hcat(convex_hull(vertices_list(Pi))...))
            if isempty(vlist)
                @warn "overapproximation during plotting was empty"
                continue
            end
            res = vlist[:, 1], vlist[:, 2]
            # add first vertex to "close" the polygon
            push!(res[1], vlist[1, 1])
            push!(res[2], vlist[1, 2])
        end
        if isempty(res)
            continue
        else
            x_new, y_new = res
        end
        if first
            first = false
        else
            push!(x, N(NaN))
            push!(y, N(NaN))
        end
        append!(x, x_new)
        append!(y, y_new)
    end
    x, y
end

#=
# TODO: default plotting without vars definition
# plot each variable vs time
@recipe function plot_reachset(R::AbstractLazyReachSet{N};
                               ε::N=N(PLOT_PRECISION)
                               ) where {N<:Real, D, M<:Integer}
    error("not implemented yet")
end

@recipe function plot_list(list::AbstractVector{RN};
                           ε::N=N(PLOT_PRECISION),
                           Nφ::Int=PLOT_POLAR_DIRECTIONS,
                           fast::Bool=true
                          ) where {N<:Real, RN<:AbstractReachSet{N}}
    error("not implemented yet")
end
=#

# solution from computation without bloating and singleton initial condition
# (eg. with ORBIT) is presened as a scatter plot
@recipe function plot_list(sol::ReachSolution{FT, <:ORBIT};
                           vars=nothing) where {N, RT<:ReachSet{N, <:Singleton{N}}, FT<:Flowpipe{N, RT}}
   label --> DEFAULT_LABEL
   grid --> DEFAULT_GRID
   if DEFAULT_ASPECT_RATIO != :none
       aspect_ratio --> DEFAULT_ASPECT_RATIO
   end
   seriesalpha --> DEFAULT_ALPHA
   seriescolor --> DEFAULT_COLOR
   seriestype --> :scatter
   markershape --> :circle

   _check_vars(vars)
   vx, vy = vars[1], vars[2]
   X(k) = (vx == 0) ? tstart(sol[k]) : element(set(sol[k]))[vx]
   Y(k) = (vy == 0) ? tstart(sol[k]) : element(set(sol[k]))[vy]

   x = Vector{N}()
   y = Vector{N}()
   for k in 1:length(sol)
       x_new = X(k)
       y_new = Y(k)
       append!(x, x_new)
       append!(y, y_new)
   end
   x, y
end


text:
module AbstractDataStructure

"""
A mutable data structure addressed by type Ref, holding values of type Value.
"""
abstract type DataStructure{Ref, Value} end

"""
Predicate used to indicate that a data structure is inherently ordered.

ordered_data_structure_p(ds::DataStructure) -> Bool
"""
function ordered_data_structure_p end

"""
Insert a value in the data structure and return a reference to its location.
The reference is only valid in the current state of the data structure.

ins!(ds::DataStructure, val::Value) -> Ref
"""
function ins! end

"""
Delete an element from the data strucutre at lotation pos.
Return the value that has been deleted

del!(ds::DataStructure, pos::Ref) -> Value
raises KeyError
"""
function del! end

"""
Return the value of the element at lotation pos.

at(ds::DataStructure, pos::Ref) -> Value
raises KeyError
"""
function at end

"""
Find the locations of all elements with value val.
Return a vector of all found references.

search(ds::DataStructure, val::Value) -> Vector{Ref}
"""
function search end

"""
Return the number of elements in the data structure

len(ds::DataStructure) -> Integer
"""
function len end

Base.length(ds::DataStructure) = len(ds)

end  # module AbstractDataStructure


text:
## --- Query TC1 properties

    tc1_550 = readdlm(joinpath(moduleresourcepath,"tc1","tc1_550.csv"), ',')
    tc1_1300 = readdlm(joinpath(moduleresourcepath,"tc1","tc1_1300.csv"), ',')
    tc1_age = readdlm(joinpath(moduleresourcepath,"tc1","tc1_age.csv"), ',', Int)

    """
    `find_tc1_crust(lat::Number,lon::Number)`

    Find the depth to the 550C isotherm for the 1x1 arc degree grid cell
    containing `lat` and `lon`

    `find_tc1_crust(lat::AbstractArray,lon::AbstractArray)`

    For each pair of latitudes and longitudes given by lat and lon, find the
    depth to the 550C isotherm for the 1x1 arc degree grid cell containing
    `lat[i]` and `lon[i]`
    """
    function find_tc1_crust(lat::Number,lon::Number)
        if !isnan(lat) && !isnan(lon)
            i = round(Int, 91-lat)
            j = round(Int, lon+181)
            crust=tc1_550[i,j]
        else
            crust = NaN
        end
        return crust
    end
    function find_tc1_crust(lat::AbstractArray,lon::AbstractArray)
        # Check input dimensions match
        if length(lat) != length(lon)
            error("lat and lon must be equal length\n")
        end

        # Query the tc1_550 array for our lat and lon
        crust = fill(NaN, size(lat))
        for n=1:length(lat)
            if !isnan(lat[n]) && !isnan(lon[n])
                i = round(Int, 91-lat[n])
                j = round(Int, lon[n]+181)
                crust[n]=tc1_550[i,j]
            end
        end
        return crust
    end
    export find_tc1_crust

    """
    `find_tc1_lith(lat::Number,lon::Number)`

    Find the depth to the 1300C isotherm for the 1x1 arc degree grid cell
    containing `lat` and `lon`

    `find_tc1_lith(lat::AbstractArray,lon::AbstractArray)`

    For each pair of latitudes and longitudes given by lat and lon, find the
    depth to the 1300C isotherm for the 1x1 arc degree grid cell containing
    `lat[i]` and `lon[i]`
    """
    function find_tc1_lith(lat::Number,lon::Number)
        if !isnan(lat) && !isnan(lon)
            i = round(Int, 91-lat)
            j = round(Int, lon+181)
            lith=tc1_1300[i,j]
        else
            lith = NaN
        end
        return lith
    end
    function find_tc1_lith(lat::AbstractArray,lon::AbstractArray)
        # Check input dimensions match
        if length(lat) != length(lon)
            error("lat and lon must be equal length\n")
        end

        # Query the tc1_1300 array for our lat and lon
        lith = fill(NaN, size(lat))
        for n=1:length(lat)
            if !isnan(lat[n]) && !isnan(lon[n])
                i = round(Int, 91-lat[n])
                j = round(Int, lon[n]+181)
                lith[n]=tc1_1300[i,j]
            end
        end
        return lith
    end
    export find_tc1_lith

    """
    `find_tc1_age(lat::Number,lon::Number)`

    Return a tuple `(age, age_min, age_max)` containing the nominal, upper, and
    lower tc1 age bounds for the 1x1 arc degree grid cell containing `lat`
    and `lon`

    `find_tc1_age(lat::AbstractArray,lon::AbstractArray)`

    Return a tuple `(age, age_min, age_max)` where `age`, `age_min`, and
    `age_max` are arrays containing the nominal, upper and lower tc1 age bounds
    for each location pair `lat[i]`, `lon[i]`
    """
    function find_tc1_age(lat::Number,lon::Number)
        ages=[ NaN  NaN  NaN
                25    0   50
               150   50  250
               395  250  540
               695  540  850
               975  850 1100
              1400 1100 1700
              2100 1700 2500
              2750 2500 3000
              3250 3000 3500]

        if !isnan(lat) && !isnan(lon)
            i = round(Int, 91-lat)
            j = round(Int, lon+181)
            result = (ages[tc1_age[i,j],:]...,)
        else
            result = (NaN, NaN, NaN,)
        end
        return result
    end
    function find_tc1_age(lat::AbstractArray,lon::AbstractArray)
        # Check input dimensions match
        if length(lat) != length(lon)
            error("lat and lon must be equal length\n")
        end

        ages=[ NaN  NaN  NaN
                25    0   50
               150   50  250
               395  250  540
               695  540  850
               975  850 1100
              1400 1100 1700
              2100 1700 2500
              2750 2500 3000
              3250 3000 3500]

        # Query the tc1_1300 array for our lat and lon
        age = fill(NaN, size(lat))
        minage = fill(NaN, size(lat))
        maxage = fill(NaN, size(lat))
        for n=1:length(lat)
            if !isnan(lat[n]) && !isnan(lon[n])
                i = round(Int, 91-lat[n])
                j = round(Int, lon[n]+181)
                ageindex = tc1_age[i,j]
                age[n] = ages[ageindex,1]
                minage[n] = ages[ageindex,2]
                maxage[n] = ages[ageindex,3]
            end
        end
        return (age, minage, maxage)
    end
    export find_tc1_age

## --- End of File


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module ParForTest

using ParallelAccelerator
#ParallelAccelerator.DomainIR.set_debug_level(3)
#ParallelAccelerator.ParallelIR.set_debug_level(3)
#ParallelAccelerator.CGen.set_debug_level(3)
@acc function parfor1(n)
 A = Array(Int, n, n)
 @par for i in 1:n, j in 1:n
    A[i,j] = i * j
 end
 return sum(A)
end

@acc function parfor2(n)
 A = Array(Int, n, n)
 s::Int = 0
 m::Int = 0
 @par s(+) m(+) for i in 1:n, j = 1:n
    A[i,j] = i * j
    s = s + A[i,j]
    m = m + 1
 end
 return s * m
end

@acc function parfor3(n)
 A = Array(Int, n, n)
 s::Array{Int, 1} = zeros(n)
 m::Int = 0
 @par s(.+) m(+) for i in 1:n
    for j = 1:n
      A[j,i] = i * j
    end
    s = s .+ A[:,i]
    m = m + 1
 end
 return s .* m
end

function test1()
  parfor1(10) == @noacc parfor1(10)
end

function test2()
  parfor2(10) == @noacc parfor2(10)
end

function test3()
  parfor3(10) == @noacc parfor3(10)
end

end

using Base.Test
println("Testing parfor support via @par macro...")
@test ParForTest.test1()
@test ParForTest.test2()
@test ParForTest.test3()
println("Done testing parfor.")


text:
module HydrogenDimerMRCC

export DimerParameters
export PolarCoords
export make_mrcc_input_file

using Printf
import Base.@kwdef

struct PolarCoords
	# coordinates that describe the length and orientation of a hydrogen molecule

	bondlength::Float64  # bond length of the hydrogen molecule
	theta::Float64       # space-fixed polar angle of the hydrogen molecule
	phi::Float64         # space-fixed azimuthal angle of the hydrogen molecule

	function PolarCoords(bondlength, theta, phi)
		if bondlength <= 0
			error("The molecule's bond length must be positive. ENTERED: $bondlength")
		end

		if !(0.0 <= theta <= pi)
			error("The polar angle must be between 0 and pi. ENTERED: $theta")
		end

		if !(0.0 <= phi <= 2*pi)
			error("The azimuthal angle must be between 0 and 2*pi. ENTERED: $phi")
		end

		new(bondlength, theta, phi)
	end
end

@kwdef struct DimerParameters
	# coordinates that describe the syste of the entire dimer system

	R::Float64
	polar1::PolarCoords
	polar2::PolarCoords

	function DimerParameters(R, r1, theta1, phi1, r2, theta2, phi2)
		if (R <= 0)
			error("The distance between the two molecules must be entered as a positive number. ENTERED: $R")
		end
		# The other parameters are checked in the PolarCoords constructor
		polar1 = PolarCoords(r1, theta1, phi1)
		polar2 = PolarCoords(r2, theta2, phi2)

		new(R, polar1, polar2)
	end

	function DimerParameters(R, polar1, polar2)
		if (R <= 0)
			error("The distance between the two molecules must be entered as a positive number. ENTERED: $R")
		end
		new(R, polar1, polar2)
	end
end

@kwdef struct MRCCInputFileInfo
	titlecomment::String
	calctype::String
	memoryMB::Int64
	threshold::Int64
	basis::String
	geometry::String
end

struct Cartesian3D
	x::Float64
	y::Float64
	z::Float64
end

function distance(pos1::Cartesian3D, pos2::Cartesian3D)
	dx = pos1.x - pos2.x
	dy = pos1.y - pos2.y
	dz = pos1.z - pos2.z

	sqrt(dx*dx + dy*dy + dz*dz)
end

function coordinate_line(pos::Cartesian3D)
	# returns a string for the cartesian coordinates of a hydrogen atom
	# for an MRCC input file. The comment is for users looking at the input file

	xpos_str = "$(@sprintf("% .9f", pos.x))"
	ypos_str = "$(@sprintf("% .9f", pos.y))"
	zpos_str = "$(@sprintf("% .9f", pos.z))"

	hatom_str = "H    $xpos_str   $ypos_str   $zpos_str\n"
	return hatom_str
end

function atom_positions(hmol::PolarCoords, com::Cartesian3D)
	# find the positions of the atoms relative to the centre of mass (com)
	shift = 0.5 * hmol.bondlength

	dx = shift * sin(hmol.theta) * cos(hmol.phi)
	dy = shift * sin(hmol.theta) * sin(hmol.phi)
	dz = shift * cos(hmol.theta)

	# return the positions of the atoms in this hydrogen molecule
	atom_pos1 = Cartesian3D(com.x + dx, com.y + dy, com.z + dz)
	atom_pos2 = Cartesian3D(com.x - dx, com.y - dy, com.z - dz)

	return atom_pos1, atom_pos2
end

function write_input_file(filename::AbstractString, mrcc::MRCCInputFileInfo)
	# Create the file
	fout = open(filename, "w")

	# Create the contents of the file
	write(fout, "# $(mrcc.titlecomment)\n")
	write(fout, "calc=$(mrcc.calctype)\n")
	write(fout, "mem=$(mrcc.memoryMB)MB\n")
	write(fout, "cctol=$(mrcc.threshold)\n")
	write(fout, "basis=$(mrcc.basis)\n")
	write(fout, "\n")
	write(fout, "unit=angs\n")
	write(fout, "geom=xyz\n")
	write(fout, "4\n")
	write(fout, "$(mrcc.geometry)\n")

	# Close it
	close(fout)
end

function make_mrcc_input_file(dparam::DimerParameters, pathtofilename::AbstractString)
	# Create the polar coordinates and centre of masses of the two molecules
	# The first is at the origin, the second is at a distance R on the z-axis
	mol1 = dparam.polar1
	com1 = Cartesian3D(0.0, 0.0, 0.0)
	mol2 = dparam.polar2
	com2 = Cartesian3D(0.0, 0.0, dparam.R)

	# Get the 3D cartesian coordinates for the four atoms
	atomA, atomB = atom_positions(mol1, com1)
	atomC, atomD = atom_positions(mol2, com2)

	# Get the strings for the cartesian coordinates for the MRCC file
	# They will be used to create the 'geometry block' of the MRCC file
	lines = [coordinate_line(atom) for atom in [atomA, atomB, atomC, atomD]]
	geometryblock = join(lines)

	# Collect the information needed for the MRCC Input File
	# Most of this information is fixed for this set of calculations
	mrcc = MRCCInputFileInfo(titlecomment = "hydrogen dimer",
	                         calctype = "ccsd(t)",
	                         memoryMB = 8192,
	                         threshold = 10,
	                         basis = "aug-cc-pVQZ",
	                         geometry = geometryblock)

	# Create the MRCC input file
	write_input_file(pathtofilename, mrcc)
end

end # module


text:
abstract SVMSpec{TKernel <: Kernel}

typealias LinearSVMSpec{TKernel<:ScalarProductKernel} SVMSpec{TKernel}

abstract SVMRiskSpec{TKernel<:Kernel, TLoss<:Loss} <: SVMSpec{TKernel}

typealias SVCSpec{TKernel<:Kernel, TLoss<:MarginBasedLoss} SVMRiskSpec{TKernel, TLoss}
typealias SVRSpec{TKernel<:Kernel, TLoss<:DistanceBasedLoss} SVMRiskSpec{TKernel, TLoss}

abstract SVMRegRiskSpec{TKernel<:Kernel, TLoss<:Loss, TReg<:Penalty} <: SVMRiskSpec{TKernel, TLoss}

"""
`CSVM <: SVMRegRiskSpec <: SVMRiskSpec <: SVMSpec`

Description
============

Specification for a soft-margin Support Vector Machine using the C formulation.

Fields
=======

- **`.kernel`** : The utilized kernel function.

- **`.loss`** : The utilized loss function.

- **`.reg`** : The utilized regularizer, including penality parameter.
Note that the stored parameter is usually not C, but λ (= C⁻¹)

- **`.C`** : The penalty parameter. It plays a role similar to λ⁻¹ in Empirical Risk Minimization.
Conceptually, it denotes the trade off between model complexity and the training error.
Larger C will increase the penalization of training errors and thus lead to behavior
more similar to a hard margin classifier.

See also
=========

`PrimalSVM`, `DualSVM`

"""
immutable CSVM{TKernel<:Kernel, TLoss<:Loss, TReg<:Penalty} <: SVMRegRiskSpec{TKernel, TLoss, TReg}
  kernel::TKernel
  loss::TLoss
  reg::TReg
  C::Float64

  function CSVM(kernel::TKernel, loss::TLoss, reg::TReg)
    new(kernel, loss, reg, Float64(1 / reg.λ))
  end
end

function CSVM(kernel::Kernel, loss::Loss, reg::Penalty)
  CSVM{typeof(kernel), typeof(loss), typeof(reg)}(kernel, loss, reg)
end

function CSVM{TReg<:Penalty}(;
              kernel::Kernel = ScalarProductKernel(),
              loss::Loss = L2HingeLoss(),
              regtype::Type{TReg} = L2Penalty,
              C::Real = 1)
  reg = TReg(Float64(1 / C))
  CSVM{typeof(kernel), typeof(loss), typeof(reg)}(kernel, loss, reg)
end

# ==========================================================================
# Base.show, Base.print

function print(io::IO, model::CSVM)
  print(io, "C-SVM ($(typeof(model.kernel).name.name), $(typeof(model.loss).name.name), $(typeof(model.reg).name.name), C = $(model.C))")
end

function show(io::IO, model::CSVM)
  println(io, typeof(model))
  _printvariable(io, 9, ".kernel", model.kernel)
  _printvariable(io, 9, ".loss", model.loss)
  _printvariable(io, 9, ".reg", model.reg)
  _printvariable(io, 9, ".C", model.C, newline = false)
end


text:
using MotionPlanningProject
using Test

@testset "MotionPlanningProject.jl" begin
    # Write your tests here.
end


text:
@doc """OKUMURA-HATA URBAN model
 #  freq: signal frequency(500Mhz e 1500Mhz);
 #  tyArea: area type (1-rural, 2-suburban e 3-urban);
 #  areaKind : cyte type (1-small, 2-medium e 3-large);
 #  hb: base station's height;
 #  hm: mobile's height;
"""
type OkumuraHataModel<:Model
  freq::Float64
  txH::Float64
  rxH::Float64
  cityKind::Int
  areaKind::Int
  checkFreqRange::Bool
  OkumuraHataModel() = new(800,50,1.5,CityKind.Medium,AreaKind.Urban,true)
end

@doc "Calculate the path loss with the Okumura-Hata Model"
function pathloss_op(model::OkumuraHataModel,dist::AbstractFloat)
    # frequency
   f::Float64  = model.freq;
   # Checking the frequency range
   if (model.checkFreqRange)
     if (f < 500 || f > 1500)
      error("The frequency range for OkumuraHata is 500MHz-1500Mhz")
     end
   end
   #Bts and MS Heights
   Hm::Float64 = model.rxH
   Hb::Float64 = model.txH

   cityKind::Int = model.cityKind
   areaKind::Int = model.areaKind
          #a(hm) from  hata model, hm = altura móvel
          a::Float64=0.0;
          if (f <= 200 && cityKind==CityKind.Large)
                # Large cities and f<=200 Mhz
                a = 8.29*(log10(1.54*Hm))^2- 1.1;

          elseif (f>=400 && cityKind==CityKind.Large)
               #Large cities and f>= 400 MHz
                a = 3.2*((Math.log10(11.75*Hm)^2))- 4.97;
          else
              #a(hm) for small and medium cities, and large cities where f<200Mhz and f>400Mhz
               a = (1.1*log10(f-0.7))*Hm - (1.56*log10(f-0.8));
          end
         # Path loos for urban area
           lossUrban::Float64 = 69.55 + (26.16)*log10(f)-13.82*log10(Hb) - a + (44.9-6.55*log10(Hb))*log10(dist);
          if (areaKind== AreaKind.Open)
            lossOpen::Float64 = lossUrban - 4.78 *((log10(f))^2)+18.33*log10(f)-40.94;
            return lossOpen;
          elseif (areaKind==AreaKind.SubUrban)
              #Loss for open are
              lossSubUrban::Float64 =  lossUrban  - 2*(log10(f/28.0))^2 - 5.4# //#Loss for suburban area
              return lossSubUrban;
          end
          return lossUrban;
end


text:
# This is where we do some type pyracy on the PSY types
get_n_buses(sys::PSY.System) = length(sys.bus_numbers)

function _filter_function(x::T) where {T <: PSY.StaticInjection}
    if PSY.get_dynamic_injector(x) === nothing
        return false
    end

    if hasfield(T, :status)
        return PSY.get_status(x) * PSY.get_available(x)
    else
        return PSY.get_available(x)
    end
    error("Filtering of $(PSY.get_name(x)) failed")
end

function get_injectors_with_dynamics(sys::PSY.System)
    return PSY.get_components(PSY.StaticInjection, sys, x -> _filter_function(x))
end

function get_injection_without_dynamics(sys::PSY.System)
    return PSY.get_components(
        PSY.StaticInjection,
        sys,
        x ->
            PSY.get_dynamic_injector(x) === nothing &&
                PSY.get_available(x) &&
                !isa(x, PSY.ElectricLoad),
    )
end

function get_dynamic_branches(sys::PSY.System)
    return PSY.get_components(PSY.DynamicBranch, sys, x -> PSY.get_available(x))
end

function _transform_all_lines!(sys::PSY.System)
    for br in PSY.get_components(PSY.DynamicBranch, sys)
        dyn_br = DynamicBranch(br)
        @debug "Converted $(PSY.get_name(dyn_br)) to DynamicBranch"
        add_component!(sys, dyn_br)
    end
end

function transform_ybus_to_rectangular(
    ybus::SparseArrays.SparseMatrixCSC{Complex{Float64}, Int},
)
    # TODO: Improve performance here
    return hcat(vcat(real(ybus), -imag(ybus)), vcat(imag(ybus), real(ybus)))
end

function transform_branches_to_dynamic(sys::PSY.System, ::Type{T}) where {T <: PSY.ACBranch}
    for b in PSY.get_components(T, sys)
        dyn_branch = PSY.DynamicBranch(b)
        PSY.add_component!(sys, dyn_branch)
    end
    return
end


text:
verbose = true
 using Dates
 if verbose println(now()) end
 if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 using RvSpectMLBase
 using EchelleInstruments#, EchelleInstruments.NEID
 using RvSpectML
 if verbose println("# Loading NeidSolarScripts")    end
 using NeidSolarScripts
 println("# Loading other packages 1/2")
 using ArgParse

println("# Parsing arguments...")
   lsf_width = 3.0e3
   function parse_commandline()
     s = ArgParseSettings( description = "Calculate order CCFs from NEID L1 FITS files.")
     #import_settings!(s, s_files_only, args_only=false)
     add_arg_group!(s, "Files", :argg_files)
     @add_arg_table! s begin
         "manifest"
             help = "Manifest file (CVS) containing FITS files to analyze.\nExpects columns named Filename, bjd, target, and continuum_filename."
             arg_type = String
             default = "manifest.csv"
             #required = true
         "output"
             help = "Filename for output CCFs (jld2)"
             arg_type = String
             default = "daily_ccfs.jld2"
         "--param"
             help = "Parameters file"
             arg_type = String
         "--overwrite"
            help = "Specify it's ok to overwrite the output file."
            action = :store_true
      end
      add_arg_group!(s, "CCF parameters", :argg_ccf_param)
      @add_arg_table! s begin
        "--orders_to_use"
           help = "First and last order _index_ to compute CCFs for"
           nargs = 2
           arg_type = Int64
           #default = [ min_order(NEID2D()), max_order(NEID2D()) ]
           #default = [ first(orders_to_use_default(NEID2D())), last(orders_to_use_default(NEID2D())) ]
        "--mask_scale_factor"
            help = "Specify CCF mask width scale as multiple of NEID default v width " * string(default_ccf_mask_v_width(NEID2D()))
            arg_type = Float64
            #default = round(lsf_width/default_ccf_mask_v_width(NEID2D()), sigdigits=3)
        "--line_width_50_default"
            help = "Specify default line full width half maximum"
            arg_type = Float64
            #default = 7.9e3 # m/2
        "--variable_mask_scale"
            help = "Vary width of mask to compensate for solar rotation."
            action = :store_true
        "--ccf_mid_vel"
            help = "Middle of velocity range to use for CCFs."
            arg_type = Float64
            #default = 0.0
        "--range_no_mask_change"
            help = "Avoid lines that would result in a mask change due to BC using line width times this factor (TODO: Verify definition)"
            arg_type = Float64
            #default = 5.0
        "--v_step"
            help = "Specify v step for computing CCF "
            arg_type = Float64
            #default = 155.0 # m/2
      end
      add_arg_group!(s, "Line list parameters", :argg_line_list_param)
      @add_arg_table! s begin
         "--line_list_filename"
             help = "Line list file"
             arg_type = String
             #default = joinpath(pkgdir(NeidSolarScripts),"data","solar_line_list_espresso.csv")
         "--recompute_line_weights"
             help = "Force recalculation of SNR-based line weight factors."
             action = :store_true
      end
      add_arg_group!(s, "Continuum normalization parameters", :argg_continuum_param)
      @add_arg_table! s begin
         "--sed_filename"
             help = "Filename for input SED to normalize by"
             arg_type = String
             default = "/home/eford/Code/RvSpectMLEcoSystem/NeidSolarScripts/data/neidMaster_HR_SmoothLampSED_20210101.fits"
         "--continuum_poly_half_width"
             help = "Half width for window to use for smoothing prior to polynomial fit to continuum."
             arg_type = Int64
             default = 50
         "--quantile_fit_continuum"
             help = "Quantile of rolling window to use prior to polynomial fit to continuum."
             arg_type = Float64
             default = 0.9
         "--order_poly_continuum"
             help = "Order of polynomial to fit after dividing by SED proivded."
             arg_type = Int64
             default = 4
        "--anchors_filename"
             help = "Filename for anchor locations to use in computing continuum to normalize by."
             arg_type = String
             #default = "/home/eford/Code/RvSpectMLEcoSystem/NeidSolarScripts/data/neidMaster_HR_SmoothLampSED_20210101.fits"
         "--smoothing_half_width"
             help = "Half width for window to use for smoothing prior to findind local maximum."
             arg_type = Int64
             default = 6
         "--stretch_factor"
             help = "Stretch factor to use for scaling flux in distance calculation for rolling pin continuum normalization"
             arg_type = Float64
             default = 5.0
         "--merging_threshold"
             help = "Threshold (in λ units) for merging nearby continuum anchors."
             arg_type = Float64
             default = 0.25
         "--fwhm_continuum"
            help = "Full width half-max to use for performing continuum normalization."
            arg_type = Float64
            default = Continuum.fwhm_sol/1000
         "--min_rollingpin_r"
            help = "Minimum rolling pin radius for continuum normalization in uints of FWHM."
            arg_type = Float64
            default = 100.0
         "--nu_continuum"
            help = "Exponent for rolling pin radius in continuum normalization"
            arg_type = Float64
            default = 1.0
         "--apply_continuum_normalization"
            help = "Apply continuum normalization."
            action = :store_true
        "--continuum_normalization_individually"
            help = "Calculate continuum normalization for each file rather than averaged spectra."
            action = :store_true
        "--save_continuum_normalization"
            help = "NOT IMPLEMENTED.  Save continuum normalization to file."
            action = :store_true
        "--recompute_continuum_normalization"
            help = "CURRENTLY ONLY OPTION.  Recompute continuum normalization even if continuum file exists."
            action = :store_true
      end
      add_arg_group!(s, "Filter manifest for ", :argg_filter_param)
      @add_arg_table! s begin
         "--target"
            help = "Target field"
            arg_type = String
            default = "Sun"
         "--datestr"
            help = "Filenames that contain date string."
            arg_type = String
        "--max_solar_hour_angle"
           help = "Maximum absolute value of solar hour angle."
           arg_type = Float64
           default = 6.0
         "--max_airmass"
            help = "Maximum airmass."
            arg_type = Float64
            #default = 10.0
         "--min_snr_factor"
            help = "Minimum SNR relative to max_snr."
            arg_type = Float64
            #default = 0.0
         "--max_snr"
            help = "Specify max_snr manually."
            arg_type = Float64
         "--start_time"
            help = "Specify daily start time for CCFs to be processed (HH MM)"
            nargs = 2
            arg_type = Int64
            default = [0, 0]
            #default = [18, 30]
        "--stop_time"
           help = "Specify daily stop time for CCFs to be processed (HH MM)"
           nargs = 2
           arg_type = Int64
           #default = [ min_order(NEID2D()), max_order(NEID2D()) ]
           #default = [22, 30]
           default = [23, 59]
         "--max_num_spectra"
            help = "Maximum number of spectra to process."
            arg_type = Int
            default = 250   # TODO: Increase  after finish testing
      end

     return parse_args(s)
 end
 args = parse_commandline()

if verbose   println("# Loading other packages 2/2")    end
 using CSV, DataFrames, Query, Dates
 using JLD2, FileIO, SHA
 using StatsBase, Statistics, NaNMath

 # Filename arguments
 @assert isfile(args["manifest"]) || islink(args["manifest"])
 manifest_filename = args["manifest"]
 #args["overwrite"] = true   # TODO: Comment after done testing
 if isfile(args["output"]) && !args["overwrite"]
    error("# Output file " * args["output"] * " already exists (size = " * string(filesize(args["output"])) * " ).")
 end
 @assert !isfile(args["output"]) || args["overwrite"] == true
 @assert match(r"\.jld2$",args["output"]) != nothing
 daily_ccf_filename = args["output"]
 touch(daily_ccf_filename)

 file_hashes = Dict{String,String}()
 start_processing_time = now()

 # Include parameters file if specified
 if args["param"] != nothing
     @warn "Haven't implemented parameters file for this script yet."
    if isfile(args["param"]) || islink(args["param"])
        println("# Including parameter code from ", args["param"], ".")
        include(args["param"])
        file_hashes["param.in"] = bytes2hex(args["param"])
    end
 end
 # Overide values from param file if on commandline.  If missing, set defaults
 if args["line_list_filename"] != nothing
    line_list_filename = args["line_list_filename"]
 elseif !@isdefined line_list_filename
     #line_list_filename = joinpath(pkgdir(NeidSolarScripts),"data","solar_line_list_espresso.csv")
     line_list_filename = joinpath(pkgdir(NeidSolarScripts),"data","espresso+neid_mask_97_to_108.mas")
 end
 @assert isfile(line_list_filename) || islink(line_list_filename)
 if args["sed_filename"] != nothing
    sed_filename = args["sed_filename"]
 #elseif !@isdefined sed_filename
 #     sed_filename = joinpath("/home/eford/Code/RvSpectMLEcoSystem/NeidSolarScripts/data","neidMaster_HR_SmoothLampSED_20210101.fits")
 end
 @assert !@isdefined(sed_filename) || isfile(sed_filename) || islink(sed_filename)
 @assert !@isdefined(anchors_filename) || isfile(anchors_filename) || islink(anchors_filename)

 # Extract arguments with non-standard types
 if args["orders_to_use"] != nothing && length(args["orders_to_use"]) == 2
     orders_to_use = first(args["orders_to_use"]):last(args["orders_to_use"])
 elseif !@isdefined orders_to_use
     orders_to_use = first(orders_to_use_default(NEID2D())):last(orders_to_use_default(NEID2D()))
 end
 @assert min_order(NEID2D()) <= first(orders_to_use) <= max_order(NEID2D())
 @assert min_order(NEID2D()) <= last(orders_to_use) <= max_order(NEID2D())

 if args["ccf_mid_vel"] != nothing
     ccf_mid_velocity = args["ccf_mid_vel"]
 elseif !@isdefined ccf_mid_velocity
     ccf_mid_velocity = 0.0
 end
 @assert -5e5 < ccf_mid_velocity < 5e5  # arbitrary ball park restrictions to catch obvious typeos

 if args["v_step"] != nothing
     v_step = args["v_step"]
 elseif !@isdefined v_step
     v_step = 150.0
 end
 @assert 50 <= v_step <= 2000 # m/s

 if args["mask_scale_factor"] != nothing
     mask_scale_factor = args["mask_scale_factor"]
 elseif !@isdefined mask_scale_factor
     mask_scale_factor = round(lsf_width/default_ccf_mask_v_width(NEID2D()), sigdigits=3)
 end
 @assert 1 <= mask_scale_factor <= 64 # arbitrary big number to catch obvious typeos

 if args["line_width_50_default"] != nothing
     line_width_50_default = args["line_width_50_default"]
 elseif !@isdefined line_width_50_default
     line_width_50_default = 7.9e3 # m/s
 end
 @assert 2e3 <= line_width_50_default <= 50e3   # arbitrary big number to catch obvious typeos

 if args["range_no_mask_change"] != nothing
     range_no_mask_change = args["range_no_mask_change"]
 elseif !@isdefined range_no_mask_change
     range_no_mask_change = 5.0 # m/s
 end
 @assert 0 <= range_no_mask_change <= 10
 start_time = args["start_time"] != nothing && length(args["start_time"]) == 2 ? Time(args["start_time"][1], args["start_time"][2]) : Time(0,0)
 stop_time =  args["stop_time"] != nothing && length(args["stop_time"]) == 2 ? Time(args["stop_time"][1],  args["stop_time"][2]) : Time(12,59,59)

 @assert 3 <= args["continuum_poly_half_width"] <= 200 # arbitrary limits for now
 @assert 0.5 <= args["quantile_fit_continuum"] <= 1-1/(2*args["smoothing_half_width"])
 @assert 0 <= args["order_poly_continuum"] <= 6 # arbitrary limits for now
 @assert 3 <= args["smoothing_half_width"] <= 200 # arbitrary limits for now
 @assert 1.0 <= args["stretch_factor"] <= 100    # arbitrary limits for now
 @assert 0.1 <= args["merging_threshold"] <= 2.0 # arbitrary limits for now
 @assert 1.0 <= args["fwhm_continuum"] <= 30.0   # km/s  arbitrary limits for now
 @assert 10 <= args["min_rollingpin_r"] <= 1000 # arbitrary limits for now
 @assert 0.7 <= args["nu_continuum"] <= 1.3  # Recommendations from Cretignier et al.
 #args["apply_continuum_normalization"] = true   # TODO Remove after done testing
 #args["continuum_normalization_individually"] = false #  TODO Remove after done testing
 #args["recompute_line_weights"] = true  # TODO Remove after done testing

if verbose println("# Reading manifest of files to process.")  end
  df_files  = CSV.read(manifest_filename, DataFrame)
    @assert size(df_files,1) >= 1
    @assert hasproperty(df_files,:Filename)
    @assert hasproperty(df_files,:target)
    @assert hasproperty(df_files,:bjd)
    @assert hasproperty(df_files,:ssbz)
    @assert hasproperty(df_files,:exptime)
    @assert hasproperty(df_files,:airmass)
    if args["target"] == "Sun" || args["target"] == "Solar"
      @assert hasproperty(df_files,:alt_sun)  # TODO: Compute if not avaliable?
      @assert hasproperty(df_files,:Δfwhm²)   # TODO: Compute if not avaliable?
      df_files[!,"solar_hour_angle"] = NeidSolarScripts.SolarRotation.get_solar_hour_angle(df_files.bjd,obs=:WIYN)
      if eltype(df_files[!,:order_snrs]) == String   # TODO: Compute if not avaliable?
    end
    @assert hasproperty(df_files,:order_snrs)
        df_files[!,:order_snrs] = map(i->parse.(Float64,split(df_files[i,:order_snrs][2:end-1],',')),1:size(df_files,1))
    end
    @assert eltype(df_files[!,:order_snrs]) == Vector{Float64}

  if  args["max_snr"] != nothing
      max_snr = args["max_snr"]
  elseif !@isdefined max_snr
      df_files_tmp = df_files |>
        @filter( args["target"] == nothing || _.target == args["target"] ) |> DataFrame
      max_snr = maximum(NaNMath.sum.(df_files_tmp.order_snrs))
  end
  @assert 0 < max_snr < Inf

  @assert args["max_airmass"] == nothing || 1 < args["max_airmass"] <= 10  # not tested beyond 10
  @assert args["min_snr_factor"] == nothing || 0 <= args["min_snr_factor"] < 1
  @assert args["max_solar_hour_angle"] == nothing || 0 < args["max_solar_hour_angle"] <= 6

  df_files_use = df_files |>
    @filter( args["target"] == nothing || _.target == args["target"] ) |>
    @filter( args["datestr"] == nothing || occursin(args["datestr"],_.target) ) |>
    @filter( args["max_airmass"] == nothing || _.airmass <= args["max_airmass"] ) |>
    @filter( args["max_solar_hour_angle"] == nothing || abs(_.solar_hour_angle) <= args["max_solar_hour_angle"] ) |>
    @filter( args["start_time"] == nothing || Time(julian2datetime(_.bjd)) >= start_time ) |>
    @filter( args["stop_time"] == nothing || Time(julian2datetime(_.bjd)) <= stop_time ) |> # TODO for other instruments may need to deal wtih cross end of 24 UTC
    @filter( args["min_snr_factor"] == nothing || sum(_.order_snrs) >= args["min_snr_factor"] * max_snr ) |>
    @take(args["max_num_spectra"] ) |> @orderby(_.bjd) |>
    DataFrame
  println("# Found ", size(df_files_use,1), " files of ",  size(df_files,1), " to process.")
  @assert size(df_files_use,1) >= 1

df_files_cleanest = df_files_use |>
    @filter( _.airmass <= 2.0 ) |>
    @filter( abs(_.solar_hour_angle) <= 1 ) |>
    @filter( sum(_.order_snrs) >= 0.9 * max_snr ) |>
    @take(args["max_num_spectra"] ) |>
    DataFrame
  println("# Found ", size(df_files_cleanest,1), " files considered clean.")
  #@assert size(df_files_cleanest,1) >= 1

if verbose println(now()) end

if @isdefined(sed_filename)
   sed = Continuum.read_master_sed_neid(filename=sed_filename)
   min_order_for_continuum = min_order(NEID2D()) # 12
   max_order_for_continuum = max_order(NEID2D())
   orders_to_use_for_continuum = min_order_for_continuum:max_order_for_continuum
end

pipeline_plan = PipelinePlan()
 if verbose println("# Reading data files.")  end
 # For v1.0.0
 #all_spectra = Spectra2DBasic{Float64, Float32, Float32, Matrix{Float64}, Matrix{Float32}, Matrix{Float32}, NEID2D}[]
 # For v1.1.*
 all_spectra = Spectra2DBasic{Float64, Float64, Float32, Matrix{Float64}, Matrix{Float64}, Matrix{Float32}, NEID2D}[]
 spec = NEID.read_data(first(eachrow(df_files_use)).Filename)
 mean_clean_flux = zeros(Float64,size(spec.flux));
 mean_clean_var = zeros(Float64,size(spec.flux));
 #mean_clean_weight_sum = zeros(Float64,size(spec.flux));
 mean_clean_flux_sed_normalized = zeros(Float64,size(spec.flux));
 mean_clean_var_sed_normalized = zeros(Float64,size(spec.flux));
 #mean_clean_flux_sed_normalized_weight_sum = zeros(Float64,size(spec.flux));
 mean_clean_flux_continuum_normalized = zeros(Float64,size(spec.flux));
 mean_clean_var_continuum_normalized = zeros(Float64,size(spec.flux));
 #mean_clean_flux_continuum_normalized_weight_sum = zeros(Float64,size(spec.flux));
 normalization_anchors_list = [];
 for (i,row) in enumerate(eachrow(df_files_use))
    if !(isfile(row.Filename)||islink(row.Filename))
        println("# Couldn't find input FITS file: ", row.Filename)
        continue
    end
    #= if !(isfile(row.continuum_filename)||islink(row.continuum_filename))
        println("# Couldn't find matching continuum file: ", row.continuum_filename)
        continue
    end =#
    if verbose
        println("# Reading ", row.Filename, "(",i,"/",size(df_files_use,1),")")
        flush(stdout)
    end
    local spec = NEID.read_data(row)
    file_hashes[row.Filename] = bytes2hex(sha256(row.Filename))
    weights = 1
    if row.Filename ∈ df_files_cleanest.Filename
            mean_clean_flux .+= spec.flux # .*weights
            mean_clean_var .+= spec.var # .*weights
    end

    #=
    m = match(r"(neidL1_\d+[T_]\d+)\.fits$", row.Filename)
    continuum_filename = joinpath(neid_data_path,target_subdir,output_dir,"continuum", m[1] * "_continuum=new.jld2")
    if !(isfile(continuum_filename)||islink(continuum_filename))
        println("# Couldn't find matching continuum file", continuum_filename)
        continue
    end
    continuum = load(continuum_filename, "continuum")
    =#

    if @isdefined sed
      (f_norm, var_norm) = Continuum.normalize_by_sed(spec.λ,spec.flux,spec.var, sed; poly_order = args["order_poly_continuum"], half_width = args["continuum_poly_half_width"], quantile = args["quantile_fit_continuum"], orders_to_use=orders_to_use_for_continuum)
      if row.Filename ∈ df_files_cleanest.Filename
            mean_clean_flux_sed_normalized .+= f_norm # .*weights
            mean_clean_var_sed_normalized .+= var_norm # .*weights
      end
    else
        f_norm = spec.flux
        var_norm = spec.var
    end

    if args["apply_continuum_normalization"]==true && args["continuum_normalization_individually"]==true
        local anchors, continuum, f_filtered
        if args["anchors_filename"] != nothing
            @assert isfile(args["anchors_filename"]) && filesize(args["anchors_filename"])>0
            anchors = load(args["anchors_filename"],"anchors")
            (anchors, continuum, f_filtered) = Continuum.calc_continuum(spec.λ, f_norm, var_norm, anchors;
                orders_to_use = orders_to_use_for_continuum, verbose = true )
            push!(normalization_anchors_list, anchors)
        else
            (anchors, continuum, f_filtered) = Continuum.calc_continuum(spec.λ, f_norm, var_norm; fwhm = args["fwhm_continuum"]*1000, ν = args["nu_continuum"],
                stretch_factor = args["stretch_factor"], merging_threshold = args["merging_threshold"], smoothing_half_width = args["smoothing_half_width"], min_R_factor = args["min_rollingpin_r"],
                orders_to_use = orders_to_use_for_continuum, verbose = false )
            push!(normalization_anchors_list, anchors)
        end
            #=
            continuum = load(row.continuum_filename, "continuum")
            file_hashes[row.continuum_filename] = bytes2hex(sha256(row.continuum_filename))
            =#
        f_norm ./= continuum
        var_norm ./= continuum.^2
        if row.Filename ∈ df_files_cleanest.Filename
                mean_clean_flux_continuum_normalized .+= f_norm # .*weights
                mean_clean_var_continuum_normalized .+= var_norm # .*weights
#                mean_clean_var_continuum_normalized_weight_sum .+= weights
            end
        spec.flux .= f_norm
        spec.var .= var_norm

    end
    push!(all_spectra,spec)
 end
 GC.gc()
 dont_need_to!(pipeline_plan,:read_spectra);

 #mean_clean_flux ./= mean_clean_flux_weight_sum
 #mean_clean_flux_sed_normalized ./= mean_clean_flux_continuum_normalized_weight_sum
 #mean_clean_flux_continuum_normalized ./= mean_clean_flux_continuum_normalized_weight_sum

 if args["apply_continuum_normalization"]==true && !(args["continuum_normalization_individually"] == true)
     println("# Computing continuum normalization from mean spectra.")
     local anchors, continuum, f_filtered
     if args["anchors_filename"] !=nothing
         @assert isfile(args["anchors_filename"]) && filesize(args["anchors_filename"])>0
         anchors = load(args["anchors_filename"],"anchors")
     else
         if @isdefined sed
             (anchors, continuum, f_filtered) = Continuum.calc_continuum(spec.λ, mean_clean_flux_sed_normalized, mean_clean_var_sed_normalized; fwhm = args["fwhm_continuum"]*1000, ν = args["nu_continuum"],
                stretch_factor = args["stretch_factor"], merging_threshold = args["merging_threshold"], smoothing_half_width = args["smoothing_half_width"], min_R_factor = args["min_rollingpin_r"],
                orders_to_use = orders_to_use_for_continuum, verbose = false )
        else
            (anchors, continuum, f_filtered) = Continuum.calc_continuum(spec.λ, mean_clean_flux, mean_clean_var; fwhm = args["fwhm_continuum"]*1000, ν = args["nu_continuum"],
                stretch_factor = args["stretch_factor"], merging_threshold = args["merging_threshold"], smoothing_half_width = args["smoothing_half_width"], min_R_factor = args["min_rollingpin_r"],
                orders_to_use = orders_to_use_for_continuum, verbose = false )
        end
    end
    normalization_anchors_list = anchors

    for (i,row) in enumerate(eachrow(df_files_use))
        (anchors, continuum, f_filtered) = Continuum.calc_continuum(all_spectra[i].λ, all_spectra[i].flux, all_spectra[i].var,
                    anchors, smoothing_half_width = args["smoothing_half_width"], orders_to_use=orders_to_use_for_continuum)
        all_spectra[i].flux ./= continuum
        all_spectra[i].var ./= continuum.^2
        if row.Filename ∈ df_files_cleanest.Filename
            mean_clean_flux_continuum_normalized .+= all_spectra[i].flux # .*weights
            mean_clean_var_continuum_normalized .+= all_spectra[i].var # .*weights
            #                mean_clean_var_continuum_normalized_weight_sum .+= weights
        end
    end
 end

 order_list_timeseries = extract_orders(all_spectra, pipeline_plan,  orders_to_use=orders_to_use, remove_bad_chunks=false, recalc=true )

line_width = line_width_50_default
 max_mask_scale_factor = max(mask_scale_factor, max(lsf_width,
                             line_width/sqrt(8*log(2)))/default_ccf_mask_v_width(NEID2D())) # 4.0
 max_bc = RvSpectMLBase.max_bc
 #max_bc = RvSpectMLBase.max_bc_earth_rotation
 #max_orders = min_order(NEID2D()):max_order(NEID2D())
 #good_orders = orders_to_use_default(NEID2D())
 #orders_to_use = max_orders
 if isfile(line_list_filename) && (filesize(line_list_filename)>0) && !args["recompute_line_weights"]
   println("# Reading ", line_list_filename)
   line_list_espresso = CSV.read(line_list_filename, DataFrame)
   @assert all(map(k->k ∈ names(line_list_espresso), ["lambda","weight","order"]))
   dont_need_to!(pipeline_plan,:clean_line_list_tellurics)
 else
    println("# Can't find ", line_list_filename, ".  Trying ESPRESSO line list.")
    #orders_to_use = good_orders
    #order_list_timeseries = extract_orders(all_spectra,pipeline_plan, orders_to_use=orders_to_use, recalc=true )
    touch(line_list_filename)
    #line_list_filename = "G2.espresso.mas"
    linelist_for_ccf_fn_w_path = joinpath(pkgdir(EchelleCCFs),"data","masks",line_list_filename)
    line_list_espresso = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan, v_center_to_avoid_tellurics=ccf_mid_velocity,
       Δv_to_avoid_tellurics = 2*max_bc+range_no_mask_change*line_width_50_default+max_mask_scale_factor*default_ccf_mask_v_width(NEID2D()), orders_to_use=#=orders_to_use=#56:108, recalc=true, verbose=true)
    #CSV.write(line_list_filename, line_list_espresso)
 end
 file_hashes[line_list_filename] = bytes2hex(sha256(line_list_filename))
 #outputs["line_list_espresso"] = line_list_espresso

if args["variable_mask_scale"]
     maxΔfwhm² = -0.569375
     @assert maximum(df_files_use.Δfwhm²) < maxΔfwhm²
     Δfwhm = 1000.0 .* sqrt.(maxΔfwhm².-df_files_use.Δfwhm²[1:length(all_spectra)])  # How much to increase fwhm by to acheive uniform fwhm
else
    Δfwhm = zeros(0)
end

if verbose println(now()) end
line_width_50 = line_width_50_default
  #order_list_timeseries = extract_orders(all_spectra,pipeline_plan,  orders_to_use=orders_to_use,  remove_bad_chunks=false, recalc=true )
  @time (order_ccfs, order_ccf_vars, v_grid_order_ccfs) = ccf_orders(order_list_timeseries, line_list_espresso, pipeline_plan,
    mask_type=:gaussian, Δfwhm=Δfwhm,
    mask_scale_factor=mask_scale_factor, range_no_mask_change=5*line_width_50, ccf_mid_velocity=ccf_mid_velocity, v_step=v_step,
    v_max=max(range_no_mask_change*line_width_50,2*max_bc), orders_to_use=orders_to_use, allow_nans=true, calc_ccf_var=true,
    recalc=true)

#orders_to_use2 = orders_to_use[map(i->!iszero(order_ccfs[:,i,:]),1:length(orders_to_use))]
#order_list_timeseries2 = extract_orders(all_spectra,pipeline_plan,  orders_to_use=orders_to_use2,  remove_bad_chunks=false, recalc=true )

#=
ccf_dir = joinpath(neid_data_path,target_subdir,"ccfs")
if !isdir(ccf_dir)
   mkdir(ccf_dir)
end
=#
if verbose println(now()) end
#daily_ccf_filename = joinpath(neid_data_path,target_subdir,"output","ccfs", date_str * "_ccfs=default.jld2")
#daily_ccf_filename = joinpath(neid_data_path,target_subdir,"ccfs", date_str * "_ccfs=default.jld2")
println("# Saving results to ", daily_ccf_filename, ".")
  stop_processing_time = now()
  jldopen(daily_ccf_filename, "w") do f
    f["v_grid"] = v_grid_order_ccfs
    f["order_ccfs"] = order_ccfs
    f["order_ccf_vars"] = order_ccf_vars
    f["orders_to_use"] = orders_to_use
    f["manifest"] = df_files_use
    f["calc_order_ccf_args"] = args
    f["ccf_line_list"] = line_list_espresso
    if size(df_files_cleanest,1) >= 1
      f["mean_clean_flux"] = mean_clean_flux
      f["mean_clean_var"] = mean_clean_var
      if size(mean_clean_flux_sed_normalized,1) > 0
          f["mean_clean_flux_sed_normalized"] = mean_clean_flux_sed_normalized
          f["mean_clean_var_sed_normalized"] = mean_clean_var_sed_normalized
      end
      if size(mean_clean_flux_continuum_normalized,1) > 0
          f["mean_clean_flux_continuum_normalized"] = mean_clean_flux_continuum_normalized
          f["mean_clean_var_continuum_normalized"] = mean_clean_var_continuum_normalized
      end
    end
    if args["apply_continuum_normalization"]
        f["normalization_anchors"] = normalization_anchors_list
    end
    f["start_processing_time"] = start_processing_time
    f["stop_processing_time"] = stop_processing_time
    f["daily_ccf_filename"] = daily_ccf_filename
    f["file_hashes"] = file_hashes
  end
if verbose println(now()) end

#=
for (i,row) in enumerate(eachrow(df_files_use))
    m = match(r"(neidL1_\d+[T_]\d+)\.fits$", row.Filename)
    #ccf_filename = joinpath(neid_data_path,target_subdir,"output","ccfs", m[1] * "_ccfs=default.jld2")
    ccf_filename = joinpath(neid_data_path,target_subdir,"ccfs", m[1] * "_ccfs=default.jld2")
    jldopen(ccf_filename, "w") do f
        f["v_grid"] = v_grid_order_ccfs
        f["order_ccfs"] = order_ccfs[:,:,i]
        f["order_ccf_vars"] = order_ccf_vars[:,:,i]
        f["orders_to_use"] = orders_to_use
        for k in keys(row)
            f[string(k)] = row[k]
        end
        f["ccf_filename"] = ccf_filename
    end
    if i==size(order_ccfs,3) break end
end
=#

#=

#good_orders = orders_to_use_default(NEID2D())
 #good_orders = 56:111 # orders_to_use_default(NEID2D())
 #good_orders = 57:96 # orders_to_use_default(NEID2D())
 #good_orders = 55 .+ vcat(1:14, 17:24, 27:29, 37:38,40:41)
 #good_orders = 55 .+ collect(1:14)
 #good_orders = 55 .+ collect(17:24)
 #good_orders = 55 .+ vcat(collect(17:24), collect(27:29), collect(37:38),collect(40:41))
 #good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41))
 #good_orders = [55, 56, 57]
# order_list_timeseries = extract_orders(all_spectra,pipeline_plan,  orders_to_use=good_orders, remove_bad_chunks=false, recalc=true )

outputs["times"] = order_list_timeseries.times
 alt_sun = map(i->calc_solar_alt(all_spectra[i].metadata[:bjd]),1:length(all_spectra))
 obs_idx_max_alt_sun = argmax(alt_sun)
 time_max_alt_sun = order_list_timeseries.times[obs_idx_max_alt_sun]
 outputs["time_max_alt_sun"] = time_max_alt_sun

 order_snrs = [RvSpectMLBase.calc_snr(order_list_timeseries[obsid].data[c]) for c in 1:length(order_list_timeseries[1].data),
                    obsid in 1:length(order_list_timeseries) ]
 order_sum_flux = [NaNMath.sum(order_list_timeseries[obsid].data[c].flux) for c in 1:length(order_list_timeseries[1].data),
                                        obsid in 1:length(order_list_timeseries) ]
 order_mean_flux = [NaNMath.mean(order_list_timeseries[obsid].data[c].flux) for c in 1:length(order_list_timeseries[1].data),
                                                                        obsid in 1:length(order_list_timeseries) ]

 mean_order_mean_flux = vec(mean(order_mean_flux,dims=2))

 order_weights =  (order_snrs[:,obs_idx_max_alt_sun]./order_snrs).^2
 #order_weights = (order_sum_flux[:,obs_idx_max_alt_sun]./order_sum_flux)
 outputs["order_snrs"] = order_snrs
 outputs["order_sum_flux"] = order_sum_flux
 outputs["order_mean_flux"] = order_mean_flux
 outputs["order_weights"] = order_weights

#times = order_list_timeseries.times.-time_max_alt_sun
#plot(times,order_weights', label=:none)
=#

#=
for obsid in 1:length(order_list_timeseries)
    for ch in 1:length(order_list_timeseries[1].data)
        order_list_timeseries[obsid].data[ch].flux .*= order_weights[ch,obsid] #./mean_order_mean_flux[ch]
        order_list_timeseries[obsid].data[ch].var .*= (order_weights[ch,obsid]).^2 # ./mean_order_mean_flux[ch]).^2
    end
end
=#

#=for obsid in 1:length(order_list_timeseries)
    for ch in 1:length(order_list_timeseries[1].data)
        order_list_timeseries[obsid].data[ch].flux ./= order_weights[ch,obsid] #./mean_order_mean_flux[ch]
        order_list_timeseries[obsid].data[ch].var ./= (order_weights[ch,obsid]).^2 # ./mean_order_mean_flux[ch]).^2
    end
end
=#

#=
using Plots
 plot(map(obsid->NaNMath.sum(order_list_timeseries[obsid].data[1].flux)/(NaNMath.sum(order_list_timeseries[obsid].data[end].flux)),
                    1:length(all_spectra)),label="Raw")
 #plot!(map(obsid->NaNMath.sum(order_list_timeseries[obsid].data[10].flux)*order_weights[10,obsid]/(NaNMath.sum(order_list_timeseries[obsid].data[30].flux*order_weights[30,obsid])), 1:length(all_spectra)),label="weighted")
 #plot!(map(obsid->NaNMath.sum(order_list_timeseries[obsid].data[10].flux)*order_weights[10,obsid]^2/(NaNMath.sum(order_list_timeseries[obsid].data[30].flux*order_weights[30,obsid]^2)), 1:length(all_spectra)),label="weighted^2")

# plot(order_list_timeseries[obsid].data[ch].λ,order_list_timeseries[obsid].data[ch].flux)
=#

#=
times = order_list_timeseries.times.-time_max_alt_sun
for obsid in 1:length(order_list_timeseries)
    order_snrs[:,obsid] ./= order_snrs[15,obsid]
    order_snrs_new[:,obsid] ./= order_snrs_new[15,obsid]
end

plot(times,order_snrs', label=:none)
plot!(times,order_snrs_new'./order_snrs_new[15,:]', label=:none)

plt = plot()
for obsid in 1:length(order_list_timeseries)
    for ch in 1:10:length(order_list_timeseries[1].data)
        plot(order_list_timeseries[obsid].data[ch].λ,order_list_timeseries[obsid].data[ch].flux)

    end
end
=#

#=
maxΔfwhm² = -0.569375
@assert maximum(df_files_use.Δfwhm²) < maxΔfwhm²
Δfwhm = 1000.0 .* sqrt.(maxΔfwhm².-df_files_use.Δfwhm²)  # How much to increase fwhm by to acheive uniform fwhm

#msf = 1.4; fwtf = 1.5
#msf = 6.0; fwtf = 1.5
#msf = 0.62; fwtf = 1.5  # one pixel because tophat not using NEID's default v width
#msf = 1.0; fwtf = 1.5  # using NEID's default v width
msf = lsf_width/default_ccf_mask_v_width(NEID2D()); fwtf = 0.5  # using LSF width
#msf = line_width_50_default/default_ccf_mask_v_width(NEID2D()); fwtf = 1.0  # using line width
#msf = 1.24; fwtf = 1.5   # two pixels
 println("# mask_scale_factor = ", msf, ", frac_width_to_fit =  ", fwtf, ". ")
 ((ccfs_espresso, ccf_vars_espresso), v_grid) = ccf_total(order_list_timeseries, line_list_espresso, pipeline_plan,
   mask_type=:gaussian, Δfwhm=Δfwhm,
   mask_scale_factor=msf, range_no_mask_change=5*line_width_50, ccf_mid_velocity=ccf_mid_velocity, v_step=100, #155,
   v_max=max(5*line_width_50,2*max_bc), allow_nans=true, calc_ccf_var=true, recalc=true)
  outputs["v_grid"] = v_grid
  outputs["ccfs_espresso"] = ccfs_espresso
  outputs["ccf_vars_espresso"] = ccf_vars_espresso

ccf_median = vec(median(ccfs_espresso,dims=2))
 ccf_dist_from_median = map(i->RvSpectML.earth_mover_distance(v_grid,ccf_median,ccfs_espresso[:,i]), 1:size(ccfs_espresso,2) )

fwtf = 0.5; println("# mask_scale_factor = ", msf, ", frac_width_to_fit =  ", fwtf, ". ")
 alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=fwtf)
 rvs_ccf_espresso = RvSpectML.calc_rvs_from_ccf_total(ccfs_espresso, ccf_vars_espresso, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)
 σ_rvs_espresso = deepcopy(read_cache(pipeline_plan,:σ_rvs_ccf_total))
 #scatter((order_list_timeseries.times.-first(order_list_timeseries.times)).*24, rvs_ccf_espresso, label=:none)
 outputs["rvs_ccf_espresso"] = rvs_ccf_espresso
 outputs["σ_rvs_espresso"] = σ_rvs_espresso

 delay_until_first_obs_used_for_fit = 1/24
 times = copy(order_list_timeseries.times)
 #Δt_fit = min(time_max_alt_sun-(first(times)+delay_until_first_obs_used_for_fit),last(times)-time_max_alt_sun)
 Δt_fit = 2.0/24
 idx_first_obs_to_fit = findfirst(t->time_max_alt_sun-t<=Δt_fit,times)
 @assert 1 <= idx_first_obs_to_fit length(times)
 idx_last_obs_to_fit = findlast(t->t-time_max_alt_sun<=Δt_fit,times)
 idx_to_fit = idx_first_obs_to_fit:idx_last_obs_to_fit
 idx_to_fit = idx_to_fit[ccf_dist_from_median[idx_to_fit] .< 10]
 (mean_rv_to_fit, std_rv_to_fit) = mean_and_std(rvs_ccf_espresso[idx_to_fit])
 idx_for_robust_fit = idx_to_fit[mean_rv_to_fit-3*std_rv_to_fit .<= rvs_ccf_espresso[idx_to_fit] .<= mean_rv_to_fit+3*std_rv_to_fit]

 times .-= time_max_alt_sun
 #times = order_list_timeseries.times.-first(order_list_timeseries.times)
 x = [times[idx_for_robust_fit]  ones(length(times[idx_for_robust_fit])) ]
 linfit = (x'*x)\x'*(rvs_ccf_espresso[idx_for_robust_fit])
 x = [times  ones(length(times)) ]
 v_fit = (x*linfit)
 v_t0 = [[0] [1]]*linfit
 idx_not_bad_distance = ccf_dist_from_median .< 10
 (times_binned, rvs_binned) = bin_times_and_rvs_max_Δt(times=times[idx_not_bad_distance], rvs=rvs_ccf_espresso[idx_not_bad_distance], Δt_threshold=5/(60*24))
 #idx_to_fit_binned = findfirst(t->t>=times[first(idx_to_fit)], times_binned):findlast(t->t<=times[last(idx_to_fit)], times_binned)
 #idx_to_fit_binned = findfirst(t->time_max_alt_sun-t.<=Δt_fit,times_binned):findlast(t->t-time_max_alt_sun.<=Δt_fit,times_binned)
 idx_to_fit_binned = findfirst(t->t>=-Δt_fit,times_binned):findlast(t->t<=Δt_fit,times_binned)
 x = [times_binned  ones(length(times_binned)) ]
 v_fit_binned = x*linfit
 println("# RMS to linear fit: ", std(rvs_ccf_espresso[idx_to_fit].-v_fit[idx_to_fit]), " m/s.  Binned: ", std(rvs_binned[idx_to_fit_binned].-v_fit_binned[idx_to_fit_binned]))
 println("# Mean RV = ", mean(rvs_binned[idx_to_fit_binned]), " m/s  Slope = ", linfit[1]/24, " m/s/hr")
 outputs["times_binned"] = times_binned
 outputs["rvs_binned"] = rvs_binned
 outputs["rv_rms"] = std(rvs_ccf_espresso[idx_to_fit].-v_fit[idx_to_fit])
 outputs["rv_0"] = v_t0
 outputs["rv_mean"] = mean(rvs_binned[idx_to_fit_binned])
 outputs["rv_slope"] = linfit[1]/24
 outputs["rv_rms_binned"] = std(rvs_binned[idx_to_fit_binned].-v_fit_binned[idx_to_fit_binned])

using Plots
 plt = scatter( (times[idx_not_bad_distance]).*24, (rvs_ccf_espresso.-v_fit)[idx_not_bad_distance],  ms=1.5, label=:none)
 scatter!(plt,(times_binned).*24, (rvs_binned.-v_fit_binned), ms=4.0, label=:none)
 vline!(plt,[-Δt_fit,Δt_fit].* 24, label=:none, color=:grey)
 vline!(plt,[0].* 24, label=:none, color=:orange)
 xlims!(plt,-3.5,3.5)
 #ylims!(plt,minimum(rvs_ccf_espresso.-v_fit)-0.5*std_rv_to_fit, maximum(rvs_ccf_espresso.-v_fit)+0.5*std_rv_to_fit)
 ylims!(plt, extrema((rvs_ccf_espresso.-v_fit)[idx_not_bad_distance]))
 xlabel!(plt,"Time (hr)")
 ylabel!(plt,"RV - linear fit (m/s)")
 title!(plt,"Sun " * date_str * " CCF ESPRESSO\n" *
        "rms_5m = " * string(round(std(rvs_binned[idx_to_fit_binned].-v_fit_binned[idx_to_fit_binned]),digits=3)) * "m/s" *
        " mean = " * string(round(mean(rvs_binned[idx_to_fit_binned]),digits=2)) *
        " Slope = " * string(round(linfit[1]/24,digits=2))  * "m/s/hr")
 display(plt)
 savefig(joinpath(output_dir,"rvs_" * date_str * ".png"))

(ccfs_norm, ccf_vars_norm) = EchelleCCFs.calc_normalized_ccfs( v_grid, ccfs_espresso, ccf_vars_espresso, v_mid=ccf_mid_velocity, dv_min=2.5*line_width_50,dv_max= 2*max_bc )
 ccf_template = EchelleCCFs.calc_ccf_template(ccfs_norm[:,idx_to_fit], ccf_vars_norm[:,idx_to_fit], assume_normalized=true)
 ccf_vars_template = EchelleCCFs.calc_ccf_template(ccf_vars_norm[:,idx_to_fit], ccf_vars_norm[:,idx_to_fit], assume_normalized=true)
 #plot(v_grid,ccfs_espresso, labels=:none)
 #plot(v_grid,ccfs_norm , labels=:none)
 #plt = plot(v_grid,(ccfs_norm.-ccf_template)[:,range(1,size(ccfs_norm,2),step=2)], labels=:none)
 #xlims!(plt,ccf_mid_velocity-2*line_width_50,ccf_mid_velocity+2*line_width_50)
 #heatmap(v_grid,1:size(ccfs_norm,2),(ccfs_norm.-ccf_template)',seriescolor=:balance)
 cols_to_fit = findlast(x->x<ccf_mid_velocity-3*line_width_50,v_grid):findfirst(x->x>ccf_mid_velocity+3*line_width_50,v_grid)
 maxabsz = maximum(abs.(extrema(view(ccfs_norm,cols_to_fit,idx_to_fit)'.-view(ccf_template,cols_to_fit,:)')))
 plt = heatmap(view(v_grid,cols_to_fit),(times[1:size(ccfs_norm,2)]).*24,view(ccfs_norm,cols_to_fit,:)'.-view(ccf_template,cols_to_fit,:)',clims=(-maxabsz,maxabsz),seriescolor=:balance)
 xlabel!(plt,"v (m/s)")
 ylabel!(plt,"Observation Time (hr)")
 title!(plt,"Sun " * date_str * ":\nCCF - <CCF>, ESPRESSO mask + SNR weights")
 #xlims!(ccf_mid_velocity-RvSpectMLBase.max_bc,ccf_mid_velocity+RvSpectMLBase.max_bc)
 xlims!(plt,ccf_mid_velocity-2*line_width_50,ccf_mid_velocity+2*line_width_50)
 hline!(plt,[-Δt_fit,Δt_fit].*24, label=:none, color=:grey)
 hline!(plt,[0.], label=:none, color=:orange)
 ylims!(plt,-3.5,3.0)
 #vline!([-2.1e4], color=4, label=:none)
 display(plt)
 #savefig(joinpath(output_dir,"ccf_heatmaps_sun_" * date_str * "_espresso_msf=" * string(round(msf,digits=1)) * ".png"))
 savefig(joinpath(output_dir,"ccf_heatmaps_sun_" * date_str * "_espresso_msf=" * string(round(msf,digits=1)) * "_new.png"))

 outputs["mean_ccf"] = ccf_template

gp_post = RvSpectML.TemporalGPInterpolation.construct_gp_posterior(v_grid,ccf_template,sigmasq_obs=ccf_vars_template,use_logx=false,smooth_factor=2)
 ccf_template_smooth = RvSpectML.TemporalGPInterpolation.predict_mean(gp_post(v_grid))
 ccf_template_deriv = RvSpectML.numerical_deriv(v_grid,RvSpectML.TemporalGPInterpolation.predict_mean(gp_post(v_grid)))
 norm = sum(abs2.(ccf_template_deriv[cols_to_fit]))
 rvs_proj = vec(sum((ccf_template_smooth[cols_to_fit].-ccfs_norm[cols_to_fit,:]).*ccf_template_deriv[cols_to_fit],dims=1)./norm)
 σ_rvs_proj = vec(sum((ccf_vars_norm[cols_to_fit]).*abs2.(ccf_template_deriv[cols_to_fit]),dims=1))./ norm
 #ccf_resid_minus_rv_proj = (ccfs_norm.-ccf_template).+(rvs_proj.*ccf_template_deriv')'
 ccf_resid_minus_rv_proj = (ccfs_norm.-ccf_template).+((v_fit.-mean(v_fit)).*ccf_template_deriv')'

 outputs["ccf_template_smooth"] = ccf_template_smooth
 outputs["rvs_proj"] = rvs_proj
 outputs["σ_rvs_proj"] = σ_rvs_proj
 outputs["ccf_resid_minus_rv_proj"] = ccf_resid_minus_rv_proj

#plt = heatmap(view(v_grid,cols_to_fit),(times[1:size(ccfs_norm,2)]).*24,ccf_resid_minus_rv_proj[cols_to_fit,:]')
maxabsz = maximum(abs.(extrema(ccf_resid_minus_rv_proj)))
climmax = 0.00030
 plt = heatmap(v_grid,(times[1:size(ccfs_norm,2)]).*24,ccf_resid_minus_rv_proj',clims=(-climmax,climmax),seriescolor=:balance)
 xlabel!(plt,"v (m/s)")
 ylabel!(plt,"Observation Time (hr)")
 title!(plt,"Sun " * date_str * ":\nCCF-<CCF>+v⋅d<CCF>/dv")
 #xlims!(ccf_mid_velocity-RvSpectMLBase.max_bc,ccf_mid_velocity+RvSpectMLBase.max_bc)
 xlims!(plt,ccf_mid_velocity-2*line_width_50,ccf_mid_velocity+2*line_width_50)
 hline!(plt,[-Δt_fit,Δt_fit].*24, label=:none, color=:grey)
 hline!(plt,[0.], label=:none, color=:orange)
 ylims!(plt,-3.5,3.0)
 display(plt)
 #savefig(joinpath(output_dir,"ccf_resid_heatmaps_sun_" * date_str * "_espresso_msf=" * string(round(msf,digits=1)) * ".png"))
 savefig(joinpath(output_dir,"ccf_resid_heatmaps_sun_" * date_str * "_espresso_msf=" * string(round(msf,digits=1)) * "_new.png"))
=#
@assert true

#=
 all_spectra = nothing
 order_list_timeseries = nothing
 GC.gc()
end
=#

#EchelleCCFs.write_each_ccf_fits(map(s->s.metadata,all_spectra), v_grid, ccfs_espresso, ccf_vars_espresso,   fits_hdr=EchelleCCFs.make_ccf_fits_header(all_spectra[1].metadata),
#            line_list_filename=line_list_filename,total_vel=rvs_ccf_espresso, sigma_total_vel=σ_rvs_espresso )

#=
##
max_orders = 56:111
  order_list_timeseries = extract_orders(all_spectra,pipeline_plan, orders_to_use=max_orders, recalc=true )

line_list_filename = "G2.espresso.mas"
 linelist_for_ccf_fn_w_path = joinpath(pkgdir(EchelleCCFs),"data","masks",line_list_filename)
 line_list_espresso = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan, v_center_to_avoid_tellurics=ccf_mid_velocity,
      Δv_to_avoid_tellurics = 1*RvSpectMLBase.max_bc+5*line_width_50+2*default_ccf_mask_v_width(EXPRES1D()),
      orders_to_use=max_orders, recalc=true, verbose=true)

##

good_orders = orders_to_use_default(NEID2D())

(order_ccfs, order_ccf_vars, v_grid_order_ccfs) = ccf_orders(order_list_timeseries, line_list_espresso, pipeline_plan, mask_scale_factor=msf, ccf_mid_velocity=ccf_mid_velocity,
 v_step=155, v_max=max(5*line_width_50,2*max_bc), orders_to_use=good_orders, allow_nans=true, calc_ccf_var=true,recalc=true)
 outputs["order_ccfs"] = order_ccfs
 outputs["order_ccf_vars"] = order_ccf_vars
 outputs["v_grid_order_ccfs"] =v_grid_order_ccfs

obs_id = idx_to_fit[1]
 plt = heatmap(v_grid_order_ccfs, 1:size(order_ccfs,2), (order_ccfs[:,:,obs_id]./mean(order_ccfs[:,:,obs_id],dims=1))')
 xlabel!(plt,"v (m/s)")
 ylabel!(plt,"Order (ignore #)")
 title!(plt,"Sun " * date_str  * " Obs# " * string(obs_id))
 #xlims!(ccf_mid_velocity-RvSpectMLBase.max_bc,ccf_mid_velocity+RvSpectMLBase.max_bc)
 xlims!(plt,ccf_mid_velocity-2.0*line_width_50,ccf_mid_velocity+2.0*line_width_50)
 display(plt)
 savefig(joinpath(output_dir,"order_ccf_heatmaps_sun_" * date_str * "_espresso_msf=" * string(round(msf,digits=1)) * ".png"))

ccf_order_template = zeros(size(order_ccfs,1),size(order_ccfs,2))
 ccf_order_template_vars = zeros(size(order_ccfs,1),size(order_ccfs,2))
 ccf_order_resid = zeros(size(order_ccfs,1),size(order_ccfs,2),size(order_ccfs,3))
 rvs_proj_order = zeros(length(times),size(order_ccfs,2))
 σ_rvs_proj_order = zeros(length(times),size(order_ccfs,2))
 ccf_order_resid_rv_from_order = zeros(size(order_ccfs,1),size(order_ccfs,2),size(order_ccfs,3))

 for ord in 1:size(order_ccfs,2)
     if !(sum(order_ccfs[:,ord,:]) > 0)  continue end
     (ccfs_order_norm, ccf_order_vars_norm) = EchelleCCFs.calc_normalized_ccfs( v_grid, order_ccfs[:,ord,:], order_ccf_vars[:,ord,:], v_mid=ccf_mid_velocity, dv_min=2.5*line_width_50,dv_max= 2*max_bc )
     ccf_order_template[:,ord] .= EchelleCCFs.calc_ccf_template(ccfs_order_norm[:,idx_to_fit], ccf_order_vars_norm[:,idx_to_fit], assume_normalized=true)
     ccf_order_template_vars[:,ord] .= EchelleCCFs.calc_ccf_template(ccf_order_vars_norm[:,idx_to_fit], ccf_order_vars_norm[:,idx_to_fit], assume_normalized=true)
     println("# ord = ", ord, " mean = ", NaNMath.mean(ccf_order_template[:,ord]), " var = ", NaNMath.mean(ccf_order_template_vars[:,ord]))
     gp_post = RvSpectML.TemporalGPInterpolation.construct_gp_posterior(v_grid,ccf_order_template[:,ord],sigmasq_obs=ccf_order_template_vars[:,ord],use_logx=false,smooth_factor=2)
     ccf_template_smooth = RvSpectML.TemporalGPInterpolation.predict_mean(gp_post(v_grid))
     ccf_template_deriv = RvSpectML.numerical_deriv(v_grid,RvSpectML.TemporalGPInterpolation.predict_mean(gp_post(v_grid)))
     ccf_order_resid[:,ord,:] .= (ccfs_order_norm.-ccf_order_template[:,ord]).+(rvs_proj.*ccf_template_deriv')'
     norm = sum(abs2.(ccf_template_deriv[cols_to_fit]))
     rvs_proj_order[:,ord] .= vec(sum((ccf_template_smooth[cols_to_fit].-ccfs_order_norm[cols_to_fit,:]).*ccf_template_deriv[cols_to_fit],dims=1))./norm
     σ_rvs_proj_order[:,ord] .= vec(sum((ccf_order_vars_norm[cols_to_fit,:]).*abs2.(ccf_template_deriv[cols_to_fit]),dims=1))./ norm
     ccf_order_resid_rv_from_order[:,ord,:] .= (ccfs_order_norm.-ccf_template).+(rvs_proj_order[:,ord].*ccf_template_deriv')'
 end
 outputs["ccf_order_template"] = ccf_order_template
 outputs["ccf_order_template_vars"] = ccf_order_template_vars
 outputs["rvs_proj_order"] = rvs_proj_order
 outputs["σ_rvs_proj_order"] = σ_rvs_proj_order
 outputs["ccf_order_resid"] = ccf_order_resid
 outputs["ccf_order_resid_rv_from_order"] = ccf_order_resid_rv_from_order
=#
#=
obs_id = idx_to_fit[6]
 #maxabsz = maximum(abs.(extrema(ccf_order_resid[:,:,obs_id])))
 #plt = heatmap(v_grid_order_ccfs, 1:size(order_ccfs,2), ccf_order_resid[:,:,obs_id]',clims=(-maxabsz,maxabsz),seriescolor=:balance)
 maxabsz = maximum(abs.(extrema(reshape(sum(ccf_order_resid,dims=3),size(ccf_order_resid)[1:2]))))
 plt = heatmap(v_grid_order_ccfs, 1:size(order_ccfs,2), reshape(sum(ccf_order_resid,dims=3),size(ccf_order_resid)[1:2])',clims=(-maxabsz,maxabsz),seriescolor=:balance)
 xlabel!(plt,"v (m/s)")
 ylabel!(plt,"Order (ignore #)")
 title!(plt,"Sun " * date_str  )
 #xlims!(plt,ccf_mid_velocity-RvSpectMLBase.max_bc,ccf_mid_velocity+RvSpectMLBase.max_bc)
 xlims!(plt,ccf_mid_velocity-2.0*line_width_50,ccf_mid_velocity+2.0*line_width_50)
 display(plt)
 #savefig(joinpath(output_dir,"order_ccf_resid_heatmaps_sun_" * date_str * "_espresso_msf=" * string(round(msf,digits=1)) * ".png"))
=#
#=
alg_fit_rv2 =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=fwtf)
order_rvs_g = zeros(length(all_spectra), size(order_ccfs,2))
order_rv_std_g = zeros(length(all_spectra), size(order_ccfs,2))
order_rvs_t = zeros(length(all_spectra), size(order_ccfs,2))
order_rv_std_t = zeros(length(all_spectra), size(order_ccfs,2))
for i in 1:size(order_ccfs,2)
   println("# Index for order: ", i)
   if sum(view(order_ccfs,:,i,1)) > 0
      rvs_order_ccf = RvSpectML.calc_rvs_from_ccf_total(view(order_ccfs,:,i,:), view(order_ccf_vars,:,i,:), pipeline_plan, v_grid=v_grid_order_ccfs, times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv2)
      alg_fit_rv3 =  EchelleCCFs.MeasureRvFromCCFTemplate(v_grid=v_grid_order_ccfs, frac_of_width_to_fit=fwtf, template = vec(mean(view(order_ccfs,:,i,:),dims=2)) )
      order_rvs_g[:,i] .= read_cache(pipeline_plan, :rvs_ccf_total )
      order_rv_std_g[:,i] .= read_cache(pipeline_plan, :σ_rvs_ccf_total)
      rvs_order_ccf = RvSpectML.calc_rvs_from_ccf_total(view(order_ccfs,:,i,:), view(order_ccf_vars,:,i,:), pipeline_plan, v_grid=v_grid_order_ccfs, times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv3)
      order_rvs_t[:,i] .= read_cache(pipeline_plan, :rvs_ccf_total )
      order_rv_std_t[:,i] .= read_cache(pipeline_plan, :σ_rvs_ccf_total)
   end
end
outputs["order_rvs_g"] = order_rvs_g
outputs["order_rvs_t"] = order_rvs_t
=#

#=
outputs[""] =
outputs[""] =
outputs[""] =
outputs[""] =

=#
#=
order_snr = mapreduce(obsid->map(ord->RvSpectMLBase.calc_snr(all_spectra[obsid],787:6214,ord)  ,1:size(first(all_spectra).flux,2)), hcat, 1:length(all_spectra))
outputs["order_snr"] = order_snr

save(outputs_filename, outputs)
=#
#=

EchelleCCFs.write_each_ccf_fits(map(s->s.metadata,all_spectra), #= 100.0 .* EXPRESS USED cm/s =# v_grid, ccfs_espresso, ccf_vars_espresso,
                  fits_hdr=EchelleCCFs.make_ccf_fits_header(all_spectra[1].metadata),
                  line_list_filename=line_list_filename,# orders=161 .-good_orders,ccf_orders=order_ccfs,ccf_orders_var=order_ccf_vars,
                  total_vel=100.0 .* rvs_ccf_espresso, sigma_total_vel= 100.0 .* σ_rvs_espresso,
                  order_vels = 100.0 .* order_rvs_g, sigma_order_vels = 100.0 .* order_rv_std_g )
=#


text:
using Polynomials: Polynomial,fit,coeffs
using Plots
using Distributions: Normal
using Random
using Statistics: mean, std

p = Polynomial([1,2,3])
@show coeffs(p)

errstd = 3
errNormal = Normal(0,errstd)

x = 0:10
y = 2 .+ 3 .* x .+rand(errNormal,11)

scatter(x,y,yerror=errstd,label=false)

fit_line = fit(x,y,1)
@show coeffs(fit_line)
plot!(x,fit_line.(x),label=false)

b_list = []
m_list = []
N = 1000
for i in 1:N
    y = 2 .+ 3 .* x .+rand(errNormal,11)
    fit_c = coeffs(fit(x,y,1))
    push!(b_list,fit_c[1])
    push!(m_list,fit_c[2])
end

@show mean(b_list)
@show mean(m_list)
@show std(b_list)
@show std(m_list)

histogram(b_list,bins=30)
histogram(m_list,bins=30)