text:
export measure_corr

"""
measure correlator.

e.g. measure_corr(chem, 1=>X, 3=>X) will measure <σₓ¹σₓ³> from a quantum MPS.
"""
function measure_corr(chem::QuantumMPS, si::Pair{Int, <:PauliGate}, sj::Pair{Int, <:PauliGate})
    si.first > sj.first && return measure_corr(chem, sj, si)
    si.first == sj.first && throw(ArgumentError("Qubit address conflict error!"))
    T = datatype(chem.initial_reg)

    input_state = chem.input_state
    reg = chem.initial_reg |> copy
    nv = chem.nbit_virtual + chem.nbit_ancilla
    nrep = nrepeat(chem)

    for i = nrep+1:nrep+nv
        input_state[i] == 1 && apply!(reg, put(nv+1, (i-nrep+1)=>XGate{T}()))
    end
    local res = nothing
    input_state[1] == 1 && apply!(reg, put(nv+1, 1=>XGate{T}()))
    for i=1:nrep
        reg |> getblock(chem, i)
        if i!=nrep
            res_i = _measure!(reg, i, [si, sj], input_state, true)
            if res_i != nothing
                if res != nothing
                    return mean(res.*res_i)
                else
                    res = res_i
                end
            end
        end
    end
    for i=1:nv+1-chem.nbit_ancilla
        res_i = _measure!(reg, i+nrep-1, [si, sj], input_state, false)
        if res_i != nothing
            if res != nothing
                return mean(res.*res_i)
            else
                res = res_i
            end
        end
    end
    throw()
end

function _measure!(reg::AbstractRegister{B, T}, i, sl, input_state, reset) where {B, T}
    for si in sl
        if i==si.first
            op_i = eigen!(si.second |> mat |>Matrix)
            reg |> put(nqubits(reg), 1=>matblock(T.(op_i.vectors' |> Matrix)))
            return @inbounds 1 .- 2 .* (reset ? measure_collapseto!(reg, 1; config=input_state[i+1]) : measure_remove!(reg, 1))
        end
    end
    reset ? measure_collapseto!(reg, 1; config=input_state[i+1]) : measure_remove!(reg, 1)
    nothing
end

max_stars_repo_name: GiggleLiu/QuantumMPS

max_issues_repo_name: GiggleLiu/QuantumMPS

max_forks_repo_name: GiggleLiu/QuantumMPS


text:
# Noise Map
using Plots
wnl = 12.0
anl = 0.15
oofnl = 0.2
oofnl2 = 0.0001

white_map = strip2.white_noise(NN, pix_size, wnl)
atmos_map = strip2.atmospheric_noise(NN, pix_size, anl)
one_ove_f = strip2.one_over_f(NN, pix_size, oofnl)
one_ove_f2 = strip2.one_over_f2(NN, pix_size, oofnl2)

noise_map = white_map .+ atmos_map .+ one_ove_f

sky = conv_map .+ noise_map

N_mask = 3.0

filtered_map = strip2.filter_noise(NN, N_mask, sky)

# white noise
strip2.Plot_Sky(NN, pix_size, white_map, cmin=-30, cmax=30)

# instrumental 1/f
strip2.Plot_Sky(NN, pix_size, one_ove_f, cmin=-80, cmax=80)
strip2.Plot_Sky(NN, pix_size, one_ove_f2, cmin=-80, cmax=80)

ell_max = 5000.0
delta_ell = 30.0

a1, b1 = strip2.get_power_spectrum(one_ove_f, one_ove_f,
                                        ell_max, delta_ell, pix_size, NN)

a2, b2 = strip2.get_power_spectrum(one_ove_f2, one_ove_f2,
                                        ell_max, delta_ell, pix_size, NN)

a2, b2 = strip2.get_power_spectrum(Map_E, Map_E,
                                        ell_max, delta_ell, pix_size, NN)

# atmospheric noise
grad = ColorGradient([:blue, :white])
strip2.Plot_Sky(NN, pix_size, strip2.convolve_beam(NN, pix_size, beam_waist, atmos_map), cmin=-100, cmax=100, gradient=grad,
                lab_title="Atmospheric Correlation Map")

# Raw map
strip2.Plot_Sky(NN, pix_size, conv_map .+ atmos_map .+ white_map .+ strip2.convolve_beam(NN, pix_size, beam_waist, one_ove_f), cmin=-500, cmax=500)

# Filtered Map
strip2.Plot_Sky(NN, pix_size, filtered_map .-conv_map, cmin=-400, cmax=400)

# @gif for iter in 1:10
#   atmos_map = strip2.atmospheric_noise(NN, pix_size, anl, iter)
#   grad = ColorGradient([:blue, :white])
#   strip2.Plot_Sky(NN, pix_size, atmos_map, cmin=-100, cmax=100, gradient=grad,
#                   lab_title="Atmospheric Correlation Map")
# end every 1
#

ext: jl


text:
module Issue770

using Gridap
using Test

domain = (0,1,0,1)
cells = (2,2)
model = CartesianDiscreteModel(domain,cells) |> simplexify
Ω = Interior(model)
dΩ = Measure(Ω,4)
m = 2
u(x) = VectorValue(sum(x)^m, sum(x)^m)
reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},m)
V = FESpace(Ω,reffe)
uh = interpolate(u,V)
t = Δ(u) - Δ(uh)
@test sum(∫(t⋅t)dΩ) < 1.0e-12

s(x) = sum(x)^m
reffe = ReferenceFE(lagrangian,Float64,m)
V = FESpace(Ω,reffe)
sh = interpolate(s,V)
t = Δ(s) - Δ(sh)
@test sum(∫(t⋅t)dΩ) < 1.0e-12

end # module


text:
using LaserLab
using Test
# using Unitful
# import Unitful:
#     nm, μm, mm, cm, m, km,
#     mg, g, kg,
#     ps, ns, μs, ms, s, minute, hr, d, yr, Hz, kHz, MHz, GHz,
#     eV,
#     μJ, mJ, J,
# 	μW, mW, W

# function hnu_eV(lambda_nm)
#     return 1240.0 * (1.0/lambda_nm)
# end

# dspot(λ::Real, NA::Real) = 1.83*λ/(2*NA)

#include("dffunctions_test.jl")
#include("setup_test.jl")
#include("glaser_test.jl")
#include("pmt_test.jl")

function pmttest()
	#@testset "pmt_test.jl" begin
        fword  = "/Users/jj/JuliaProjects/LaserLab/test/C1--Trace--00010.trc"
        fbyte  = "/Users/jj/JuliaProjects/LaserLab/test/C1--Trace--00108.trc"

        iob = open(fbyte, "r")
        iow = open(fword, "r")

        @test typeof(iob) == IOStream
        @test typeof(iow) == IOStream

        WAVEDESCb = LaserLab.wavedesc(iob)
        WAVEDESCw = LaserLab.wavedesc(iow)
        @test WAVEDESCb == WAVEDESCw

        @test LaserLab.readword(iob, WAVEDESCb + 32) == 0  #Byte
        @test LaserLab.readword(iow, WAVEDESCw + 32) == 1  #Word

        @test LaserLab.readlong(iow, WAVEDESCw + 60) == 2 * LaserLab.readlong(iob, WAVEDESCb + 60)

        @test LaserLab.readstring(iow, WAVEDESCw + 76, 14) == "LECROYWS4104HD"
        @test LaserLab.readstring(iob, WAVEDESCb + 76, 14) == "LECROYWS4104HD"

        @test LaserLab.readfloat(iow, WAVEDESCw + 160) == LaserLab.readfloat(iob, WAVEDESCb + 160)

        @test LaserLab.readdouble(iow, WAVEDESCw+ 180) ≈ LaserLab.readdouble(iob, WAVEDESCb + 180)

        tsw = LaserLab.readtimestamp(iow, WAVEDESCw + 296)
        @test tsw.year == 2022

        xadd = LaserLab.scope_rawdata(iow, WAVEDESCw)
        @test xadd["NOMINAL_BITS"] == 12

        wvf   = LaserLab.xydata(iow, xadd)
        stats = LaserLab.wstats(wvf; nsigma=3.0)

        @test stats.mean + 3.0*stats.std ≈ stats.thrp
        @test stats.mean - 3.0*stats.std ≈ stats.thrn

    #end
end

max_stars_repo_name: jjgomezcadenas/LaserLab

max_issues_repo_name: jjgomezcadenas/LaserLab

max_forks_repo_name: jjgomezcadenas/LaserLab


text:
using BinDeps
using Compat

# ensure the requisite libraries are available

noNemo = "Nemo.jl is not found:\n  Pkg.rm(\"Nemo\"); Pkg.add(\"Nemo\"); quit()\n  Pkg.rm(\"ArbFloats\");Pkg.add(\"ArbFloats\");"
reNemo = "Nemo.jl is not as expected:\n  Pkg.rm(\"Nemo\"); Pkg.add(\"Nemo\"); quit()\n  Pkg.rm(\"ArbFloats\");Pkg.add(\"ArbFloats\");"

function package_directory(pkgName::String)
    pkgdir = Base.find_in_path(pkgName)
    nothing == pkgdir && throw(ErrorException(noNemo))
    return abspath(joinpath(split(pkgdir, pkgName)[1], pkgName))
end

function library_filepath(libsdir::String, filenames::Vector{String}, libname::String)
    libfile = filenames[ findfirst([startswith(x,libname) for x in filenames]) ]
    return joinpath( libsdir, libfile )
end

NemoLibsDir = abspath(joinpath( package_directory("Nemo"), "local/lib"))
libFiles = readdir(NemoLibsDir)

libarb   = library_filepath( NemoLibsDir, libFiles, "libarb"  )
libflint = library_filepath( NemoLibsDir, libFiles, "libflint")

isfile(libarb)   || throw(ErrorException(reNemo))
isfile(libflint) || throw(ErrorException(reNemo))


text:
type Monomial <: LsmBasisSystemPolynomType end

type MonomialFunction{I <: Integer} <: LSMBasisSystemFunction
  order::I
end

function call(m::MonomialFunction, x::Float64)
  ret = 1.0
  for i = 1:m.order
    ret *= x
  end

  return ret
end

get_type(::Monomial) = MonomialFunction{Int}

function path_basis_system!(::Monomial, order::Int, v::Vector)
  for i = 1:order + 1
    v[i] = MonomialFunction(i - 1)
  end

  return v
end


text:
## discretize!

"""
    discretize!(a::AbstractArray{<:Real}, n::Integer=30)

Sort `a`, divide its values into (roughly) `n` groups corresponding to
equal-length intervals, and replace each value in a group by the `mean`
of this group. This can be useful to transform samples from a continuous
distribution in order to plot an histogram.

# Examples
```julia-repl
julia> pmf(discretize!(randn(1000), 10))
pmf for [-3.2977741401356235, -2.804784333001702, -2.804784333001702, -2.804784333001702, -2.0873488466856744, -2.0873488466856744, -2.0873488466856744, -2.0873488466856744, -2.0873488466856744, -2.0873488466856744  …  2.0792314833569865, 2.0792314833569865, 2.0792314833569865, 2.0792314833569865, 2.0792314833569865, 2.9112746357518313, 2.9112746357518313, 2.9112746357518313, 2.9112746357518313, 2.9112746357518313] with support of length 10:
  -3.29777  => 0.001
  -2.80478  => 0.003
  -2.08735  => 0.026
  -1.36503  => 0.104
  -0.679956 => 0.237
  0.0211377 => 0.273
  0.748744  => 0.232
  1.46951   => 0.101
  2.07923   => 0.018
  2.91127   => 0.005
```
"""
function discretize!(a::AbstractArray{<:Real}, n::Integer=30)
    require_one_based_indexing(a)
    length(a) < 2n && return a
    n = min(n, length(a) ÷ 2)
    sort!(a)
    w = (a[end] - a[1]) / n
    i = 1
    for x in LinRange(a[1], a[end], n)
        j = searchsortedlast(a, x)
        y = mean(view(a, i:j))
        a[i:j] .= y
        i = j+1
    end
    a
end


text:
@timedtestset "dense" begin

    @testset "Array ∘ mortar" begin
        @testset "BlockVector" begin

            # Generate some blocks.
            Ps = [5, 6, 7]
            x = BlockArray(randn(sum(Ps)), Ps).blocks

            # Verify the pullback.
            ȳ = randn(sum(Ps))
            adjoint_test(Array ∘ mortar, ȳ, x)
        end
        @testset "BlockMatrix" begin

            # Generate some blocks.
            Ps = [3, 4, 5]
            Qs = [6, 7, 8, 9]
            X = BlockArray(randn(sum(Ps), sum(Qs)), Ps, Qs).blocks
            Ȳ = randn(sum(Ps), sum(Qs))

            # Verify pullback.
            adjoint_test(Array ∘ mortar, Ȳ, X)
        end
    end
end


text:
#### ApplyBCs

export apply_bcs!

"""
    apply_bcs!

Apply boundary conditions to the state vector.
"""
function apply_bcs! end

function apply_bcs!(grid::Grid, q::StateVec, aux::StateVec, params, case::BOMEX)
    @unpack obukhov_length, SurfaceModel, wstar, UpdVar = params
    ustar = SurfaceModel.ustar
    ρq_tot_flux = SurfaceModel.ρq_tot_flux
    ρθ_liq_flux = SurfaceModel.ρθ_liq_flux
    gm, en, ud, sd, al = allcombinations(q)
    n_updrafts = length(ud)
    k_1 = first_interior(grid, Zmin())
    zLL = grid.zc[k_1]
    θ_liq_1 = q[:θ_liq, k_1, gm]
    q_tot_1 = q[:q_tot, k_1, gm]
    alpha0LL = aux[:α_0, k_1]
    cv_q_tot = surface_variance(
        ρq_tot_flux * alpha0LL,
        ρq_tot_flux * alpha0LL,
        ustar,
        zLL,
        obukhov_length,
    )
    cv_θ_liq = surface_variance(
        ρθ_liq_flux * alpha0LL,
        ρθ_liq_flux * alpha0LL,
        ustar,
        zLL,
        obukhov_length,
    )
    @inbounds for i in ud
        UpdVar[i].surface_bc.a = SurfaceModel.area / n_updrafts
        UpdVar[i].surface_bc.w = 0.0
        UpdVar[i].surface_bc.θ_liq =
            (θ_liq_1 + UpdVar[i].surface_scalar_coeff * sqrt(cv_θ_liq))
        UpdVar[i].surface_bc.q_tot =
            (q_tot_1 + UpdVar[i].surface_scalar_coeff * sqrt(cv_q_tot))
        q[:a, k_1, i] = UpdVar[i].surface_bc.a
        q[:θ_liq, k_1, i] = UpdVar[i].surface_bc.θ_liq
        q[:q_tot, k_1, i] = UpdVar[i].surface_bc.q_tot
    end
    @inbounds for i in sd
        apply_Neumann!(q, :θ_liq, grid, 0.0, (Zmin(), Zmax()), i)
        apply_Neumann!(q, :q_tot, grid, 0.0, (Zmin(), Zmax()), i)
    end
    apply_Dirichlet!(q, :w, grid, 0.0, (Zmin(), Zmax()), en)
    apply_Neumann!(q, :tke, grid, 0.0, (Zmin(), Zmax()), en)
    q[:tke, k_1, en] = surface_tke(ustar, wstar, zLL, obukhov_length)
end


text:
using Statistics
using ProgressMeter
using MATLAB
using MAT
using Dates
include("atn_coupled_module.jl")
@everywhere begin
    using Statistics
    using ProgressMeter
    using MATLAB
    using MAT
    using Dates
    include("atn_coupled_module.jl")
    N = 10_000          # number of nodes
    m = 5               # number of active links per active node
    tmax = 1000         # maximum duration of the simulation
    Ni = Int(0.01 * N)  # initial number of infects
    gamma = -2.1        # exponent of the power law
    lambda_size = 200
    lambda_min = 0.0
    lambda_max = 1.0
    n_trials = 5
    Lambda = range(lambda_min, stop = lambda_max, length = lambda_size)
    beta_size = 200
    beta_min = 0.0
    beta_max = 1.0
    Beta = range(beta_min, stop = beta_max, length = beta_size)
    ActMin = 0.001      # minimum activation probability
    ActMax = 1
    agent_act_1 =
        (
            (ActMax .^ (gamma + 1) - ActMin .^ (gamma + 1)) .* rand(N) .+
            ActMin .^ (gamma + 1)
        ) .^ (1.0 / (gamma + 1))

    agent_act_2 =
        (
            (ActMax .^ (gamma + 1) - ActMin .^ (gamma + 1)) .* rand(N) .+
            ActMin .^ (gamma + 1)
        ) .^ (1.0 / (gamma + 1))
    activity_1 = [agent_act_1]
    activity_2 = [agent_act_2]

    N_iters = beta_size * lambda_size
    for k = 1:(N_iters-1)
        push!(activity_1, agent_act_1)
        push!(activity_2, agent_act_2)
    end
    activity_1 = reshape(activity_1, lambda_size, beta_size)
    activity_2 = reshape(activity_2, lambda_size, beta_size)

    mu_1 = 0.1
    mu_2 = 0.1

    lambda_1 = 1
    # lambda_2 = 0.3

    beta_1 = 0.9
    # beta_2 = beta_1
end

grid = [(beta, lambda) for lambda in Lambda, beta in Beta]

Iend_2 = @showprogress pmap(
    atn_coupled_module.parameter_sweep,
    n_trials * ones(Int, lambda_size, beta_size),
    activity_1,
    activity_2,
    mu_1 * ones(Int, lambda_size, beta_size),
    mu_2 * ones(Int, lambda_size, beta_size),
    lambda_1 * ones(Int, lambda_size, beta_size),
    # beta_1 * ones(Int, lambda_size, beta_size),
    grid,
)

location_mat = "/mnt/storage/MAT/"
right_now = replace(
    replace(replace(string(Dates.now()), "." => "_"), ":" => "_"),
    "-" => "_",
)

file_name = "ATN_SIS_COUPLED" * right_now
mat_name = location_mat * file_name * ".mat"

matwrite(mat_name, Dict("Iend_2" => Iend_2); compress = true)
Iend_2 = Float64.(Iend_2)
mat"figure();set(gcf, 'Position',  [0, 0, 2500, 1500]); hold on;"
mat"axis([ -0,1,-0,10 ])"
mat"imagesc([0,1],[0,10],$Iend_2)"
mat"colormap jet"


text:
include("../src/fdfd_local_field.jl")
using Plots
pyplot()

### reference simulation
refractive_indexes = ones(3) * 1.45
ps = [0.75, 0.5, 0.66, 0.33, 0.66, 0.5, 0.75, 0.5, 0.66, 0.33]
x, y, Ez, dpml, dsource, resolution = simulation_hole_layers_unit_cell(ps, refractive_indexes=refractive_indexes, frequency=0.25)

heatmap(real.(ϵ_hole_layers(x, y, ps, refractive_indexes=refractive_indexes)))
gui()

heatmap(real.(Ez), c=ColorGradient([:red, :white, :blue]))
gui()

### structure simulation
refractive_indexes = Float64[1.0, 1.0, 1.45]
x, y, Ez, dpml, dsource, resolution = simulation_hole_layers_unit_cell(ps, refractive_indexes=refractive_indexes)

heatmap(real.(ϵ_hole_layers(x, y, ps, refractive_indexes=refractive_indexes)))
gui()

heatmap(real.(Ez), c=ColorGradient([:red, :white, :blue]))
gui()


text:
 verbose = true
 make_plots = false
 if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 using RvSpectML
if verbose   println("# Loading other packages")    end
 using CSV, DataFrames, Query, StatsBase, Statistics, Dates
 if verbose   println("# Loading unpackaged code from this repo")    end
 include(joinpath("..","..","src","io.jl"))

 function calc_epoch_offset_vector(rvs_ccf, σ_rvs, tidx_epoch4, tidx_epoch5)
     vcat(rvs_ccf[tidx_epoch4].-mean(rvs_ccf[tidx_epoch4], weights(1.0 ./σ_rvs[tidx_epoch4].^2)), (rvs_ccf[tidx_epoch5].-mean(rvs_ccf[tidx_epoch5],weights(1.0 ./σ_rvs[tidx_epoch5].^2))))
 end

 using Scalpels
 function compute_scalpels_output(mnb::Integer; sort_by_responce::Bool=true)
   (rv_scalpels, scores, ) = calc_clean_rvs_scores_basis_scalpels(rvs_ccf, ccfs, num_basis=mnb, σ_rvs=σ_rvs, sort_by_responce=sort_by_responce)
   rv_scalpels_offsets = (rv_scalpels .- vcat(fill(mean(rv_scalpels[tidx_epoch4], weights(1.0 ./σ_rvs[tidx_epoch4].^2)), sum(tidx_epoch4)), fill(mean(rv_scalpels[tidx_epoch5],weights(1.0 ./σ_rvs[tidx_epoch5].^2)), sum(tidx_epoch5)) ) )

   println("max_num_basis = ", mnb, "  std RV = ", std(rv_scalpels_offsets) )

   output_df = DataFrame("Time [MJD]"=>order_list_timeseries.times, "RV_C [m/s]"=>rv_scalpels_offsets, "e_RV_C [m/s]"=>σ_rvs,
                       "RV_D [m/s]"=>rvs_ccf.-rv_scalpels_offsets, "e_RV_D [m/s]"=>σ_rvs )
   for i in 1:mnb
      output_df[!,"Score_" * string(i) * " [nondim]"] = vec(scores[:,i])
      output_df[!,"e_Score_" * string(i) * " [nondim]"] = fill(NaN,length(rv_scalpels_offsets))
   end
   return output_df
 end

all_spectra = include(joinpath(pkgdir(EchelleInstruments),"examples/read_expres_data_" * string(starid) * ".jl"))
 GC.gc()
 ref_obs_idx = RvSpectMLBase.choose_obs_idx_for_init_guess(df_files_use,get_inst(all_spectra))
 tidx_epoch4 = df_files_use.expres_epoch.==4
 tidx_epoch5 = df_files_use.expres_epoch.==5

# Something buggy if include order 1.  Order 86 essentially ignored due to tellurics. First and last orders have NaN issues
lfc_orders = 43:72
max_orders = 12:83
  order_list_timeseries = extract_orders(all_spectra,pipeline_plan, orders_to_use=lfc_orders, recalc=true )
  #order_list_timeseries = extract_orders(all_spectra,pipeline_plan, orders_to_use=max_orders, recalc=true )
  order_list_timeseries_ref_obs_only = extract_chunklist_timeseries_with_subset_obs(order_list_timeseries, ref_obs_idx)

linelist_for_ccf_fn_w_path = joinpath(pkgdir(EchelleCCFs),"data","masks",espresso_mask_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 = RvSpectMLBase.max_bc, recalc=true, verbose=false)
 ((ccfs, ccf_vars), v_grid) = ccf_total(order_list_timeseries_ref_obs_only, line_list_espresso, pipeline_plan,  mask_scale_factor=2.0, ccf_mid_velocity=ccf_mid_velocity, v_step=100, recalc=true, calc_ccf_var=true)
 line_width_50 = RvSpectMLBase.calc_line_width(v_grid,view(ccfs,:,1),frac_depth=0.5)
 line_sigma = line_width_50/sqrt(2*log(2))
 line_width_05 = RvSpectMLBase.calc_line_width(v_grid,view(ccfs,:,1),frac_depth=0.05)
 σ_rvs_formal_espresso = map(obsid->calc_formal_rv_precission(all_spectra[obsid],RvSpectMLBase.make_chunk_list_from_loc_df(all_spectra[obsid],get_inst(all_spectra), line_list_espresso),smooth_factor=2.0),1:length(all_spectra) )
 alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=2)
 rvs_ccf_ref = RvSpectML.calc_rvs_from_ccf_total( reshape(ccfs,size(ccfs,1),1), reshape(ccf_vars,size(ccf_vars,1),1), pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, bin_nightly=false, recalc=true, alg_fit_rv=alg_fit_rv)
 mean(σ_rvs_formal_espresso), median(σ_rvs_formal_espresso), extrema(σ_rvs_formal_espresso), std(rvs_ccf_ref)

((ccfs, ccf_vars), v_grid) = ccf_total(order_list_timeseries_ref_obs_only, line_list_espresso, pipeline_plan,  mask_scale_factor=4.0, range_no_mask_change=line_width_05, ccf_mid_velocity=ccf_mid_velocity, v_step=100, recalc=true, calc_ccf_var=true)
 alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=2)
 rvs_ccf_ref = RvSpectML.calc_rvs_from_ccf_total( reshape(ccfs,size(ccfs,1),1), reshape(ccf_vars,size(ccf_vars,1),1), pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, bin_nightly=false, recalc=true, alg_fit_rv=alg_fit_rv)
 σ_rv_fit_ref = deepcopy(read_cache(pipeline_plan,:σ_rvs_ccf_total))
 σ_rvs_formal_espresso
 ccf_var_scale_factor = σ_rvs_formal_espresso[ref_obs_idx]/σ_rv_fit_ref[1]

((ccfs_espresso, ccf_vars_espresso), v_grid) = ccf_total(order_list_timeseries, line_list_espresso, pipeline_plan,  mask_scale_factor=4.0, range_no_mask_change=line_width_05, ccf_mid_velocity=ccf_mid_velocity, v_step=100, ccf_var_scale=ccf_var_scale_factor, calc_ccf_var=true, recalc=true)
 # ccf_vars_espresso .*= ccf_var_scale_factor  # If wanted to scale ccf_var manually, rather than passing ccf_var_scale=ccf_var_scale_factor above
 alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=2)
 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=false, alg_fit_rv=alg_fit_rv)
 σ_rvs_ccf_espresso = deepcopy(read_cache(pipeline_plan,:σ_rvs_ccf_total))
 rvs_ccf_espresso, σ_rvs_ccf_espresso
 σ_rvs_ccf_espresso./σ_rvs_formal_espresso
 println("CCF ESPRESSO, LFC: ")
 rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=rvs_ccf_espresso.-mean(rvs_ccf_espresso)))
 rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=rvs_ccf_espresso.-mean(rvs_ccf_espresso))
 println("# RMS of RVs: ", std(rvs_ccf_espresso), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night, " σ_rvs_formal = ", median(σ_rvs_ccf_espresso))

#=
# If you wanted to compute the full covariance matrix for the CCFs.  VERY slow.  And not that different.
((ccfs_espresso2, ccf_covars_espresso2), v_grid2) = ccf_total(order_list_timeseries_ref_obs_only, line_list_espresso, pipeline_plan,  mask_scale_factor=8.0, range_no_mask_change=line_width_05, ccf_mid_velocity=ccf_mid_velocity, calc_ccf_covar=true, recalc=true)
alg_fit_rv2 =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=2)
rvs_ccf_espresso2 = RvSpectML.calc_rvs_from_ccf_total(ccfs_espresso2, ccf_covars_espresso2, pipeline_plan, v_grid=v_grid2, times = order_list_timeseries.times, recalc=true, bin_nightly=false, alg_fit_rv=alg_fit_rv2)
rvs_ccf_espresso2, σ_rvs_ccf_espresso2
=#

num_basis_scalpels = 2
  println("# ESPRESSO mask, LFC, Scalpels, num_basis_vectors = ", num_basis_scalpels)
  ccfs = ccfs_espresso; ccf_vars = ccf_vars_espresso; rvs_ccf = rvs_ccf_espresso; σ_rvs = σ_rvs_ccf_espresso; println("# Setting ccfs to ESPRESSO mask ", num_basis_scalpels, " basis vectors");  @time df_scalpels0 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
  CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), df_scalpels0)
  rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels0[!,Symbol("RV_C [m/s]")].-mean(df_scalpels0[!,Symbol("RV_C [m/s]")])))
  rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels0[!,Symbol("RV_C [m/s]")].-mean(df_scalpels0[!,Symbol("RV_C [m/s]")]))
  println("# RMS of RVs: ", std(df_scalpels0[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

#= =#

chunk_list_timeseries = RvSpectMLBase.make_chunk_list_timeseries_telluric_free(all_spectra, line_list_espresso, min_pixels_in_chunk=5, verbose=false)
 dont_need_to!(pipeline_plan,:extract_orders)
 rv_offset_epoch4_minus_epoch5 = calc_epoch_offset_vector(rvs_ccf_espresso, σ_rvs_formal_espresso, tidx_epoch4, tidx_epoch5 )
 set_rv_est!(chunk_list_timeseries, rvs_ccf_espresso .+ tidx_epoch5.*rv_offset_epoch4_minus_epoch5 )
 chunk_grids = RvSpectML.make_grids_for_chunklist_timeseries(chunk_list_timeseries)
 λmatrix = reduce(vcat,chunk_grids)
 idx_use_pca = trues(length(all_spectra));
 if starid == 101501 idx_use_pca[2] = false  end
 size(λmatrix)

 #using Plots
 #plt = scatter(chunk_list_timeseries.times,rvs_ccf_espresso.-mean(rvs_ccf_espresso),yerr=σ_rvs_ccf_espresso, label=:none, xlabel="Time (d)", ylabel="RV (m/s)", title="RVs CCF, ESPRESSO, LFC")
 #plt = scatter(rvs_ccf_espresso.-mean(rvs_ccf_espresso), label="RVs CCF, ESPRESSO LFC")
 #savefig("rvs_" * string(starid) * "_espresso_lfc.png")

 num_basis_dcpca = 2
 println("# ESPRESSO mask, LFC, DCPCA, num_basis_vectors = ", num_basis_dcpca)
 #for num_basis_dcpca in 1:6
   dcpca_result_espresso = RvSpectML.DCPCA.clean_rvs_dcpca(chunk_list_timeseries, rvs_ccf_espresso, num_basis_vectors_pca=num_basis_dcpca, idx_use_pca=idx_use_pca, chunk_grids=chunk_grids, num_samples_σ_scores=40)
   #print("DCPCA ", num_basis_dcpca, ": " ) #std(rvs) = ", std(dcpca_result_espresso.rvs_clean))
   rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=dcpca_result_espresso.rvs_clean.-mean(dcpca_result_espresso.rvs_clean)))
   rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=dcpca_result_espresso.rvs_clean.-mean(dcpca_result_espresso.rvs_clean))
   println("# RMS of RVs: ", std(dcpca_result_espresso.rvs_clean), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

   #scatter!(plt,dcpca_result_espresso.rvs_dirty, label=string(num_basis_dcpca))
 #end
 #display(plt)
 σ_rvs_clean_espresso = σ_rvs_formal_espresso # sqrt.(σ_rvs_formal_espresso.^2 .+ dcpca_result_espresso.σ_rvs_dirty.^2 )
 df_espresso_dcpca0 = make_results_df(order_list_timeseries.times,dcpca_result_espresso.rvs_clean,σ_rvs_clean_espresso,dcpca_result_espresso.rvs_dirty,NaN.+dcpca_result_espresso.σ_rvs_dirty,dcpca_result_espresso.scores,NaN.+dcpca_result_espresso.σ_scores)
 CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_dcpca0_results.csv"), df_espresso_dcpca0)

#dcpca_result_espresso = RvSpectML.DCPCA.clean_rvs_dcpca(chunk_list_timeseries, rvs_ccf_espresso, num_basis_vectors_pca=6, idx_use_pca=idx_use_pca, chunk_grids=chunk_grids, num_samples_σ_scores=10)
#dcpca_result_espresso.σ_scores

#=
using Plots
scatter(chunk_list_timeseries.times)
plt = scatter(dcpca_result_espresso.rvs_dirty)
plt = scatter(dcpca_result_espresso.scores[:,1])
scatter!(plt,dcpca_result_espresso.scores[:,1])
scatter!(plt,dcpca_result_espresso.scores[:,2])
scatter!(plt,dcpca_result_espresso.scores[:,3])
scatter!(plt,dcpca_result_espresso.scores[:,4])
scatter!(plt,dcpca_result_espresso.scores[:,5])
=#

#scatter(chunk_list_timeseries.times,dcpca_result_espresso.rvs_clean.-mean(dcpca_result_espresso.rvs_clean))
#scatter!(chunk_list_timeseries.times,dcpca_result_espresso.rvs_dirty)

#=
plt = scatter(λmatrix,dcpca_result_espresso.model.mean,ms=1)
xlims!(5200,5400.0)
scatter!(plt,λmatrix,dcpca_result_espresso.model.mean,ms=2)
scatter!(plt,λmatrix,1.0.+dcpca_result_espresso.model.proj[:,1],ms=2)
plt =scatter(λmatrix,1.0.+dcpca_result_espresso.model.proj[:,1],ms=1)
xlims!(5200,5210.0)
scatter!(plt,λmatrix,1.0.+dcpca_result_espresso.model.proj[:,2],ms=1.5)
plot!(plt,λmatrix,1.0.+dcpca_result_espresso.model.proj[:,3],ms=1.5)
xlims!(5202.5,5204.16)
xlims!(5298,5302)
xlims!(5590,5593)
xlims!(5299,5300)
xlims!(5340,5400)
xlims!(5347,5348)
xlims!(5393,5395)
xlims!(6159,6162)
xlims!(6158,6165)
xlims!(5211.0,5213.)
xlims!(5445.0,5455.)
xlims!(5540.0,5555.)
xlims!(5546.0,5550.)
xlims!(5250.0,5450.)
xlims!(5440.0,5450.)
xlims!(5300.0,5310.)
xlims!(5430.0,5450.)

idxstart = sum(length.(chunk_grids)[1:94]); idxstop = idxstart + 30; scatter!(plt,λmatrix[idxstart:idxstop],1.0.+dcpca_result_espresso.model.proj[idxstart:idxstop,3],ms=2)

plt = scatter(λmatrix,1.0.+dcpca_result_espresso.model.proj[:,2],ms=1.5)
xlims!(5400,5580)
plt = plot(λmatrix,dcpca_result_espresso.model.mean,ms=1)
xlims!(5500,5580)
scatter!(plt,λmatrix,1.0.+10.0 .* dcpca_result_espresso.model.proj[:,2],ms=1.5)
ylims!(0.9,1.1)
xlims!(5500,5510)
=#

println("CCF custom mask common, LFC: ")
linelist_for_ccf_fn_w_path = joinpath("..","..","data",string(starid),psu_mask3_filename)
 line_list_clean1 = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan, v_center_to_avoid_tellurics=0*ccf_mid_velocity, Δv_to_avoid_tellurics = RvSpectMLBase.max_bc, recalc=true, verbose=false)
 σ_rvs_formal_clean1 = map(obsid->calc_formal_rv_precission(all_spectra[obsid],RvSpectMLBase.make_chunk_list_from_loc_df(all_spectra[obsid],get_inst(all_spectra), line_list_clean1),smooth_factor=2.0),1:length(all_spectra) )
 mean(σ_rvs_formal_clean1), median(σ_rvs_formal_clean1), extrema(σ_rvs_formal_clean1)
t_start_1 = now()
 ((ccfs_clean1, ccf_vars_clean1), v_grid) = ccf_total(order_list_timeseries, line_list_clean1, pipeline_plan,  mask_scale_factor=4.0, range_no_mask_change=line_width_05, ccf_mid_velocity=0*ccf_mid_velocity, v_step=100, calc_ccf_var=true, recalc=true)

 alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=2.0)
 rvs_ccf_clean1 = calc_rvs_from_ccf_total(ccfs_clean1, ccf_vars_clean1, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)
 σ_rvs_clean1 = deepcopy(read_cache(pipeline_plan,:σ_rvs_ccf_total))
 #std(rvs_ccf_clean1)
 rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=rvs_ccf_clean1.-mean(rvs_ccf_clean1)))
 rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=rvs_ccf_clean1.-mean(rvs_ccf_clean1))
 println("# RMS of RVs: ", std(rvs_ccf_clean1), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night, " σ_rvs_formal = ", median(σ_rvs_formal_clean1))

 num_basis_scalpels = 2
   println("# CCF custom mask common, LFC, Scalpels, num_basis_vectors = ", num_basis_scalpels)
   ccfs = ccfs_clean1; ccf_vars = ccf_vars_clean1; rvs_ccf = rvs_ccf_clean1; σ_rvs = σ_rvs_clean1; println("# Setting ccfs to clean1 mask ", num_basis_scalpels, " basis vectors");
   @time df_scalpels1 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
   CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels3_results.csv"), df_scalpels1)
   #plt = scatter(df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")]), yerr=df_scalpels1[!,Symbol("e_RV_C [m/s]")], label="1")
   rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")])))
   rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")]))
   println("# RMS of RVs: ", std(df_scalpels1[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

println("# CCF custom mask common, LFC, DCPCA, num_basis_vectors = ", num_basis_dcpca)
chunk_list_timeseries = RvSpectMLBase.make_chunk_list_timeseries_telluric_free(all_spectra, line_list_clean1, min_pixels_in_chunk=5, verbose=false)
 dont_need_to!(pipeline_plan,:extract_orders)
 rv_offset_epoch4_minus_epoch5 = calc_epoch_offset_vector(rvs_ccf_clean1, σ_rvs_formal_clean1, tidx_epoch4, tidx_epoch5 )
 set_rv_est!(chunk_list_timeseries, rvs_ccf_clean1 .+ tidx_epoch5.*rv_offset_epoch4_minus_epoch5 )
 chunk_grids = RvSpectML.make_grids_for_chunklist_timeseries(chunk_list_timeseries)
 idx_use_pca = trues(length(all_spectra)); if starid == 101501 idx_use_pca[2] = false  end
 dcpca_result_clean1 = RvSpectML.DCPCA.clean_rvs_dcpca(chunk_list_timeseries, rvs_ccf_clean1, num_basis_vectors_pca=2, idx_use_pca=idx_use_pca, chunk_grids=chunk_grids, num_samples_σ_scores=40)
 σ_rvs_clean_clean1 = σ_rvs_formal_clean1 # sqrt.(σ_rvs_formal_clean1.^2 .+ dcpca_result_clean1.σ_rvs_dirty.^2 )
 df_clean1 = make_results_df(order_list_timeseries.times,dcpca_result_clean1.rvs_clean,σ_rvs_clean_clean1,dcpca_result_clean1.rvs_dirty,NaN.+dcpca_result_clean1.σ_rvs_dirty,dcpca_result_clean1.scores,NaN.+dcpca_result_clean1.σ_scores)
 t_stop_1 = now()
 CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_dcpca3_results.csv"), df_clean1)
 t_stop_1-t_start_1
 std(rvs_ccf_clean1), std(dcpca_result_clean1.rvs_clean)
 rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_clean1[!,Symbol("RV_C [m/s]")].-mean(df_clean1[!,Symbol("RV_C [m/s]")])))
 rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_clean1[!,Symbol("RV_C [m/s]")].-mean(df_clean1[!,Symbol("RV_C [m/s]")]))
 println("# RMS of RVs: ", std(df_clean1[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

#=
println("CCF custom mask common: std(rvs) = ", std(rvs_ccf_clean1))
 plt = scatter(rvs_ccf_clean1.-mean(rvs_ccf_clean1), label="RV CCF")
  #for num_basis_dcpca in 1:6
  println("# Custom mask common, LFC, DCPCA, num_basis_vectors = ", num_basis_dcpca)
    dcpca_result_clean1 = RvSpectML.DCPCA.clean_rvs_dcpca(chunk_list_timeseries, rvs_ccf_clean1, num_basis_vectors_pca=num_basis_dcpca, idx_use_pca=idx_use_pca, chunk_grids=chunk_grids, num_samples_σ_scores=0)
    #print("DCPCA ", num_basis_dcpca, ": " ) #std(rvs) = ", std(dcpca_result_espresso.rvs_clean))
    rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=dcpca_result_clean1.rvs_clean.-mean(dcpca_result_clean1.rvs_clean)))
    rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=dcpca_result_clean1.rvs_clean.-mean(dcpca_result_clean1.rvs_clean))
    println("# RMS of RVs: ", std(dcpca_result_clean1.rvs_clean), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

    scatter!(plt,dcpca_result_clean1.rvs_dirty, label=string(num_basis_dcpca))
  #end
  σ_rvs_clean_espresso = σ_rvs_formal_espresso # sqrt.(σ_rvs_formal_espresso.^2 .+ dcpca_result_espresso.σ_rvs_dirty.^2 )
  df_clean1_dcpca1 = make_results_df(order_list_timeseries.times,dcpca_result_espresso.rvs_clean,σ_rvs_clean_espresso,dcpca_result_espresso.rvs_dirty,NaN.+dcpca_result_espresso.σ_rvs_dirty,dcpca_result_espresso.scores,NaN.+dcpca_result_espresso.σ_scores)
  CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_dcpca3_results.csv"), df_clean1_dcpca1)
  display(plt)
=#

#=
using Plots
scatter(order_list_timeseries.times,rvs_ccf_espresso.-mean(rvs_ccf_espresso),yerr=σ_rvs_ccf_espresso)
scatter!(order_list_timeseries.times,rvs_ccf_clean1.-mean(rvs_ccf_clean1),yerr=σ_rvs_clean_clean1)
std(rvs_ccf_clean1)
scatter(order_list_timeseries.times,rvs_ccf_clean1.-mean(rvs_ccf_clean1) .- (rvs_ccf_espresso.-mean(rvs_ccf_espresso)),yerr=σ_rvs_ccf_espresso, xlabel="Time (d)", ylabel="ΔRV (m/s)")
=#

#=
num_basis_scalpels = 1
 println("# Custom mask common, LFC, Scalpels, num_basis_vectors = ", num_basis_scalpels)
 ccfs = ccfs_clean1; ccf_vars = ccf_vars_clean1; rvs_ccf = rvs_ccf_clean1; σ_rvs = σ_rvs_clean_clean1; println("# Setting ccfs to clean1 mask ", num_basis_scalpels, " basis vectors");
 @time df_scalpels1 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
 #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), df_scalpels1)
 plt = scatter(df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")]), yerr=df_scalpels1[!,Symbol("e_RV_C [m/s]")], label="1")
 rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")])))
 rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")]))
 println("# RMS of RVs: ", std(df_scalpels1[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

num_basis_scalpels = 2
  ccfs = ccfs_clean1; ccf_vars = ccf_vars_clean1; rvs_ccf = rvs_ccf_clean1; σ_rvs = σ_rvs_clean_clean1; println("# Setting ccfs to clean1 mask ", num_basis_scalpels, " basis vectors");
  @time df_scalpels1 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
  CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels3_results.csv"), df_scalpels1)
  plt = scatter(df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")]), yerr=df_scalpels1[!,Symbol("e_RV_C [m/s]")], label="1")
  rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")])))
  rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")]))
  println("# RMS of RVs: ", std(df_scalpels1[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

 num_basis_scalpels = 3
   ccfs = ccfs_clean1; ccf_vars = ccf_vars_clean1; rvs_ccf = rvs_ccf_clean1; σ_rvs = σ_rvs_clean_clean1; println("# Setting ccfs to clean1 mask ", num_basis_scalpels, " basis vectors");
   @time df_scalpels1 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
   #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), df_scalpels1)
   plt = scatter(df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")]), yerr=df_scalpels1[!,Symbol("e_RV_C [m/s]")], label="1")
   rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")])))
   rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")]))
   println("# RMS of RVs: ", std(df_scalpels1[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

 num_basis_scalpels = 4
    ccfs = ccfs_clean1; ccf_vars = ccf_vars_clean1; rvs_ccf = rvs_ccf_clean1; σ_rvs = σ_rvs_clean_clean1; println("# Setting ccfs to clean1 mask ", num_basis_scalpels, " basis vectors");
    @time df_scalpels1 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
    #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), df_scalpels1)
    plt = scatter(df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")]), yerr=df_scalpels1[!,Symbol("e_RV_C [m/s]")], label="1")
    rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")])))
    rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")]))
    println("# RMS of RVs: ", std(df_scalpels1[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)
 num_basis_scalpels = 5

     ccfs = ccfs_clean1; ccf_vars = ccf_vars_clean1; rvs_ccf = rvs_ccf_clean1; σ_rvs = σ_rvs_clean_clean1; println("# Setting ccfs to clean1 mask ", num_basis_scalpels, " basis vectors");
     @time df_scalpels1 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
     #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), df_scalpels1)
     plt = scatter(df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")]), yerr=df_scalpels1[!,Symbol("e_RV_C [m/s]")], label="1")
     rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")])))
     rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")]))
     println("# RMS of RVs: ", std(df_scalpels1[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

 num_basis_scalpels = 6
      ccfs = ccfs_clean1; ccf_vars = ccf_vars_clean1; rvs_ccf = rvs_ccf_clean1; σ_rvs = σ_rvs_clean_clean1; println("# Setting ccfs to clean1 mask ", num_basis_scalpels, " basis vectors");
      @time df_scalpels1 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
      #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), df_scalpels1)
      plt = scatter(df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")]), yerr=df_scalpels1[!,Symbol("e_RV_C [m/s]")], label="1")
      rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")])))
      rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels1[!,Symbol("RV_C [m/s]")].-mean(df_scalpels1[!,Symbol("RV_C [m/s]")]))
      println("# RMS of RVs: ", std(df_scalpels1[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)
=#

println("CCF custom mask best, LFC: ")
linelist_for_ccf_fn_w_path = joinpath("..","..","data",string(starid),psu_mask4_filename)
  line_list_clean2 = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan, v_center_to_avoid_tellurics=0*ccf_mid_velocity, Δv_to_avoid_tellurics = RvSpectMLBase.max_bc, recalc=true, verbose=false)
  σ_rvs_formal_clean2 = map(obsid->calc_formal_rv_precission(all_spectra[obsid],RvSpectMLBase.make_chunk_list_from_loc_df(all_spectra[obsid],get_inst(all_spectra), line_list_clean2),smooth_factor=2.0),1:length(all_spectra) )
  mean(σ_rvs_formal_clean2), median(σ_rvs_formal_clean2), extrema(σ_rvs_formal_clean2)
  t_start_2 = now()
  ((ccfs_clean2, ccf_vars_clean2), v_grid) = ccf_total(order_list_timeseries, line_list_clean2, pipeline_plan,  mask_scale_factor=4.0, range_no_mask_change=line_width_05, ccf_mid_velocity=0*ccf_mid_velocity, v_step=100, calc_ccf_var=true, recalc=true)
  alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=2)
  rvs_ccf_clean2 = calc_rvs_from_ccf_total(ccfs_clean2, ccf_vars_clean2, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)
  σ_rvs_clean2 = deepcopy(read_cache(pipeline_plan,:σ_rvs_ccf_total))
  rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=rvs_ccf_clean2.-mean(rvs_ccf_clean2)))
  rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=rvs_ccf_clean2.-mean(rvs_ccf_clean2))
  println("# RMS of RVs: ", std(rvs_ccf_clean2), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night, " σ_rvs_formal = ", median(σ_rvs_formal_clean2))

  num_basis_scalpels = 2
    println("# CCF custom mask best, LFC, Scalpels, num_basis_vectors = ", num_basis_scalpels)
    ccfs = ccfs_clean2; ccf_vars = ccf_vars_clean2; rvs_ccf = rvs_ccf_clean2; σ_rvs = σ_rvs_clean2; #println("# Setting ccfs to clean2 mask ", num_basis_scalpels, " basis vectors");
    @time df_scalpels2 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
    CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels4_results.csv"), df_scalpels1)
    #plt = scatter(df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")]), yerr=df_scalpels2[!,Symbol("e_RV_C [m/s]")], label="1")
    rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")])))
    rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")]))
    println("# RMS of RVs: ", std(df_scalpels2[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

println("# Custom mask best, DCPCA: ")
  chunk_list_timeseries = RvSpectMLBase.make_chunk_list_timeseries_telluric_free(all_spectra, line_list_clean2, min_pixels_in_chunk=5, verbose=false)
  dont_need_to!(pipeline_plan,:extract_orders)
  rv_offset_epoch4_minus_epoch5 = calc_epoch_offset_vector(rvs_ccf_clean2, σ_rvs_formal_clean2, tidx_epoch4, tidx_epoch5 )
  set_rv_est!(chunk_list_timeseries, rvs_ccf_clean2 .+ tidx_epoch5.*rv_offset_epoch4_minus_epoch5 )
  chunk_grids = RvSpectML.make_grids_for_chunklist_timeseries(chunk_list_timeseries)
  idx_use_pca = trues(length(all_spectra)); if starid == 101501 idx_use_pca[2] = false  end
  dcpca_result_clean2 = RvSpectML.DCPCA.clean_rvs_dcpca(chunk_list_timeseries, rvs_ccf_clean2, num_basis_vectors_pca=2, idx_use_pca=idx_use_pca, chunk_grids=chunk_grids, num_samples_σ_scores=40)
  σ_rvs_clean_clean2 = σ_rvs_formal_clean2 # sqrt.(σ_rvs_formal_clean2.^2 .+ dcpca_result_clean2.σ_rvs_dirty.^2 )
  df_clean2 = make_results_df(order_list_timeseries.times,dcpca_result_clean2.rvs_clean,σ_rvs_clean_clean2,dcpca_result_clean2.rvs_dirty,NaN.+dcpca_result_clean2.σ_rvs_dirty,dcpca_result_clean2.scores,NaN.+dcpca_result_clean2.σ_scores)
  t_stop_2 = now()
  rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_clean2[!,Symbol("RV_C [m/s]")].-mean(df_clean2[!,Symbol("RV_C [m/s]")])))
  rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_clean2[!,Symbol("RV_C [m/s]")].-mean(df_clean2[!,Symbol("RV_C [m/s]")]))
  println("# RMS of RVs: ", std(df_clean2[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)
  CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_dcpca4_results.csv"), df_clean2)

  t_stop_2-t_start_2
#std(dcpca_result_espresso.rvs_clean)
#std(dcpca_result_clean1.rvs_clean), std(dcpca_result_clean2.rvs_clean)
#df_clean2

#=
num_basis_scalpels = 1
 ccfs = ccfs_clean2; ccf_vars = ccf_vars_clean2; rvs_ccf = rvs_ccf_clean2; σ_rvs = σ_rvs_clean_clean2; println("# Setting ccfs to clean2 mask ", num_basis_scalpels, " basis vectors");
 @time df_scalpels2 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
 #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels4_results.csv"), df_scalpels2)
 plt = scatter(df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")]), yerr=df_scalpels2[!,Symbol("e_RV_C [m/s]")], label="1")
 rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")])))
 rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")]))
 println("# RMS of RVs: ", std(df_scalpels2[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

 num_basis_scalpels = 2
  ccfs = ccfs_clean2; ccf_vars = ccf_vars_clean2; rvs_ccf = rvs_ccf_clean2; σ_rvs = σ_rvs_clean_clean2; println("# Setting ccfs to clean2 mask ", num_basis_scalpels, " basis vectors");
  @time df_scalpels2 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
  #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels4_results.csv"), df_scalpels2)
  plt = scatter(df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")]), yerr=df_scalpels2[!,Symbol("e_RV_C [m/s]")], label="1")
  rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")])))
  rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")]))
  println("# RMS of RVs: ", std(df_scalpels2[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

 num_basis_scalpels = 3
   ccfs = ccfs_clean2; ccf_vars = ccf_vars_clean2; rvs_ccf = rvs_ccf_clean2; σ_rvs = σ_rvs_clean_clean2; println("# Setting ccfs to clean2 mask ", num_basis_scalpels, " basis vectors");
   @time df_scalpels2 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
   #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels4_results.csv"), df_scalpels2)
   plt = scatter(df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")]), yerr=df_scalpels2[!,Symbol("e_RV_C [m/s]")], label="1")
   rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")])))
   rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")]))
   println("# RMS of RVs: ", std(df_scalpels2[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

 num_basis_scalpels = 4
    ccfs = ccfs_clean2; ccf_vars = ccf_vars_clean2; rvs_ccf = rvs_ccf_clean2; σ_rvs = σ_rvs_clean_clean2; println("# Setting ccfs to clean2 mask ", num_basis_scalpels, " basis vectors");
    @time df_scalpels2 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
    #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels4_results.csv"), df_scalpels2)
    plt = scatter(df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")]), yerr=df_scalpels2[!,Symbol("e_RV_C [m/s]")], label="1")
    rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")])))
    rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")]))
    println("# RMS of RVs: ", std(df_scalpels2[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

 num_basis_scalpels = 5
     ccfs = ccfs_clean2; ccf_vars = ccf_vars_clean2; rvs_ccf = rvs_ccf_clean2; σ_rvs = σ_rvs_clean_clean2; println("# Setting ccfs to clean2 mask ", num_basis_scalpels, " basis vectors");
     @time df_scalpels2 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
     #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels4_results.csv"), df_scalpels2)
     plt = scatter(df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")]), yerr=df_scalpels2[!,Symbol("e_RV_C [m/s]")], label="1")
     rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")])))
     rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")]))
     println("# RMS of RVs: ", std(df_scalpels2[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)

 num_basis_scalpels = 6
      ccfs = ccfs_clean2; ccf_vars = ccf_vars_clean2; rvs_ccf = rvs_ccf_clean2; σ_rvs = σ_rvs_clean_clean2; println("# Setting ccfs to clean2 mask ", num_basis_scalpels, " basis vectors");
      @time df_scalpels2 = compute_scalpels_output(num_basis_scalpels, sort_by_responce=true)
      #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels4_results.csv"), df_scalpels2)
      plt = scatter(df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")]), yerr=df_scalpels2[!,Symbol("e_RV_C [m/s]")], label="1")
      rms_rv_nightly = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")])))
      rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=df_scalpels2[!,Symbol("RV_C [m/s]")].-mean(df_scalpels2[!,Symbol("RV_C [m/s]")]))
      println("# RMS of RVs: ", std(df_scalpels2[!,Symbol("RV_C [m/s]")]), "  nightly RVs: ", rms_rv_nightly, "  within night: ",rms_rv_within_night)
=#

#=
println("std rvs = ", std(rvs_ccf_espresso) )
num_basis_scalpels = 1
 ccfs = ccfs_espresso; ccf_vars = ccf_vars_espresso; rvs_ccf = rvs_ccf_espresso; σ_rvs = σ_rvs_formal_espresso; println("# Setting ccfs to ESPRESSO mask ", num_basis_scalpels, " basis vectors");
 @time output_df1 = compute_scalpels_output(num_basis_scalpels,sort_by_responce=true)
 #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), output_df)
 plt = scatter(output_df1[!,Symbol("RV_C [m/s]")].-mean(output_df1[!,Symbol("RV_C [m/s]")]), yerr=output_df1[!,Symbol("e_RV_C [m/s]")], label="1")

num_basis_scalpels = 2
 @time output_df = compute_scalpels_output(num_basis_scalpels)
 #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), output_df)
 plt = scatter!(output_df[!,Symbol("RV_C [m/s]")].-mean(output_df[!,Symbol("RV_C [m/s]")])#= , yerr=output_df[!,Symbol("e_RV_C [m/s]")] =#, label="2")

num_basis_scalpels = 3
  @time output_df3 = compute_scalpels_output(num_basis_scalpels)
  #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), output_df)
  plt = scatter(output_df3[!,Symbol("RV_C [m/s]")].-mean(output_df3[!,Symbol("RV_C [m/s]")]),  #=yerr=output_df3[!,Symbol("e_RV_C [m/s]")],=# label="3")

num_basis_scalpels = 4
  @time output_df4 = compute_scalpels_output(num_basis_scalpels)
  #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), output_df)
  plt = scatter!(output_df4[!,Symbol("RV_C [m/s]")].-mean(output_df4[!,Symbol("RV_C [m/s]")]),  #=yerr=output_df4[!,Symbol("e_RV_C [m/s]")], =# label="4")

num_basis_scalpels = 5
  @time output_df5 = compute_scalpels_output(num_basis_scalpels)
  #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), output_df)
  plt = scatter!(output_df5[!,Symbol("RV_C [m/s]")].-mean(output_df5[!,Symbol("RV_C [m/s]")]), #=yerr=output_df5[!,Symbol("e_RV_C [m/s]")], =# label="5")

num_basis_scalpels = 6
   @time output_df6 = compute_scalpels_output(num_basis_scalpels)
   #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), output_df)
   plt = scatter!(output_df6[!,Symbol("RV_C [m/s]")].-mean(output_df6[!,Symbol("RV_C [m/s]")]), #=yerr=output_df6[!,Symbol("e_RV_C [m/s]")], =# label="6")

num_basis_scalpels = 7
     @time output_df6 = compute_scalpels_output(num_basis_scalpels)
     #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), output_df)
     plt = scatter!(output_df6[!,Symbol("RV_C [m/s]")].-mean(output_df6[!,Symbol("RV_C [m/s]")]), #=yerr=output_df6[!,Symbol("e_RV_C [m/s]")], =# label="6")

num_basis_scalpels = 8
          @time output_df6 = compute_scalpels_output(num_basis_scalpels)
          #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), output_df)
          plt = scatter!(output_df6[!,Symbol("RV_C [m/s]")].-mean(output_df6[!,Symbol("RV_C [m/s]")]), #=yerr=output_df6[!,Symbol("e_RV_C [m/s]")], =# label="6")

num_basis_scalpels = 9
                    @time output_df6 = compute_scalpels_output(num_basis_scalpels)
                    #CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels0_results.csv"), output_df)
                    plt = scatter!(output_df6[!,Symbol("RV_C [m/s]")].-mean(output_df6[!,Symbol("RV_C [m/s]")]), #=yerr=output_df6[!,Symbol("e_RV_C [m/s]")], =# label="6")

scatter!(dcpca_result_espresso.rvs_clean.-mean(dcpca_result_espresso.rvs_clean), yerr=σ_rvs_clean_espresso)
# scatter!(df_espresso[!,Symbol("RV_C [m/s]")], yerr=df_espresso[!,Symbol("e_RV_C [m/s]")])

(scalpels_score,scalpels_basis) = Scalpels.calc_basis_scores_scalpels(rvs_ccf, ccfs, σ_rvs = σ_rvs, num_basis=12, assume_centered=false, sort_by_responce=false )

scatter(scalpels_score[:,1])
scatter!(scalpes_score[:,2])
scatter!(scalpels_score[:,3])
scatter(scalpels_score[:,4])
scatter!(scalpels_score[:,5])
scatter(scalpels_score[:,6])

scatter(chunk_list_timeseries.times,scalpels_score[:,1])
scatter!(chunk_list_timeseries.times,scalpels_score[:,2])
scatter!(chunk_list_timeseries.times,scalpels_score[:,3])
scatter(chunk_list_timeseries.times,scalpels_score[:,4])
scatter!(chunk_list_timeseries.times,scalpels_score[:,5])
scatter(chunk_list_timeseries.times,scalpels_score[:,6])

plot(scalpels_basis[:,1])
plot!(scalpels_basis[:,2])
plot!(scalpels_basis[:,3])
plot!(scalpels_basis[:,4])
plot!(scalpels_basis[:,5])
plot!(scalpels_basis[:,6])
scatter!(scalpels_basis[:,7])
scatter!(scalpels_basis[:,8])
scatter!(scalpels_basis[:,9])
scatter!(scalpels_basis[:,10])
scatter!(scalpels_basis[:,11])
scatter!(scalpels_basis[:,12])

scatter(scalpels_score[:,4])

output_df[tidx_epoch4,Symbol("RV_C [m/s]")]
=#

 #=
 ccfs = ccfs_clean1; ccf_vars = ccf_vars_clean1; rvs_ccf = rvs_ccf_clean1; σ_rvs = σ_rvs_formal_clean1;  println("# Setting ccfs to clean1 mask");
 output_df = compute_scalpels_output(num_basis_scalpels)
 CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels3_results.csv"), output_df)

 ccfs = ccfs_clean2; ccf_vars = ccf_vars_clean2; rvs_ccf = rvs_ccf_clean2; σ_rvs = σ_rvs_formal_clean2;  println("# Setting ccfs to clean2 mask");
 output_df = compute_scalpels_output(num_basis_scalpels)
 CSV.write(joinpath("..","..","pennstate",string(starid),string(starid) * "_pennstate_scalpels4_results.csv"), output_df)
 =#

println("# Done")
#= =#


text:
abstract type CoordinateDifferentiableFunction end

"""
  Set internal parameters of the function f at the point x.
"""
initialize!(f, x) = error("initialize! not implemented for $(typeof(f))")

"""
  Coordinate k of the gradient of f evaluated at x.
"""
gradient(f, x, k) = error("gradient not implemented for $(typeof(f))")

"""
  This should return number of coordinates or blocks of coordinates over
  which the coordinate descent iterates.
"""
numCoordinates(f) = error("numCoordinates not implemented for $(typeof(f))")

"""
  Arguments:

  * f is CoordinateDifferentiableFunction
  * g is a prox function

  This function does two things:

  * It finds h such that f(x + e_k⋅h) + g(x_k + h) decreses f(x) + g(x).
    Often h = arg_min f(x + e_k⋅h) + g(x_k + h), but it could also
    minimize a local quadratic approximation.

  * The function also updates its internals. This is done by expecting
    that the algorithm to call this function again is coordinate descent.
    In future, we may want to implement other variants of coordinate descent.
"""
descendCoordinate!(f, g, x, k) = error("descendCoordinate not implemented for $(typeof(f))")

####################################
#
# loss |Y - X⋅β|^2 / (2⋅n)
#
####################################
struct CDLeastSquaresLoss{T<:AbstractFloat, S, U} <: CoordinateDifferentiableFunction
  y::S
  X::U
  r::Vector{T}

  CDLeastSquaresLoss{T, S, U}(y::AbstractVector{T}, X::AbstractMatrix{T}, r::Vector{T}) where {T,S,U} =
    new(y,X,r)
end

function CDLeastSquaresLoss(y::AbstractVector{T}, X::AbstractMatrix{T}) where {T<:AbstractFloat}
  length(y) == size(X, 1) || throw(DimensionMismatch())
  CDLeastSquaresLoss{T, typeof(y), typeof(X)}(y,X,copy(y))
end

numCoordinates(f::CDLeastSquaresLoss) = size(f.X, 2)

function initialize!(f::CDLeastSquaresLoss{T}, x::SparseIterate{T}) where {T<:AbstractFloat}
  # compute residuals for the loss

  X = f.X
  y = f.y
  r = f.r

  n = size(X, 1)

  @simd for i=1:n
    @inbounds r[i] = y[i] - A_mul_B_row(X, x, i)
  end
  nothing
end

gradient(f::CDLeastSquaresLoss{T}, x::SparseIterate{T}, j::Int64) where {T<:AbstractFloat} =
  - At_mul_B_row(f.X, f.r, j) / size(f.X, 1)

# a = X[:, k]' X[:, k]
# b = X[:, k]' r
#
# h        = arg_min a/(2n) (h-b/a)^2 + λ_k⋅|x_k + h|
# xnew[k]  = arg_min a/(2n) (xnew_k - (x_k + b/a))^2 + λ_k⋅|xnew_k|
function descendCoordinate!(
  f::CDLeastSquaresLoss{T},
  g::ProxL1{T},
  x::SparseIterate{T},
  k::Int64) where {T<:AbstractFloat}

  y = f.y
  X = f.X
  r = f.r
  n = length(f.y)

  a = zero(T)
  b = zero(T)
  @inbounds @simd for i=1:n
    a += X[i, k] * X[i, k]
    b += r[i] * X[i, k]
  end

  oldVal = x[k]
  x[k] += b / a
  newVal = cdprox!(g, x, k, n / a)
  h = newVal - oldVal

  # update internals -- residuls = y - X * xnew
  @inbounds @simd for i=1:n
    r[i] -= X[i, k] * h
  end
  h
end

################################################################################
#
# loss ∑_i w_i ⋅ |Y_i - X_i⋅β|^2 / (2⋅n)
#
################################################################################
struct CDWeightedLSLoss{T<:AbstractFloat, S, U} <: CoordinateDifferentiableFunction
  y::S
  X::U
  w::S
  r::Vector{T}

  CDWeightedLSLoss{T, S, U}(y::AbstractVector{T}, X::AbstractMatrix{T}, w::AbstractVector{T}, r::Vector{T}) where {T,S,U} =
    new(y,X,w,r)
end

function CDWeightedLSLoss(y::AbstractVector{T}, X::AbstractMatrix{T}, w::AbstractVector{T}) where {T<:AbstractFloat}
  length(y) == size(X, 1) == length(w) || throw(DimensionMismatch())
  CDWeightedLSLoss{T, typeof(y), typeof(X)}(y,X,w,copy(y))
end

numCoordinates(f::CDWeightedLSLoss) = size(f.X, 2)

function initialize!(f::CDWeightedLSLoss{T}, x::SparseIterate{T}) where {T<:AbstractFloat}
  # compute residuals for the loss

  X = f.X
  y = f.y
  r = f.r

  n = size(X, 1)

  @simd for i=1:n
    @inbounds r[i] = y[i] - A_mul_B_row(X, x, i)
  end
  nothing
end

function gradient(f::CDWeightedLSLoss{T}, x::SparseIterate{T}, j::Int64) where {T<:AbstractFloat}
  out = zero(T)

  n = length(f.r)
  @inbounds @simd for i=1:n
    out += f.w[i] * f.X[i, j] * f.r[i]
  end
  - out / n
end

# a = X[:, k]' X[:, k]
# b = X[:, k]' r
#
# h        = arg_min a/(2n) (h-b/a)^2 + λ_k⋅|x_k + h|
# xnew[k]  = arg_min a/(2n) (xnew_k - (x_k + b/a))^2 + λ_k⋅|xnew_k|
function descendCoordinate!(
  f::CDWeightedLSLoss{T},
  g::ProxL1{T},
  x::SparseIterate{T},
  k::Int64) where {T<:AbstractFloat}

  y = f.y
  X = f.X
  r = f.r
  w = f.w
  n = length(f.y)

  a = zero(T)
  b = zero(T)
  @inbounds @simd for i=1:n
    a += X[i, k] * X[i, k] * w[i]
    b += r[i] * X[i, k] * w[i]
  end

  oldVal = x[k]
  x[k] += b / a
  newVal = cdprox!(g, x, k, n / a)
  h = newVal - oldVal

  # update internals -- residuls = y - X * xnew
  @inbounds @simd for i=1:n
    r[i] -= X[i, k] * h
  end
  h
end

####################################
#
# loss |Y - X⋅β|_2 / sqrt(n)
#
####################################
struct CDSqrtLassoLoss{T<:AbstractFloat, S, U} <: CoordinateDifferentiableFunction
  y::S
  X::U
  r::Vector{T}

  CDSqrtLassoLoss{T, S, U}(y::AbstractVector{T}, X::AbstractMatrix{T}, r::Vector{T}) where {T,S,U} =
    new(y,X,r)
end

function CDSqrtLassoLoss(y::AbstractVector{T}, X::AbstractMatrix{T}) where {T<:AbstractFloat}
  length(y) == size(X, 1) || throw(DimensionMismatch())
  CDSqrtLassoLoss{T, typeof(y), typeof(X)}(y,X,copy(y))
end

numCoordinates(f::CDSqrtLassoLoss) = size(f.X, 2)

function initialize!(f::CDSqrtLassoLoss{T}, x::SparseIterate{T}) where {T<:AbstractFloat}
  # compute residuals for the loss

  X = f.X
  y = f.y
  r = f.r

  n, p = size(f.X)

  @simd for i=1:n
    @inbounds r[i] = y[i] - A_mul_B_row(X, x, i)
  end
  nothing
end

gradient(f::CDSqrtLassoLoss{T}, x::SparseIterate{T}, j::Int64) where {T<:AbstractFloat} =
  - At_mul_B_row(f.X, f.r, j) / norm(f.r)

# a = X[:, k]' X[:, k]
# b = X[:, k]' r
#
# h        = arg_min a/(2n) (h-b/a)^2 + λ_k⋅|x_k + h|
# xnew[k]  = arg_min a/(2n) (xnew_k - (x_k + b/a))^2 + λ_k⋅|xnew_k|
function descendCoordinate!(
  f::CDSqrtLassoLoss{T},
  g::ProxL1{T},
  x::SparseIterate{T},
  k::Int64) where {T<:AbstractFloat}

  y = f.y
  X = f.X
  r = f.r
  n = length(f.y)

  # residuls = y - X * x + X[:, k] * x[k]
  @inbounds @simd for i=1:n
    r[i] += X[i, k] * x[k]
  end
  # r = y - X*x + X[:,k]*x[k]

  s = zero(T)
  xsqr = zero(T)
  rsqr = zero(T)
  @inbounds @simd for i=1:n
    xsqr += X[i, k] * X[i, k]
    s += r[i] * X[i, k]
    rsqr += r[i] * r[i]
  end
  # s = dot(r, X[:,k])
  # xsqr = dot(X[:,k], X[:, k])
  # rsqr = dot(r, r)

  λ = g.λ0
  if !isa(g, ProxL1{T, Nothing})
    λ *= g.λ[k]
  end

  oldVal = x[k]
  if abs(s) <= λ * sqrt(rsqr)
    x[k] = zero(T)
  elseif s > λ * sqrt(rsqr)
    x[k] = ( s - λ / sqrt(1 - λ^2 / xsqr) * sqrt(rsqr - s^2/xsqr) ) / xsqr
  else
    x[k] = ( s + λ / sqrt(1 - λ^2 / xsqr) * sqrt(rsqr - s^2/xsqr) ) / xsqr
  end

  # update internals -- residuls = y - X * xnew
  @inbounds @simd for i=1:n
    r[i] -= X[i, k] * x[k]
  end

  x[k] - oldVal
end

####################################
#
# quadratic x'Ax/2 + x'b
#
####################################
struct CDQuadraticLoss{T<:AbstractFloat, S, U} <: CoordinateDifferentiableFunction
  A::S
  b::U
  Ax::Vector{T}
end

function CDQuadraticLoss(A::AbstractMatrix{T}, b::AbstractVector{T}) where {T<:AbstractFloat}
  (issymmetric(A) && length(b) == size(A, 2)) || throw(ArgumentError())
  CDQuadraticLoss{T, typeof(A), typeof(b)}(A, b, zeros(length(b)))
end

numCoordinates(f::CDQuadraticLoss) = length(f.b)
function initialize!(f::CDQuadraticLoss, x::SparseIterate)
  Ax = f.Ax
  A = f.A
  p = length(Ax)

  @simd for i=1:p
    @inbounds Ax[i] = A_mul_B_row(A, x, i)
  end
  nothing
end
gradient(f::CDQuadraticLoss{T}, ::SparseIterate{T}, j::Int64) where {T<:AbstractFloat} =
  f.Ax[j] + f.b[j]

function descendCoordinate!(
  f::CDQuadraticLoss{T},
  g::ProxL1{T},
  x::SparseIterate{T},
  k::Int64) where {T<:AbstractFloat}

  a = f.A[k,k]
  b = gradient(f, x, k)

  oldVal = x[k]
  a = one(T) / a
  x[k] -= b * a
  newVal = cdprox!(g, x, k, a)
  h = newVal - oldVal

  # update Ax
  Ax = f.Ax
  A  = f.A
  p = length(f.Ax)
  @simd for i=1:p
    @inbounds Ax[i] += A[i,k] * h
  end

  h
end


text:
module Revise

using FileWatching, REPL, Base.CoreLogging, Distributed
import LibGit2

using OrderedCollections: OrderedDict

export revise, includet

"""
    Revise.watching_files[]

Returns `true` if we watch files rather than their containing directory.
FreeBSD and NFS-mounted systems should watch files, otherwise we prefer to watch
directories.
"""
const watching_files = Ref(Sys.KERNEL == :FreeBSD)

"""
    Revise.polling_files[]

Returns `true` if we should poll the filesystem for changes to the files that define
loaded code. It is preferable to avoid polling, instead relying on operating system
notifications via `FileWatching.watch_file`. However, NFS-mounted
filesystems (and perhaps others) do not support file-watching, so for code stored
on such filesystems you should turn polling on.

See the documentation for the `JULIA_REVISE_POLL` environment variable.
"""
const polling_files = Ref(false)
function wait_changed(file)
    polling_files[] ? poll_file(file) : watch_file(file)
    return nothing
end

"""
    Revise.tracking_Main_includes[]

Returns `true` if files directly included from the REPL should be tracked.
The default is `false`. See the documentation regarding the `JULIA_REVISE_INCLUDE`
environment variable to customize it.
"""
const tracking_Main_includes = Ref(false)

include("relocatable_exprs.jl")
include("types.jl")
include("parsing.jl")
include("exprutils.jl")
include("pkgs.jl")
include("git.jl")
include("recipes.jl")

### Globals to keep track of state

"""
    Revise.watched_files

Global variable, `watched_files[dirname]` returns the collection of files in `dirname`
that we're monitoring for changes. The returned value has type [`WatchList`](@ref).

This variable allows us to watch directories rather than files, reducing the burden on
the OS.
"""
const watched_files = Dict{String,WatchList}()

"""
    Revise.revision_queue

Global variable, `revision_queue` holds the names of files that we need to revise, meaning
that these files have changed since we last processed a revision.
This list gets populated by callbacks that watch directories for updates.
"""
const revision_queue = Set{String}()

"""
    Revise.fileinfos

`fileinfos` is the core information that tracks the relationship between source code
and julia objects, and allows re-evaluation of code in the proper module scope.
It is a dictionary indexed by absolute paths of files;
`fileinfos[filename]` returns a value of type [`Revise.FileInfo`](@ref).
"""
const fileinfos = Dict{String,FileInfo}()

"""
    Revise.module2files

`module2files` holds the list of filenames used to define a particular
module. This is only used by `revise(MyModule)` to "refresh" all the definitions in a module.
"""
const module2files = Dict{Symbol,Vector{String}}()

"""
    Revise.included_files

Global variable, `included_files` gets populated by callbacks we register with `include`.
It's used to track non-precompiled packages and, optionally, user scripts (see docs on
`JULIA_REVISE_INCLUDE`).
"""
const included_files = Tuple{Module,String}[]  # (module, filename)

"""
    Revise.basesrccache

Full path to the running Julia's cache of source code defining `Base`.
"""
const basesrccache = joinpath(Sys.BINDIR, Base.DATAROOTDIR, "julia", "base.cache")

"""
    Revise.juliadir

Constant specifying full path to julia top-level directory from which julia was built.
This is reliable even for cross-builds.
"""
const juliadir = begin
    basefiles = map(x->x[2], Base._included_files)
    sysimg = filter(x->endswith(x, "sysimg.jl"), basefiles)[1]
    dirname(dirname(sysimg))
end

"""
    Revise.dont_watch_pkgs

Global variable, use `push!(Revise.dont_watch_pkgs, :MyPackage)` to prevent Revise
from tracking changes to `MyPackage`. You can do this from the REPL or from your
`.julia/config/startup.jl` file.

See also [`Revise.silence`](@ref).
"""
const dont_watch_pkgs = Set{Symbol}()
const silence_pkgs = Set{Symbol}()
const depsdir = joinpath(dirname(@__DIR__), "deps")
const silencefile = Ref(joinpath(depsdir, "silence.txt"))  # Ref so that tests don't clobber

function use_compiled_modules()
    return Base.JLOptions().use_compiled_modules != 0
end

"""
    fmrep = eval_revised(fmnew::FileModules, fmref::FileModules)

Implement the changes from `fmref` to `fmnew`, returning a replacement [`FileModules`](@ref)
`fmrep`.
"""
function eval_revised(fmnew::FileModules, fmref::FileModules)
    fmrep = FileModules(first(keys(fmref)))  # replacement for fmref
    for (mod, fmmnew) in fmnew
        fmrep[mod] = fmm = FMMaps()
        if haskey(fmref, mod)
            eval_revised!(fmm, mod, fmmnew, fmref[mod])
            for p in workers()
                p == myid() && continue
                remotecall(eval_revised_dummy!, p, mod, fmmnew, fmref[mod])
            end
        else  # a new submodule (see #43)
            eval_and_insert_all!(fmm, mod, fmmnew.defmap)
            for p in workers()
                p == myid() && continue
                remotecall(eval_and_insert_all_dummy!, p, mod, fmmnew.defmap)
            end
        end
    end
    fmrep
end

function eval_revised!(fmmrep::FMMaps, mod::Module,
                       fmmnew::FMMaps, fmmref::FMMaps)
    # Update to the state of fmmnew, preventing any unnecessary evaluation
    for (def,val) in fmmnew.defmap
        @assert def != nothing
        defref = getkey(fmmref.defmap, def, nothing)
        if defref != nothing
            # The same expression is found in both, only update the lineoffset
            if val !== nothing
                sigtref = fmmref.defmap[defref][1]
                lnref = firstlineno(defref)
                lnnew = firstlineno(def)
                lineoffset = (isa(lnref, Integer) && isa(lnnew, Integer)) ? lnnew-lnref : 0
                fmmrep.defmap[defref] = (sigtref, lineoffset)
                for sigt in sigtref
                    fmmrep.sigtmap[sigt] = defref
                end
            else
                fmmrep.defmap[defref] = nothing
            end
        else
            eval_and_insert!(fmmrep, mod, def=>val)
        end
    end
    # Delete any methods missing in fmmnew
    for (sigt,_) in fmmref.sigtmap
        if !haskey(fmmrep.sigtmap, sigt)
            m = get_method(sigt)
            if isa(m, Method)
                Base.delete_method(m)
            else
                mths = Base._methods_by_ftype(sigt, -1, typemax(UInt))
                io = IOBuffer()
                println(io, "Extracted method table:")
                println(io, mths)
                info = String(take!(io))
                @warn "Revise failed to find any methods for signature $sigt\n  Perhaps it was already deleted.\n$info"
            end
        end
    end
    return fmmrep
end

eval_revised_dummy!(mod::Module, fmmnew::FMMaps, fmmref::FMMaps) =
    eval_revised!(FMMaps(), mod, fmmnew, fmmref)

function eval_and_insert!(fmm::FMMaps, mod::Module, pr::Pair)
    def, val = pr.first, pr.second
    ex = convert(Expr, def)
    try
        if isdocexpr(ex) && mod == Base.__toplevel__
            Core.eval(Main, ex)
        else
            Core.eval(mod, ex)
        end
        if val isa RelocatableExpr
            instantiate_sigs!(fmm, def, val, mod)
        else
            fmm.defmap[def] = val
            if val !== nothing
                for sigt in val[1]
                    fmm.sigtmap[sigt] = def
                end
            end
        end
    catch err
        @error "failure to evaluate changes in $mod"
        showerror(stderr, err)
        println_maxsize(stderr, "\n", ex; maxlines=20)
    end
    return fmm
end

function eval_and_insert_all!(fmm::FMMaps, mod::Module, defmap::DefMap)
    for pr in defmap
        eval_and_insert!(fmm, mod, pr)
    end
end

eval_and_insert_all_dummy!(mod::Module, defmap::DefMap) =
    eval_and_insert_all!(FMMaps(), mod, defmap)

"""
    Revise.init_watching(files)

For every filename in `files`, monitor the filesystem for updates. When the file is
updated, either [`revise_dir_queued`](@ref) or [`revise_file_queued`](@ref) will
be called.
"""
function init_watching(files)
    udirs = Set{String}()
    for file in files
        dir, basename = splitdir(file)
        haskey(watched_files, dir) || (watched_files[dir] = WatchList())
        push!(watched_files[dir], basename)
        if watching_files[]
            @async revise_file_queued(file)
        else
            push!(udirs, dir)
        end
    end
    for dir in udirs
        updatetime!(watched_files[dir])
        @async revise_dir_queued(dir)
    end
    return nothing
end

"""
    revise_dir_queued(dirname)

Wait for one or more of the files registered in `Revise.watched_files[dirname]` to be
modified, and then queue the corresponding files on [`Revise.revision_queue`](@ref).
This is generally called within an `@async`.
"""
function revise_dir_queued(dirname)
    if !isdir(dirname)
        sleep(0.1)   # in case git has done a delete/replace cycle
        if !isfile(dirname)
            with_logger(SimpleLogger(stderr)) do
                @warn "$dirname is not an existing directory, Revise is not watching"
            end
            return nothing
        end
    end
    latestfiles = watch_files_via_dir(dirname)  # will block here until file(s) change
    for file in latestfiles
        push!(revision_queue, file)
    end
    @async revise_dir_queued(dirname)
end

# See #66.
"""
    revise_file_queued(filename)

Wait for modifications to `filename`, and then queue the corresponding files on [`Revise.revision_queue`](@ref).
This is generally called within an `@async`.

This is used only on platforms (like BSD) which cannot use [`revise_dir_queued`](@ref).
"""
function revise_file_queued(file)
    if !isfile(file)
        sleep(0.1)  # in case git has done a delete/replace cycle
        if !isfile(file)
            with_logger(SimpleLogger(stderr)) do
                @error "$file is not an existing file, Revise is not watching"
            end
            return nothing
        end
    end

    wait_changed(file)  # will block here until the file changes
    push!(revision_queue, file)
    @async revise_file_queued(file)
end

"""
    Revise.revise_file_now(file)

Process revisions to `file`. This parses `file` and computes an expression-level diff
between the current state of the file and its most recently evaluated state.
It then deletes any removed methods and re-evaluates any changed expressions.

`file` must be a key in [`Revise.fileinfos`](@ref)
"""
function revise_file_now(file)
    if !haskey(fileinfos, file)
        println("Revise is currently tracking the following files: ", keys(fileinfos))
        error(file, " is not currently being tracked.")
    end
    fi = fileinfos[file]
    maybe_parse_from_cache!(fi, file)
    fmref = fi.fm
    topmod = first(keys(fi.fm))
    fmnew = parse_source(file, topmod)
    if fmnew != nothing
        fmrep = eval_revised(fmnew, fmref)
        fileinfos[file] = FileInfo(fmrep, fi)
    end
    nothing
end

function instantiate_sigs!(fm::FileModules)
    for (mod, fmm) in fm
        instantiate_sigs!(fmm, mod)
    end
    return fm
end

function instantiate_sigs!(fmm::FMMaps, mod::Module)
    for (def, sig) in fmm.defmap
        if sig isa RelocatableExpr
            instantiate_sigs!(fmm, def, sig, mod)
        end
    end
    return fmm
end

function instantiate_sigs!(fmm::FMMaps, def::RelocatableExpr, sig::RelocatableExpr, mod::Module)
    # Generate the signature-types
    sigtexs = sig_type_exprs(sig)
    local sigts
    try
        sigts = Any[Core.eval(mod, s) for s in sigtexs]
    catch err
        sigwarn(mod, sigtexs, def)
        rethrow(err)
    end
    # Insert into the maps
    fmm.defmap[def] = (sigts, 0)
    for sigt in sigts
        fmm.sigtmap[sigt] = def
    end
    return def
end
@noinline sigwarn(mod, sigtexs, def) = @warn "error processing module $mod signature expressions $sigtexs from $def"

"""
    revise()

`eval` any changes in the revision queue. See [`Revise.revision_queue`](@ref).
"""
function revise()
    sleep(0.01)  # in case the file system isn't quite done writing out the new files
    for file in revision_queue
        revise_file_now(file)
    end
    empty!(revision_queue)
    tracking_Main_includes[] && queue_includes(Main)
    nothing
end

# This variant avoids unhandled task failures
function revise(backend::REPL.REPLBackend)
    sleep(0.01)  # in case the file system isn't quite done writing out the new files
    for file in revision_queue
        try
            revise_file_now(file)
        catch err
            put!(backend.response_channel, (err, catch_backtrace()))
        end
    end
    empty!(revision_queue)
    tracking_Main_includes[] && queue_includes(Main)
    nothing
end

"""
    revise(mod::Module)

Reevaluate every definition in `mod`, whether it was changed or not. This is useful
to propagate an updated macro definition, or to force recompiling generated functions.
"""
function revise(mod::Module)
    for file in module2files[Symbol(mod)]
        for (mod,fmm) in fileinfos[file].fm
            for def in keys(fmm.defmap)
                Core.eval(mod, convert(Expr, def))
            end
        end
    end
    return true  # fixme try/catch?
end

"""
    Revise.track(mod::Module, file::AbstractString)
    Revise.track(file::AbstractString)

Watch `file` for updates and [`revise`](@ref) loaded code with any
changes. `mod` is the module into which `file` is evaluated; if omitted,
it defaults to `Main`.
"""
function track(mod::Module, file::AbstractString)
    isfile(file) || error(file, " is not a file")
    file = normpath(abspath(file))
    fm = parse_source(file, mod)
    if fm != nothing
        instantiate_sigs!(fm)
        fileinfos[file] = FileInfo(fm)
    end
    init_watching((file,))
end

track(file::AbstractString) = track(Main, file)

"""
    includet(filename)

Include and track `filename`.
"""
includet(file::AbstractString) = (Base.include(Main, file); track(Main, file))

"""
    Revise.silence(pkg)

Silence warnings about not tracking changes to package `pkg`.
"""
function silence(pkg::Symbol)
    push!(silence_pkgs, pkg)
    if !isdir(depsdir)
        mkpath(depsdir)
    end
    open(silencefile[], "w") do io
        for p in silence_pkgs
            println(io, p)
        end
    end
    nothing
end
silence(pkg::AbstractString) = silence(Symbol(pkg))

## Utilities

"""
    method = get_method(sigt)

Get the method `method` with signature-type `sigt`. This is used to provide
the method to `Base.delete_method`. See also [`get_signature`](@ref).

If `sigt` does not correspond to a method, returns `nothing`.

# Examples

```jldoctest; setup = :(using Revise), filter = r"in Main at.*"
julia> mymethod(::Int) = 1
mymethod (generic function with 1 method)

julia> mymethod(::AbstractFloat) = 2
mymethod (generic function with 2 methods)

julia> Revise.get_method(Tuple{typeof(mymethod), Int})
mymethod(::Int64) in Main at REPL[0]:1

julia> Revise.get_method(Tuple{typeof(mymethod), Float64})
mymethod(::AbstractFloat) in Main at REPL[1]:1

julia> Revise.get_method(Tuple{typeof(mymethod), Number})

```
"""
function get_method(@nospecialize(sigt))
    mths = Base._methods_by_ftype(sigt, -1, typemax(UInt))
    length(mths) == 1 && return mths[1][3]
    if !isempty(mths)
        # There might be many methods, but the one that should match should be the
        # last one, since methods are ordered by specificity
        i = lastindex(mths)
        while i > 0
            m = mths[i][3]
            m.sig == sigt && return m
            i -= 1
        end
    end
    return nothing
end

"""
    rex = get_def(method::Method)

Return the RelocatableExpr defining `method`.
The source-file defining `method` must be tracked.
If it is in Base, this will execute `track(Base)` if necessary.
"""
function get_def(method::Method)
    filename = String(method.file)
    startswith(filename, "REPL") && error("methods defined at the REPL are not yet supported")
    if !haskey(fileinfos, filename)
        # See whether it's in Base
        basefile = Base.find_source_file(filename)
        basefile == nothing && error("file $filename must be tracked by Revise")
        filename = realpath(basefile)
    end
    if !haskey(fileinfos, filename)
        @info "tracking Base"
        track(Base)
    end
    fi = fileinfos[filename]
    maybe_parse_from_cache!(fi, filename)
    map = fi.fm[method.module].sigtmap
    haskey(map, method.sig) && return map[method.sig]
    @warn "$(method.sig) was not found, perhaps it was generated by code"
    nothing
end

function printf_maxsize(f::Function, io::IO, args...; maxchars::Integer=500, maxlines::Integer=20)
    # This is dumb but certain to work
    iotmp = IOBuffer()
    for a in args
        print(iotmp, a)
    end
    print(iotmp, '\n')
    seek(iotmp, 0)
    str = read(iotmp, String)
    if length(str) > maxchars
        str = first(str, (maxchars+1)÷2) * "…" * last(str, maxchars - (maxchars+1)÷2)
    end
    lines = split(str, '\n')
    if length(lines) <= maxlines
        for line in lines
            f(io, line)
        end
        return
    end
    half = (maxlines+1) ÷ 2
    for i = 1:half
        f(io, lines[i])
    end
    maxlines > 1 && f(io, ⋮)
    for i = length(lines) - (maxlines-half) + 1:length(lines)
        f(io, lines[i])
    end
end
println_maxsize(args...; kwargs...) = println_maxsize(stdout, args...; kwargs...)
println_maxsize(io::IO, args...; kwargs...) = printf_maxsize(println, stdout, args...; kwargs...)

function fix_line_statements!(ex::Expr, file::Symbol, line_offset::Int=0)
    if ex.head == :line
        ex.args[1] += line_offset
        ex.args[2] = file
    else
        for (i, a) in enumerate(ex.args)
            if isa(a, Expr)
                fix_line_statements!(a::Expr, file, line_offset)
            elseif isa(a, LineNumberNode)
                ex.args[i] = file_line_statement(a::LineNumberNode, file, line_offset)
            end
        end
    end
    ex
end

file_line_statement(lnn::LineNumberNode, file::Symbol, line_offset) =
    LineNumberNode(lnn.line + line_offset, file)

function update_stacktrace_lineno!(trace)
    for i = 1:length(trace)
        t, n = trace[i]
        if t.linfo isa Core.MethodInstance
            sigt = t.linfo.def.sig
            file = String(t.file)
            if haskey(fileinfos, file)
                fm = fileinfos[file].fm
                for (mod, fmm) in fm
                    if haskey(fmm.sigtmap, sigt)
                        def = fmm.sigtmap[sigt]
                        lineoffset = fmm.defmap[def][2]
                        if lineoffset != 0
                            t = StackTraces.StackFrame(t.func, t.file, t.line+lineoffset, t.linfo, t.from_c, t.inlined, t.pointer)
                            trace[i] = (t, n)
                            break
                        end
                    end
                end
            end
        end
    end
    return trace
end

function macroreplace!(ex::Expr, filename)
    for i = 1:length(ex.args)
        ex.args[i] = macroreplace!(ex.args[i], filename)
    end
    if ex.head == :macrocall
        m = ex.args[1]
        if m == Symbol("@__FILE__")
            return String(filename)
        elseif m == Symbol("@__DIR__")
            return dirname(String(filename))
        end
    end
    return ex
end
macroreplace!(s, filename) = s

"""
    steal_repl_backend(backend = Base.active_repl_backend)

Replace the REPL's normal backend with one that calls [`revise`](@ref) before executing
any REPL input.
"""
function steal_repl_backend(backend = Base.active_repl_backend)
    @async begin
        # terminate the current backend
        put!(backend.repl_channel, (nothing, -1))
        fetch(backend.backend_task)
        # restart a new backend that differs only by processing the
        # revision queue before evaluating each user input
        backend.backend_task = @async begin
            while true
                tls = task_local_storage()
                tls[:SOURCE_PATH] = nothing
                ast, show_value = take!(backend.repl_channel)
                if show_value == -1
                    # exit flag
                    break
                end
                # Process revisions
                revise(backend)
                # Now eval the input
                REPL.eval_user_input(ast, backend)
            end
        end
    end
    nothing
end

"""
    Revise.async_steal_repl_backend()

Wait for the REPL to complete its initialization, and then call [`steal_repl_backend`](@ref).
This is necessary because code registered with `atreplinit` runs before the REPL is
initialized, and there is no corresponding way to register code to run after it is complete.
"""
function async_steal_repl_backend()
    mode = get(ENV, "JULIA_REVISE", "auto")
    if mode == "auto"
        atreplinit() do repl
            @async begin
                iter = 0
                # wait for active_repl_backend to exist
                while !isdefined(Base, :active_repl_backend) && iter < 20
                    sleep(0.05)
                    iter += 1
                end
                if isdefined(Base, :active_repl_backend)
                    steal_repl_backend(Base.active_repl_backend)
                else
                    @warn "REPL initialization failed, Revise is not in automatic mode. Call `revise()` manually."
                end
            end
        end
    end
    return nothing
end

function __init__()
    # Base.register_root_module(Base.__toplevel__)
    if isfile(silencefile[])
        pkgs = readlines(silencefile[])
        for pkg in pkgs
            push!(silence_pkgs, Symbol(pkg))
        end
    end
    push!(Base.package_callbacks, watch_package)
    push!(Base.include_callbacks,
        (mod::Module, fn::AbstractString) -> push!(included_files, (mod, normpath(abspath(fn)))))
    mode = get(ENV, "JULIA_REVISE", "auto")
    if mode == "auto"
        if isdefined(Base, :active_repl_backend)
            steal_repl_backend(Base.active_repl_backend)
        elseif isdefined(Main, :IJulia)
            Main.IJulia.push_preexecute_hook(revise)
        end
        if isdefined(Main, :Atom)
            for x in ["eval", "evalall", "evalrepl"]
                old = Main.Atom.handlers[x]
                Main.Atom.handle(x) do data
                    revise()
                    old(data)
                end
            end
        end
    end
    polling = get(ENV, "JULIA_REVISE_POLL", "0")
    if polling == "1"
        polling_files[] = watching_files[] = true
    end
    rev_include = get(ENV, "JULIA_REVISE_INCLUDE", "0")
    if rev_include == "1"
        tracking_Main_includes[] = true
    end
    # Correct line numbers for code moving around
    Base.update_stackframes_callback[] = update_stacktrace_lineno!
end

## WatchList utilities
function systime()
    tv = Libc.TimeVal()
    tv.sec + tv.usec/10^6
end
function updatetime!(wl::WatchList)
    wl.timestamp = systime()
end
Base.push!(wl::WatchList, filename) = push!(wl.trackedfiles, filename)
WatchList() = WatchList(systime(), Set{String}())
Base.in(file, wl::WatchList) = in(file, wl.trackedfiles)

end # module


text:
using Printf
using KiteModels, KitePodModels, KiteUtils

if ! @isdefined kcu;  const kcu  = KCU(se()) end
if ! @isdefined kps4; const kps4 = KPS4(kcu) end
if ! @isdefined kps3; const kps3 = KPS3(kcu) end

# the following values can be changed to match your interest
dt = 0.05
ALPHA_ZERO = 8.8
STEPS = round(600/dt*0.05)
PLOT = true
FRONT_VIEW = false
ZOOM = true
PRINT = false
STATISTIC = false
# end of user parameter section #

se().version = 2

if PLOT
    using Plots
    include("plot2d.jl")
end

function simulate(s, integrator, steps, plot=false)
    start = integrator.p.iter
    for i in 1:steps
        if PRINT
            lift, drag = KiteModels.lift_drag(s)
            @printf "%.2f: " round(integrator.t, digits=2)
            println("lift, drag  [N]: $(round(lift, digits=2)), $(round(drag, digits=2))")
        end

        KiteModels.next_step!(s, integrator, dt=dt)

        if plot
            reltime = i*dt
            if mod(i, 5) == 0
                p = plot2d(s.pos, reltime; zoom=ZOOM, front=FRONT_VIEW)
                display(p)
            end
        end
    end
    (integrator.p.iter - start) / steps
end

integrator = KiteModels.init_sim!(kps3, stiffness_factor=0.04, prn=STATISTIC)
av_steps = simulate(kps3, integrator, STEPS, true)

lift, drag = KiteModels.lift_drag(kps3)
println("KPS3")
println("lift, drag  [N]: $(round(lift, digits=2)), $(round(drag, digits=2))")
println("winch_force [N]: $(round(winch_force(kps3), digits=2))")
println("Average number of callbacks per time step: $av_steps")

kps4.set.alpha_zero = ALPHA_ZERO
integrator = KiteModels.init_sim!(kps4, stiffness_factor=0.04, prn=STATISTIC)
av_steps = simulate(kps4, integrator, STEPS, true)

lift, drag = KiteModels.lift_drag(kps4)
println("KPS4")
println("lift, drag  [N]: $(round(lift, digits=2)), $(round(drag, digits=2))")
println("winch_force [N]: $(round(winch_force(kps4), digits=2))")
println("Average number of callbacks per time step: $av_steps")


text:
module Essentials

using Reexport: @reexport

include("Spatials.jl")
include("DegreesOfFreedom.jl")
include("Terms.jl")
include("FockPackage.jl")
include("SpinPackage.jl")
include("Frameworks.jl")

@reexport using .Spatials
@reexport using .DegreesOfFreedom
@reexport using .Terms
@reexport using .FockPackage
@reexport using .SpinPackage
@reexport using .Frameworks

end # module


text:
"""
    todays_plot(; region = "")

Plot the 48 forecast for the carbon intensity starting from midnight on the day that the
function is called. By default shows the national forecast. If the `region` kwarg is a valid
region, then plot data for this region only. Valid regions are contained in
[`AVAILABLE_REGIONS`](@ref).
"""
@userplot todays_plot
@recipe function f(CI::todays_plot; region = "")

    data = if isempty(region)
        get_todays_forecast()
    else
        get_todays_forecast(;regional = true, region = region)
    end

    :seriestype --> :line
    :yguide --> "Carbon Intensity (gCO₂/kWh)"
    :title --> if isempty(region) "UK Carbon Intensity" else "$region Carbon Intensity" end
    :xrotation --> 45
    :size --> (700, 500)
    :legend --> :topleft
    :margin --> 12mm

    @series begin
        :label --> "Forecast"
        :seriescolor --> :black
        :linestyle --> :dot
        :linewidth --> 2
        (data.to, data.forecast)
    end

    if isempty(region)
        # plot handles `NaN`s, not `nothing`
        data.actual = replace(data.actual, nothing=>NaN)
        @series begin
            :label --> "Actual"
            :seriescolor --> :black
            :linewidth --> 2
            :linestyle --> :solid
            (data.to, data.actual)
        end
    end

    colors = PlotUtils.palette(:RdYlGn_10, length(CARBON_INDEX_BOUNDS), rev = true)
    for (counter, (label, boundaries)) in enumerate(CARBON_INDEX_BOUNDS)
        @series begin
            :seriestype --> :line
            :label --> label
            :seriescolor --> colors[counter]
            :fillrange --> (boundaries)
            :fillalpha --> 0.2
            :seriesalpha --> 0.8
            (data.to, fill(last(boundaries), length(data.to)))
        end
    end
end

"""
    carbon_intensity(start_date, end_date)

Plots the recorded carbon intensity over the period given by the `start_date`
and `end_date`.
"""
@userplot carbon_intensity
@recipe function f(CI::carbon_intensity)
    start_date, end_date = CI.args

    data = get_carbon_intensity(start_date, end_date)

    # plot handles `NaN`s, not `nothing`
    data.actual = replace(data.actual, nothing=>NaN)

    :yguide --> "Carbon Intensity (gCO₂/kWh)"
    :title --> "UK Carbon Intensity"
    :size --> (700, 500)
    :margin --> 10mm

    @series begin
        :seriestype --> :line
        :label --> "actual"
        :seriescolor --> :black
        (data.to, data.actual)
    end

    colors = PlotUtils.palette(:RdYlGn_10, length(CARBON_INDEX_BOUNDS), rev = true)
    for (counter, (label, boundaries)) in enumerate(CARBON_INDEX_BOUNDS)
        @series begin
            :seriestype --> :line
            :label --> label
            :seriescolor --> colors[counter]
            :fillrange --> (boundaries)
            :fillalpha --> 0.2
            :seriesalpha --> 0.8
            (data.to, fill(last(boundaries), length(data.to)))
        end
    end
end


text:
# # Code:
# abstract Robot
#
# abstract Manipulator <: Robot
#
# abstract Movile <: Robot

#
# robots_types.jl ends here

max_stars_repo_name: elsuizo/Robot

max_issues_repo_name: elsuizo/Robot

max_forks_repo_name: elsuizo/Robot