text:
#  _____ _____ ____ _____ ____        ____        ____  ___   ___  ____
# |_   _| ____/ ___|_   _/ ___|      |  _ \      / ___|/ _ \ / _ \|  _ \
#   | | |  _| \___ \ | | \___ \ _____| |_) |____| |  _| | | | | | | | | |
#   | | | |___ ___) || |  ___) |_____|  _ <_____| |_| | |_| | |_| | |_| |
#   |_| |_____|____/ |_| |____/      |_| \_\     \____|\___/ \___/|____/

#                _       _ _
#  ___  ___ _ __(_) __ _| (_)_______
# / __|/ _ \ '__| |/ _` | | |_  / _ \
# \__ \  __/ |  | | (_| | | |/ /  __/
# |___/\___|_|  |_|\__,_|_|_/___\___|

@test DynamoDB.null_or_val(nothing) == Dict("NULL" => true)

@test DynamoDB.null_or_val(true) == Dict("BOOL" => true)
@test DynamoDB.null_or_val(false) == Dict("BOOL" => false)

@test DynamoDB.null_or_val("foo") == Dict("S" => "foo")

@test DynamoDB.null_or_val(21) == Dict("N" => "21")
@test DynamoDB.null_or_val(21.1) == Dict("N" => "21.1")

@test DynamoDB.null_or_val([1, 2, 3]) == Dict("L" => [Dict("N" => "1"), Dict("N" => "2"), Dict("N" => "3")])

# result lists is (unsurprisingly) unordered...
res = DynamoDB.null_or_val(Set([1, 2, 3]))
@test length(keys(res)) == 1
@test length(res["NS"]) == 3
@test contains(==, res["NS"], "1")
@test contains(==, res["NS"], "2")
@test contains(==, res["NS"], "3")

res = DynamoDB.null_or_val(Set(["1", "2", "3"]))
@test length(keys(res)) == 1
@test length(res["SS"]) == 3
@test contains(==, res["SS"], "1")
@test contains(==, res["SS"], "2")
@test contains(==, res["SS"], "3")

@test DynamoDB.null_or_val(Dict("four" => 3, 7 => 11)) == Dict("M"=>Dict("four"=>Dict("N"=>"3"),"7"=>Dict("N"=>"11")))

# see runtests.jl for Foo's definition

@test DynamoDB.null_or_val(Foo(3, "fourty")) == Dict("M"=>Dict("a"=>Dict("N"=>"3"),"b"=>Dict("S"=>"fourty")))

#      _                     _       _ _
#   __| | ___  ___  ___ _ __(_) __ _| (_)_______
#  / _` |/ _ \/ __|/ _ \ '__| |/ _` | | |_  / _ \
# | (_| |  __/\__ \  __/ |  | | (_| | | |/ /  __/
#  \__,_|\___||___/\___|_|  |_|\__,_|_|_/___\___|

function check_round_trip(val)
    @test DynamoDB.value_from_attributes(DynamoDB.null_or_val(val)) == val
end

@test DynamoDB.real_val(21.0) == 21.0
@test DynamoDB.bool_val("true") == true
@test DynamoDB.bool_val("false") == false
@test_throws ErrorException DynamoDB.bool_val("not a bool")

@test_throws ErrorException DynamoDB.value_from_attributes(Dict("M" => 1, "L" => 2)) # two type decls
@test_throws ErrorException DynamoDB.value_from_attributes(Dict("pants" => "frowny")) # unknown type decl

# base types
check_round_trip(nothing)
check_round_trip(true)
check_round_trip(false)
check_round_trip("foo")
check_round_trip("false")
check_round_trip("true")
check_round_trip("nothing")
check_round_trip("21")
check_round_trip(21)
check_round_trip(21.0)

# compound types
check_round_trip([1, 2, 3])
check_round_trip(Set([1, 2, 3]))
check_round_trip(Set(["1", "2", "3"]))

@test DynamoDB.value_from_attributes(DynamoDB.null_or_val(Dict("four" => 3, 7 => 11))) == Dict("four" => 3, "7" => 11)
@test DynamoDB.value_from_attributes(DynamoDB.null_or_val(Foo(3, "fourty"))) == Dict("a"=>3,"b"=>"fourty")

res = DynamoDB.value_from_attributes(Foo, DynamoDB.null_or_val(Foo(3, "fourty")))
@test res.a == 3
@test res.b == "fourty"


text:
@enum VType realv intv

"""
Random variable.

``Var`` is a parametric type with a type parameter N, which indicates
the number of dimensions for each sample.

Each instance of ``Var`` represents a random variable that can take
different values in different trials. The instance contains four
public fields:

- ``id::String``:       The identifier of the variable. In a joint model,
                        each variable should have a unique id.
- ``space::S``:         The underlying value space (for each element).
- ``dims::Dims{N}``:    The dimensions of each sample.
- ``len::Int``:         The number of elements in each sample
                        (i.e. ``prod(size)``).
"""
mutable struct Var{N}
    id::String          # identifier
    vtype::VType        # variable type (realv or intv)
    card::Int           # cardinality (0 when it is infinity)
    dims::Dims{N}       # dimensions for each sample
    len::Int            # number of elements of each sample
end

"""Scalar random variable."""
const SVar = Var{0}

dvar(id::String, c::Int) = Var{0}(id, intv, c, (), 1)
dvar(id::String, c::Int, n::Int) = Var{1}(id, intv, c, (n,), n)
dvar{N}(id::String, c::Int, dims::Dims{N}) = Var{N}(id, intv, c, dims, prod(dims))

rvar(id::String) = Var{0}(id, realv, 0, (), 1)
rvar(id::String, n::Int) = Var{1}(id, realv, 0, (n,), n)
rvar{N}(id::String, dims::Dims{N}) = Var{N}(id, realv, 0, dims, prod(dims))

ndims{N}(v::Var{N}) = N
size(v::Var) = v.dims
length(v::Var) = v.len
vtype(v::Var) = v.vtype
cardinality(v::Var) = v.card

==(x::Var, y::Var) = false
=={N}(x::Var{N}, y::Var{N}) =
    (x.id == y.id && x.vtype == y.vtype && x.card == y.card && x.dims == y.dims)
!=(x::Var, y::Var) = !(x == y)

function show(io::IO, v::Var)
    print(io, "Var($(v.id))[vtype=$(v.vtype), dims=$(v.dims)]")
end

const VarList{N} = NamedList{Var{N}}
const SVarList = VarList{0}

varlist{N}(vars::AbstractVector{Var{N}}) =
    namedlist(vars, v->v.id)::VarList{N}

ext: jl


text:
using TypeDB_tutorial
using Test

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


text:
using ProgressMeter
using Base.Threads

"""
    filter_positions(snap_file::String, corner_lowerleft::Array{<:Real}, corner_upperright::Array{<:Real},
                          parttype::Integer)

Reads positions from `snap_file` and returns the indices of particles contained in a box defined by `corner_lowerleft` and `corner_upperright`.
"""
function filter_positions(snap_file::String, corner_lowerleft::Array{<:Real}, corner_upperright::Array{<:Real},
                          parttype::Integer)

    pos = read_snap(snap_file, "POS", parttype)
    sel = findall( @views ( ( corner_lowerleft[1] .<= pos[1,:] .< corner_upperright[1] ) .&
                            ( corner_lowerleft[2] .<= pos[2,:] .< corner_upperright[2] ) .&
                            ( corner_lowerleft[3] .<= pos[3,:] .< corner_upperright[3] )) )
    return sel
end

"""
    filter_ids(snap_file::String, selected_ids::Vector{<:Integer})

Reads IDs from `snap_file` and returns the indices of particles matching the `selected_ids`
"""
function filter_ids(snap_file::String, selected_ids::Vector{<:Integer}, parttype::Integer)

    # read ids in file
    ids = read_snap(snap_file, "ID", parttype)

    # check which IDs match
    return filter_by_ids( selected_ids, ids)
end

"""
    read_blocks_over_all_files( snap_base::String, blocks::Array{String};
                                filter_function::Union{Function, Nothing}=nothing,
                                read_positions::Union{Dict, Nothing}=nothing,
                                parttype::Integer=0, verbose::Bool=true )

Reads the specified blocks from all distributed files where particles pass the `filter_function`, or are given by a `Dict` of `read_positions`.
For `read_positions` please see [`find_read_positions`](@ref).
"""
function read_blocks_over_all_files(snap_base::String, blocks::Array{String};
                                    filter_function::Union{Function, Nothing}=nothing,
                                    read_positions::Union{Dict, Nothing}=nothing,
                                    parttype::Integer=0, verbose::Bool=true )

    # default behaviour is to find the particles to read.
    # If this number is known in advance it can be given as an input parameter.
    if isnothing(read_positions)

        # this only works if a filter function is provided!
        if isnothing(filter_function)
            error("Please provide either a dictionary with read positions or a filter function!")
        end

        # check how many particles fulfill the filter criterion and find the positions in the blocks.
        # -> needed for array pre-allocation to speed up read-in.
        read_positions = find_read_positions(snap_base, filter_function, verbose=verbose)
    end

    if verbose
        @info "Allocating storage dictionary..."
        t1 = time_ns()
    end

    # select current file
    filename = select_file(snap_base, 0 )

    # read info block
    snap_info = read_info(filename)

    # pre-allocate all data arrays in a dictionary
    d = allocate_data_dict(blocks, read_positions["N_part"], snap_info, false)

    if verbose
        t2 = time_ns()
        @info "  elapsed: $(output_time(t1,t2)) s"
    end

    # store the number of particles that have been read
    N_read = 0

    # get the relevant files from the dictionary
    key_entries = collect(keys(read_positions))
    files = key_entries[ key_entries .!= "N_part"]

    sort!(files)

    if verbose
        @info "Reading $(size(files,1)) snapshots..."
        t1 = time_ns()
    end

    for file ∈ files

        # select current file
        filename = select_file(snap_base, file )

        # read header block
        h = read_header(filename)

        # read info block
        snap_info = read_info(filename)

        # find the location of the blocks in the current file
        file_block_positions = get_block_positions(filename)

        # read the blocks
        @threads for block ∈ blocks

            # read the info ot the curent block
            block_info = snap_info[getfield.(snap_info, :block_name) .== block][1]

            # add offset of particle types that should not be read
            offset = 0
            for i = 1:size(h.npart,1)
                if (block_info.is_present[i] > 0) & (h.npart[i] > 0) & ( i < parttype + 1)
                    offset += h.npart[i]
                end
            end

            # store the relevant data in the `d[block]` array
            read_block_with_offset!(d[block], N_read, filename,
                                    file_block_positions[block],
                                    block_info, offset,
                                    read_positions[file]["index"],
                                    read_positions[file]["n_to_read"])

        end
        # count the number of particles already read
        N_read += sum(read_positions[file]["n_to_read"])

        if verbose
            @info "Read $N_read / $(read_positions["N_part"]) particles"
        end

    end

    if verbose
        t2 = time_ns()
        @info "  elapsed: $(output_time(t1,t2)) s"
    end

    return d
end


text:
verbose = true
 if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 using RvSpectMLBase
 using EchelleInstruments, EchelleInstruments.NEID
 using RvSpectML
 if verbose println("# Loading NeidSolarScripts")    end
 using NeidSolarScripts
 #using NeidSolarScripts.SolarRotation
 #using NeidSolarScripts.DifferentialExtinction
 if verbose   println("# Loading other packages")    end
 using CSV, DataFrames, Query, Dates
 using JLD2, FileIO
 #using StatsBase, Statistics, Dates

fits_target_str = "Sun"
  #target_subdir = "/mnt/data_simons/NEID/DRPv0.7-fixedflatfielding2"   # USER: Replace with directory of your choice
  #output_dir = "output/"
  target_subdir = "/gpfs/group/ebf11/default/ebf11/neid_solar/data/outputs/20210330"
  output_dir = ""
 if !@isdefined(target_subdir)
      target_subdir = ""
 end
 if !@isdefined(create_missing_continuum_files)
     create_missing_continuum_files = false
 end
 paths_to_search_for_param = [pwd(),joinpath(dirname(pathof(NeidSolarScripts)),"..")]

 if verbose println("# Finding what data files are avaliable.")  end
 eval(read_data_paths(paths_to_search=paths_to_search_for_param))
 @assert isdefined(Main,:neid_data_path)
 @assert isdefined(Main,:output_dir)
 output_dir = ""
 output_path = joinpath(neid_data_path, target_subdir,output_dir)
 manifest_filename = joinpath(output_path,"manifest.csv")
 manifest_calib_filename = joinpath(output_path,"manifest_calib.csv")

if isfile(manifest_filename)
    df_files  = CSV.read(manifest_filename, DataFrame)
    @assert size(df_files,1) >= 1
    @assert hasproperty(df_files,:Filename)
    @assert hasproperty(df_files,:target)
    @assert hasproperty(df_files,:bjd)
    @assert hasproperty(df_files,:ssbz)
    @assert hasproperty(df_files,:exptime)
    @assert hasproperty(df_files,:alt_sun)
    @assert hasproperty(df_files,:Δfwhm²)
    @assert hasproperty(df_files,:order_snrs)
    if eltype(df_files[!,:order_snrs]) == String
        df_files[!,:order_snrs] = map(i->parse.(Float64,split(df_files[i,:order_snrs][2:end-1],',')),1:size(df_files,1))
    end
    @assert eltype(df_files[!,:order_snrs]) == Vector{Float64}
else
    @error "Need to make manifest at $manifest_filename"
    #df_files = NEID.make_manifest(joinpath(neid_data_path, target_subdir))
    #CSV.write(manifest_filename, df_files)
end

paths_to_search_for_param = [pwd(),"..","examples",joinpath(pkgdir(RvSpectMLBase),"..","RvSpectML","examples"), "/gpfs/scratch/jpn23/"]
  # NOTE: make_manifest does not update its paths_to_search when default_paths_to_search is defined here, so if you change the line above, you must also include "paths_to_search=default_paths_to_search" in the make_manifest() function call below
  #outputs = Dict{String,Any}()
  pipeline_plan = PipelinePlan()
  dont_make_plot!(pipeline_plan, :movie)
  if verbose println("# Reading in customized parameters from param.jl.")  end
   if !@isdefined(idx_day_to_use)
       idx_day_to_use = 1
   end
   eval(code_to_include_param_jl(paths_to_search=paths_to_search_for_param))
   #outputs["df_files_use"] = df_files_use

if match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1] ==  match(r"neidL1_(\d+)[T_](\d+)\.fits$", last(df_files_use.Filename))[1]
      date_str = match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1]
    else
      date_str = string(match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1]) * "-" * string(match(r"neidL1_(\d+)[T_](\d+)\.fits$", last(df_files_use.Filename))[1])
   end

#=
   outputs_filename = joinpath(output_dir,"solar_" * date_str * "_new.jld2")
   if isfile(outputs_filename) && false
     times_already_processed = load(outputs_filename, "times")
     files_in_day_to_process = size(df_files_solar_by_day.data[idx_day_to_use],1)
      if files_in_day_to_process == length(times_already_processed)
         println("# Already processed all ", length(times_already_processed), " files for ", date_str)
         exit()
      end
   end
=#

#for idx_day_to_use in 1:size(df_files_solar_by_day,1)
   #df_files_use = df_files_solar_by_day[idx_day_to_use,:data] |> @orderby(_.bjd) |> @take(max_spectra_to_use) |> DataFrame

#=
if verbose println("# Reading in ", size(df_files_use,1), " FITS files for ", date_str, ".")  end
    @time all_spectra = map(row->NEID.read_data(row), eachrow(df_files_use))
    GC.gc()
    dont_need_to!(pipeline_plan,:read_spectra)
=#

all_spectra = Spectra2DBasic{Float64, Float32, Float32, Matrix{Float64}, Matrix{Float32}, Matrix{Float32}, NEID2D}[]
for (i,row) in enumerate(eachrow(df_files_use))
    m = match(r"(neidL1_\d+[T_]\d+)\.fits$", row.Filename)
    #continuum_filename = joinpath(neid_data_path,target_subdir,"output","continuum", m[1] * "_continuum=afs.jld2")
    continuum_filename = joinpath(neid_data_path,target_subdir,output_dir,"continuum", m[1] * "_continuum=afs.jld2")
    if !(isfile(continuum_filename)||islink(continuum_filename))
        println("# Couldn't find matching continuum file", continuum_filename)
        continue
    end
    spec = NEID.read_data(row)
    continuum = load(continuum_filename, "continuum")
    spec.flux ./= continuum
    spec.var ./= continuum.^2
    push!(all_spectra,spec)
    #if i>=20  break   end
end
GC.gc()
dont_need_to!(pipeline_plan,:read_spectra)

line_width_50_default = 7.9e3
 lsf_width = 3.0e3
 max_mask_scale_factor = max(lsf_width,line_width_50_default/sqrt(8*log(2)))/default_ccf_mask_v_width(NEID2D()) # 4.0
 max_bc = RvSpectMLBase.max_bc
 #max_bc = RvSpectMLBase.max_bc_earth_rotation
 line_list_filename = joinpath("data","solar_line_list_espresso.csv")
 max_orders = min_order(NEID2D()):max_order(NEID2D())
 good_orders = orders_to_use_default(NEID2D())
 orders_to_use = max_orders
 if isfile(line_list_filename)
   println("# Reading ", line_list_filename)
   line_list_espresso = CSV.read(line_list_filename, DataFrame)
   dont_need_to!(pipeline_plan,:clean_line_list_tellurics)
 else
    println("# Can't find ", line_list_filename, ".  Trying ESPRESSO line list.")
    orders_to_use = good_orders
    order_list_timeseries = extract_orders(all_spectra,pipeline_plan, orders_to_use=orders_to_use, recalc=true )

    line_list_filename = "G2.espresso.mas"
    linelist_for_ccf_fn_w_path = joinpath(pkgdir(EchelleCCFs),"data","masks",line_list_filename)

    line_list_espresso = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan, v_center_to_avoid_tellurics=ccf_mid_velocity,
       Δv_to_avoid_tellurics = 2*max_bc+5*line_width_50_default+max_mask_scale_factor*default_ccf_mask_v_width(NEID2D()), orders_to_use=orders_to_use, recalc=true, verbose=true)

    #CSV.write(custom_line_list_filename, line_list_espresso)
 end
 #outputs["line_list_espresso"] = line_list_espresso

line_width_50 = line_width_50_default
 maxΔfwhm² = -0.569375
 @assert maximum(df_files_use.Δfwhm²) < maxΔfwhm²
 Δfwhm = 1000.0 .* sqrt.(maxΔfwhm².-df_files_use.Δfwhm²[1:length(all_spectra)])  # How much to increase fwhm by to acheive uniform fwhm

msf = lsf_width/default_ccf_mask_v_width(NEID2D()); fwtf = 0.5  # using LSF width
 order_list_timeseries = extract_orders(all_spectra,pipeline_plan,  orders_to_use=orders_to_use,  remove_bad_chunks=false, recalc=true )

@time (order_ccfs, order_ccf_vars, v_grid_order_ccfs) = ccf_orders(order_list_timeseries, line_list_espresso, pipeline_plan,
    mask_type=:gaussian, Δfwhm=Δfwhm,
    mask_scale_factor=msf, range_no_mask_change=5*line_width_50, ccf_mid_velocity=ccf_mid_velocity, v_step=155,
    v_max=max(5*line_width_50,2*max_bc), orders_to_use=orders_to_use, allow_nans=true, calc_ccf_var=true,
    recalc=true)

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

ccf_dir = joinpath(neid_data_path,target_subdir,"ccfs")
if !isdir(ccf_dir)
   mkdir(ccf_dir)
end

#daily_ccf_filename = joinpath(neid_data_path,target_subdir,"output","ccfs", date_str * "_ccfs=default.jld2")
daily_ccf_filename = joinpath(neid_data_path,target_subdir,"ccfs", date_str * "_ccfs=default.jld2")
  jldopen(daily_ccf_filename, "w") do f
    f["v_grid"] = v_grid_order_ccfs
    f["order_ccfs"] = order_ccfs
    f["order_ccf_vars"] = order_ccf_vars
    f["orders_to_use"] = orders_to_use
    f["manifest"] = df_files_use
    f["daily_ccf_filename"] = daily_ccf_filename
  end

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

#=

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

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

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

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

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

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

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

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

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

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

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

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

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

    end
end
=#

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

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

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

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

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

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

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

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

 outputs["mean_ccf"] = ccf_template

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

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

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

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

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

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

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

##

good_orders = orders_to_use_default(NEID2D())

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

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

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

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

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

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

save(outputs_filename, outputs)
=#
#=

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


text:
using LinearAlgebra

@testset "Test the tensor network utlities" begin
    # test single qubit gates
    h = [[1., 1.] [1., -1.]]
    @test QXTns.find_hyper_edges(h) == []
    z = [[1., 0.] [0., -1.]]
    @test QXTns.find_hyper_edges(z) == [[1, 2]]

    # test 2 qubit identity
    id = collect(Diagonal(ones(4)))
    @test QXTns.find_hyper_edges(id) == [[1, 2]]
    @test QXTns.find_hyper_edges(reshape(id, (2, 2, 2, 2))) == [[1, 3], [2, 4]]
    B, C = QXTns.decompose_gate(reshape(id, (2, 2, 2, 2)))
    @test size(B)[3] == 1
    @test size(C)[1] == 1
    @test QXTns.find_hyper_edges(B) == [[1, 2]]
    @test QXTns.find_hyper_edges(C) == [[2, 3]]

    # test 2 qubit cz
    id = collect(Diagonal(ones(4)))
    @test QXTns.find_hyper_edges(id) == [[1, 2]]
    @test QXTns.find_hyper_edges(reshape(id, (2, 2, 2, 2))) == [[1, 3], [2, 4]]
    B, C = QXTns.decompose_gate(reshape(id, (2, 2, 2, 2)))
    @test size(B)[3] == 1
    @test size(C)[1] == 1
    @test QXTns.find_hyper_edges(B) == [[1, 2]]
    @test QXTns.find_hyper_edges(C) == [[2, 3]]
end

@testset "Test reduce tensor function" begin
    # first a 2d diagonal tensor
    A = rand(5)
    Ar = QXTns.reduce_tensor(Diagonal(A), [[1, 2]])
    @test Ar == A

    # create a rank 5 tensor from a rank 3 tensor and reduce again
    A = rand(4, 5, 6)
    Afull = zeros(4, 5, 4, 5, 6)
    for i in CartesianIndices(size(A))
        Afull[i[1], i[2], i[1], i[2], i[3]] = A[i]
    end
    Ar = QXTns.reduce_tensor(Afull, [[1, 3], [2, 4]])
    @test Ar == A
end


text:
using UtilitiesForMRI, LinearAlgebra, Test, Random, AbstractLinearOperators

# Cartesian domain
n = (256,256,256)
h = (abs(randn()), abs(randn()), abs(randn()))
X = spatial_sampling(Float64, n; h=h)

# Cartesian sampling in k-space
phase_encoding = (1,2)
subsampling = (1:256^2)[randperm(256^2)][1:32]
K = kspace_Cartesian_sampling(X; phase_encoding=phase_encoding, subsampling=subsampling)

# Fourier operator
tol = 1e-6
F = nfft(X, K; tol=tol)

# Adjoint test (linear operator)
nt, nk = size(K)
θ = 1e-3*pi*randn(Float64, nt, 6)
Fθ = F(θ)
u = randn(ComplexF64, n)
d = randn(ComplexF64, nt, nk)
@test dot(Fθ*u, d) ≈ dot(u, Fθ'*d) rtol=1e-6

# Jacobian & consistency test
d, _, ∂Fθu = ∂(F()*u, θ)
@test d ≈ Fθ*u rtol=1e-6

# Adjoint test (Jacobian)
Δθ = randn(Float64, nt, 6); Δθ *= norm(θ)/norm(Δθ)
ΔF_ = ∂Fθu*Δθ; ΔF = randn(ComplexF64, size(ΔF_)); ΔF *= norm(ΔF_)/norm(ΔF)
@test dot(∂Fθu*Δθ, ΔF) ≈ dot(Δθ, ∂Fθu'*ΔF) rtol=1e-6

# Gradient test
Δθ = randn(Float64, nt, 6); Δθ *= norm(θ)/norm(Δθ)
t = 1e-6
Fθu_p1 = F(θ+0.5*t*Δθ)*u
Fθu_m1 = F(θ-0.5*t*Δθ)*u
@test (Fθu_p1-Fθu_m1)/t ≈ ∂Fθu*Δθ rtol=1e-3

# Gauss-Newton Hessian
w = randn(ComplexF64, size(ΔF))
W = linear_operator(ComplexF64, size(ΔF), size(ΔF), d->w.*d, d-> conj(w).*d)
h = randn(ComplexF64, size(ΔF))
H = linear_operator(ComplexF64, size(ΔF), size(ΔF), d->h.*d, d-> conj(h).*d)
HF = sparse_matrix_GaussNewton(∂Fθu; W=W, H=H)

# Consistency test
Δθ = randn(Float64, nt, 6)
@test HF*vec(Δθ) ≈ vec(real(∂Fθu'*W'*(H*(W*(∂Fθu*Δθ))))) rtol=1e-6

max_stars_repo_name: grizzuti/UtilitiesForMRI

max_issues_repo_name: grizzuti/UtilitiesForMRI

max_forks_repo_name: grizzuti/UtilitiesForMRI


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Compute the precession as in equinox-based IAU-2006 theory.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [2] Wallace, P. T., Capitaine, N (2006). Precession-nutation procedures
#       consistent with IAU 2006 resolutions. Astronomy & Astrophysics.
#
#   [3] IERS (2010). Transformation between the International Terrestrial
#       Reference System and the Geocentric Celestial Reference System. IERS
#       Technical Note No. 36, Chapter 5.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export precession_iau2006

"""
    precession_iau2006(jd_tt::Number)

Compute the precession angles [rad] according to equinox-based IAU-2006 theory
in the Julia day `jd_tt` [Terrestrial Time].

This algorithm was obtained from **[1]**(p. 49).

# References

- **[1]**: IERS (2010). Transformation between the International Terrestrial
    Reference System and the Geocentric Celestial Reference System. IERS
    Technical Note No. 36, Chapter 5.
"""
function precession_iau2006(jd_tt::Number)
    # Compute the Julian Centuries from `jd_tt`.
    t_tt = (jd_tt - JD_J2000) / 36525

    # Auxiliary variables
    # ===================

    a2d = 1 / 3600
    d2r = π / 180
    a2r = a2d * d2r

    Ψ_a = @evalpoly(
        t_tt,
        0,
        +5038.481_507,
        -1.079_006_9,
        -0.001_140_45,
        +1.328_51e-4,
        -9.51e-8
    )
    Ψ_a = mod(Ψ_a * a2r, 2π)

    ω_a = @evalpoly(
        t_tt,
        +84381.406,
        -0.025_754,
        +0.051_262_3,
        -0.007_725_03,
        -4.67e-7,
        +3.337e-7
    )
    ω_a = mod(ω_a * a2r, 2π)

    χ_a = @evalpoly(
        t_tt,
        0,
        +10.556_403,
        -2.381_429_2,
        -0.001_211_97,
        +1.706_63e-4,
        -5.60e-8
    )
    χ_a = mod(χ_a * a2r, 2π)

    return Ψ_a, ω_a, χ_a
end


text:
using Test
using CrypticCrosswords

function solve_and_explain(clue, length)
    solutions, state = solve(clue, length=length)
    derivations = derive!(state, solutions[1])
    io = IOBuffer()
    explain(io, first(derivations))
    String(take!(io))
end

@testset "Explanations" begin

    @test solve_and_explain("spin broken shingle", 7) == """
The answer is "english".
"spin" is the definition.
"broken" means to anagram "shingle" to get "english".
"english" matches "spin" with confidence 100%.
"""

    @test solve_and_explain("initially babies are naked", 4) == """
The answer is "bare".
"initially" means to take the first letter of "babies" to get "b".
Combine "b" and "are" to get "bare".
"naked" is the definition.
"bare" matches "naked" with confidence 100%.
"""

    @test solve_and_explain("hungary's leader, stuffy and bald", 8) == """
The answer is "hairless".
"leader" means to take the first letter of "hungarys" to get "h".
Take a synonym of "stuffy" to get "airless".
Combine "h" and "airless" to get "hairless".
"bald" is the definition.
"hairless" matches "bald" with confidence 100%.
"""

    @test solve_and_explain("aerial worker anne on the way up", 7) == """
The answer is "antenna".
"aerial" is the definition.
A common replacement for "worker" is "ant".
"on the way up" means to reverse "anne" to get "enna".
Combine "ant" and "enna" to get "antenna".
"antenna" matches "aerial" with confidence 100%.
"""

end


text:
using EmpiricalPseudopotential
using Test

const p = EmpiricalPseudopotential

@testset "Bulk InSb" begin
    InSb = p.Material("InSb", "fcc")
    band = p.EigenEnergy(InSb, 500)

    p.BandPlot(band, title="InSb")

    gap = p.BandGap(band)
    mass = p.EffectiveMass(band, 200, 300)
end

@testset "Strain InSb" begin
    InSb = p.Material("InSb", "fcc")
    p.strain(InSb, percent=0.5)
    band = p.EigenEnergy(InSb, 500)

    p.BandPlot(band, title="Strained InSb")

    gap = p.BandGap(band)
    mass = p.EffectiveMass(band, 200, 300)
end

@testset "GaP" begin
    GaP = p.Material("GaP", "fcc")
    band = p.EigenEnergy(GaP, 500)

    p.BandPlot(band, title="GaP")

    gap = p.BandGap(band)
    mass = p.EffectiveMass(band, 200, 300)
end

@testset "GaAs" begin
    GaAs = p.Material("GaAs", "fcc")

    band = p.EigenEnergy(GaAs, 500)

    p.BandPlot(band, title="GaAs")

    gap = p.BandGap(band)
    mass = p.EffectiveMass(band, 200, 300)
end

@testset "Strain Ratio" begin
    InSb = p.Material("InSb", "fcc")
    emassies = Array{Float64, 1}()
    ratio = Array{Float64, 1}()
    for rate in -3.0:0.5:3.0
        p.strain(InSb, percent=rate)
        band = p.EigenEnergy(InSb, 500)
        gap = p.BandGap(band)
        push!(emassies, p.EffectiveMass(band, 200, 300))
        push!(ratio, rate)
    end
    p.Plot(ratio, emassies)
end

@testset "Ternary Material" begin
    InSb = p.Material("InSb", "fcc")
    GaSb = p.Material("GaSb", "fcc")

    InGaSb = p.mix("InGaSb", InSb, 0.4, GaSb, 0.6)

    band = p.EigenEnergy(InGaSb, 500)

    p.BandPlot(band, title="InGaSb")

    gap = p.BandGap(band)
    mass = p.EffectiveMass(band, 200, 300)
end


text:
using Pkg
Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
using MLCourse, CSV, MLJ, DataFrames, MLJLinearModels, Random, Flux, MLJFlux, MLJXGBoostInterface, MLJDecisionTreeInterface
include("./Data_Processing.jl")

#Load the data
drop_std = generate(option = "drop", std = "true", valid = "false", test = "true");
med_std = generate(option = "med", std = "true", valid = "false", test = "true");
drop = generate(option = "drop", std = "false", valid = "false", test = "true");
med = generate(option = "med", std = "false", valid = "false", test = "true");

#Function to tune the model with hyper-parameters eta, num_round, and max_depth
function TunedModel_XGB(data)
    Random.seed!(2711)
    model = XGBoostClassifier()
    self_tuning_model = TunedModel(model = model,
                                resampling = CV(nfolds = 10),
                                tuning = Grid(goal = 100),
                                range = [range(model, :eta,
                                        lower = 1e-4, upper = .1, scale = :log),
                                 range(model, :num_round, lower = 50, upper = 500),
                                 range(model, :max_depth, lower = 2, upper = 6)],
                                measure = auc)
    self_tuning_mach = machine(self_tuning_model,
                            select(data, Not(:precipitation_nextday)),
                            data.precipitation_nextday) |> fit!
    self_tuning_mach
end

#Test the tuned model on the drop data set and med data set with both standardized and not-standardized data:
rep1 = report(TunedModel_XGB(drop_std.train));
rep1.best_history_entry.measurement

rep2 = report(TunedModel_XGB(med_std.train));
rep2.best_history_entry.measurement #AUC for med, standardized data 0.93; best hyper-parameters being eta = 0.1, num_round = 500, max_depth = 6

rep3 = report(TunedModel_XGB(drop.train));
rep3.best_history_entry.measurement #AUC for dropped, non standardized data 0.9176; best hyper-parameters being eta = 0.1, num_round = 400, max_depth = 4

rep4 = report(TunedModel_XGB(med.train));
rep4.best_history_entry.measurement #AUC for med, non standardized data 0.9297; best hyper-parameters being eta = 0.1, num_round = 310, max_depth = 6

#Test the best model over the test set with the best hyper-parameters to pick the best one:
best_mach1 = machine(XGBoostClassifier(eta = rep1.best_model.eta, num_round = rep1.best_model.num_round, max_depth = rep1.best_model.max_depth), select(drop_std.train[:,:], Not(:precipitation_nextday)), drop_std.train.precipitation_nextday)|> fit!
best_mach2 = machine(XGBoostClassifier(eta = rep2.best_model.eta, num_round = rep2.best_model.num_round, max_depth = rep2.best_model.max_depth), select(med_std.train[:,:], Not(:precipitation_nextday)), med_std.train.precipitation_nextday)|> fit!
best_mach3 = machine(XGBoostClassifier(eta = rep3.best_model.eta, num_round = rep3.best_model.num_round, max_depth = rep3.best_model.max_depth), select(drop.train[:,:], Not(:precipitation_nextday)), drop.train.precipitation_nextday)|> fit!
best_mach4 = machine(XGBoostClassifier(eta = rep4.best_model.eta, num_round = rep4.best_model.num_round, max_depth = rep4.best_model.max_depth), select(med.train[:,:], Not(:precipitation_nextday)), med.train.precipitation_nextday)|> fit!
losses(best_mach1, select(drop_std.test, Not(:precipitation_nextday)), drop_std.test.precipitation_nextday)
losses(best_mach2, select(med_std.test, Not(:precipitation_nextday)), med_std.test.precipitation_nextday)
losses(best_mach3, select(drop.test, Not(:precipitation_nextday)), drop.test.precipitation_nextday)
losses(best_mach4, select(med.test, Not(:precipitation_nextday)), med.test.precipitation_nextday)

#med with standardized data seems to perform the best

#Write in the submission file with the best machine trained on all data
train, test = generate(option = "med", std = "true", valid = "false", test = "false");
best_mach = machine(XGBoostClassifier(eta = rep2.best_model.eta, num_round = rep2.best_model.num_round, max_depth = rep2.best_model.max_depth), select(train[:,:], Not(:precipitation_nextday)), train.precipitation_nextday)|> fit!
sample = CSV.read(joinpath(@__DIR__, "data", "sample_submission.csv"), DataFrame);
pred = pdf.(predict(best_mach, test), true);
sample.precipitation_nextday = pred;
CSV.write(joinpath(@__DIR__,"results","XGBoost_submission.csv"), sample);

max_stars_repo_name: clara-rossignol/Will-it-rain-tomorrow

max_issues_repo_name: clara-rossignol/Will-it-rain-tomorrow

max_forks_repo_name: clara-rossignol/Will-it-rain-tomorrow


text:
module BFloat

import Base: sign_mask,
             exponent_mask, exponent_one, exponent_half, exponent_bits,
             exponent_bias, exponent_max, exponent_raw_max,
             significand_mask, significand_bits,
             uinttype,
             promote_rule, widen,
             Integer, Bool, Float16, Float32, Float64,
             unsafe_trunc,
             trunc, floor, ceil, round,
             -, +, *, /, ^, \, muladd,
             div, fld, cld, rem, mod,
             ==, <, <=, isless,
             abs, uabs,
             isnan, isfinite, precision,
             nextfloat, prevfloat,
             issubnormal,
             truncmask, truncbits,
             typemin, typemax, floatmin, floatmax, eps,
             bswap,
             iszero,
             reinterpret, bitstring, show, print

import Core.Intrinsics: bitcast, bswap_int

export
    # Constructors
    BFloat16,
    # Constants
    InfB16,
    NaNB16,
    # BFloat16 properties
    sign_mask,
    exponent_mask,
    exponent_one,
    exponent_half,
    significand_mask,
    significand_bits,
    exponent_bits,
    exponent_bias,
    exponent_max,
    exponent_raw_max,
    uinttype,
    # Conversions
    Float16,
    Float32,
    Float64,
    Bool,
    Integer,
    reinterpret,
    # Promotions
    promote_rule,
    widen,
    # Functions
    unsafe_trunc,
    trunc, floor, ceil, round,
    -, +, *, /, ^, \, muladd,
    div, fld, cld, rem, mod,
    ==, <, <=, isless,
    abs,
    isnan, isfinite, precision,
    nextfloat, prevfloat,
    issubnormal,
    truncmask, truncbits,
    typemin, typemax, floatmin, floatmax, eps,
    bswap,
    iszero,
    # Representation
    bitstring, show, print,
    # Utils
    get_scm

primitive type BFloat16 <: AbstractFloat 16 end

const InfB16 = bitcast(BFloat16, 0x7f80)
const NaNB16 = bitcast(BFloat16, 0x7fc0)

uinttype(::Type{BFloat16}) = UInt16

sign_mask(::Type{BFloat16}) =        0x8000
exponent_mask(::Type{BFloat16}) =    0x7f80
exponent_one(::Type{BFloat16}) =     0x3f80
exponent_half(::Type{BFloat16}) =    0x3f00
significand_mask(::Type{BFloat16}) = 0x007f

significand_bits(::Type{BFloat16}) = trailing_ones(significand_mask(BFloat16))
exponent_bits(::Type{BFloat16}) = sizeof(BFloat16)*8 - significand_bits(BFloat16) - 1
exponent_bias(::Type{BFloat16}) = Int(exponent_one(BFloat16) >> significand_bits(BFloat16))
exponent_max(::Type{BFloat16}) = Int(exponent_mask(BFloat16) >> significand_bits(BFloat16)) - exponent_bias(BFloat16)
exponent_raw_max(::Type{BFloat16}) = Int(exponent_mask(BFloat16) >> significand_bits(BFloat16))

include("conversions.jl")
include("functions.jl")
include("representation.jl")
include("util.jl")

end # module


text:
# Test Affine Constraints

using TrajOptSOCPs
using Test

@testset "Affine Equality Constraint Simple" begin
    # Create a simple AL_AffineEquality object
    dT = AL_AffineEquality([5], [10])

    # Check that it can determine when it is on the line
    @test TrajOptSOCPs.satisfied(dT, 2)

    # And likewise when it is not on the line
    @test !TrajOptSOCPs.satisfied(dT, 4)
    @test !TrajOptSOCPs.satisfied(dT, 0)

    # Run a simple evaluation
    @test TrajOptSOCPs.getRaw(dT, 2) == [0]
end

@testset "Affine Equality Constraint" begin
    # Create a higher dimensional AL_AffineEquality object
    dT2 = AL_AffineEquality([5 0; 0 6], [10; 12])

    # Check that it can determine when it is on the intersection
    @test TrajOptSOCPs.satisfied(dT2, [2; 2])

    # And likewise when it is not on the intersection
    ws = TrajOptSOCPs.whichSatisfied(dT2, [2; 4])
    #(should return (True, False))

    @test ws[1]
    @test !ws[2]

    # Test Projections and Violations
    testpt = [10; 12]
    # This point is not on the intersection:
    @test !TrajOptSOCPs.satisfied(dT2, testpt)

    # Get the projections and norm
    pv = TrajOptSOCPs.getProjVecs(dT2, testpt)
    gN = TrajOptSOCPs.getNormToProjVals(dT2, testpt)

    @test (gN'gN) != 0.0

    # Test the Calculus
    gradDT2 = TrajOptSOCPs.getGradC(dT2, testpt)
    hessDT2 = TrajOptSOCPs.getHessC(dT2)

    # Here we simply check the size.
    @test size(gradDT2, 1) == 2
    @test size(hessDT2) == (2, 2)

    # Test the the augmented lagrangian term also works
    alHess = TrajOptSOCPs.getHessC_ALTerm(dT2, testpt)

    @test size(alHess) == (2, 2)

end

@testset "Affine Inequality Constraint Simple" begin
    # Create a simple AL_AffineEquality object
    diT = AL_AffineInequality([5], [10])

    # Check that it can determine when it is on the line
    @test TrajOptSOCPs.satisfied(diT, 2)

    # And likewise when it is not on the line.
    # Above the line (false)
    @test !TrajOptSOCPs.satisfied(diT, 4)
    # Under the line (true)
    @test TrajOptSOCPs.satisfied(diT, 0)

    # Run a simple evaluation
    @test TrajOptSOCPs.getRaw(diT, 2) == [0]
end

@testset "Affine Inequality Constraint" begin
    # Create a higher dimensional AL_AffineInequality object
    diT2 = AL_AffineInequality([5 0; 0 6], [10; 12])

    # Check that it can determine when it is on the intersection
    @test TrajOptSOCPs.satisfied(diT2, [2; 2])

    # And likewise when it is not on the intersection
    ws = TrajOptSOCPs.whichSatisfied(diT2, [2; 4])
    #(should return (True, False))

    @test ws[1]
    @test !ws[2]

    # try a different set
    ws = TrajOptSOCPs.whichSatisfied(diT2, [2; 0])

    @test ws[1]
    @test ws[2]

    # Test Projections and Violations
    testpt = [10; 12]
    # This point is not on the intersection:
    @test !TrajOptSOCPs.satisfied(diT2, testpt)

    # Get the projections and violation
    pv = TrajOptSOCPs.getProjVecs(diT2, testpt)
    gN = TrajOptSOCPs.getNormToProjVals(diT2, testpt)

    @test (gN'gN) != 0.0

    # Try another test point
    testpt2 = [-10; -12]

    # This one should satisfy the constraint
    @test TrajOptSOCPs.satisfied(diT2, testpt2)

    # Get the projections and violation
    pv2 = TrajOptSOCPs.getProjVecs(diT2, testpt2)
    gN2 = TrajOptSOCPs.getNormToProjVals(diT2, testpt2)

    @test (gN2'gN2) == 0.0

    # Test the Calculus
    gradDT2 = TrajOptSOCPs.getGradC(diT2, testpt)
    gradDT22 = TrajOptSOCPs.getGradC(diT2, testpt2)
    hessDT2 = TrajOptSOCPs.getHessC(diT2)

    # Here we simply check the size.
    @test size(gradDT2, 1) == 2
    @test size(gradDT22, 1) == 2
    @test size(hessDT2) == (2, 2)

    # Test the the augmented lagrangian term also works
    alHess = TrajOptSOCPs.getHessC_ALTerm(diT2, testpt)

    @test size(alHess) == (2, 2)

end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
# Themelis, Stella, Patrinos, "Douglas-Rachford splitting and ADMM
# for nonconvex optimization: Accelerated and Newton-type linesearch algorithms",
# arXiv:2005.10230, 2020.
#
# https://arxiv.org/abs/2005.10230

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

function drls_default_gamma(f, muf, Lf, alpha, lambda)
    if muf !== nothing && muf > 0
        return 1 / (alpha * muf)
    end
    if ProximalOperators.is_convex(f)
        return alpha / Lf
    else
        return alpha * (2 - lambda) / (2 * Lf)
    end
end

function drls_C(f, muf, Lf, gamma, lambda)
    a = muf === nothing || muf <= 0 ? gamma * Lf : 1 / (gamma * muf)
    m = ProximalOperators.is_convex(f) ? max(a - lambda / 2, 0) : 1
    return (lambda / ((1 + a)^2) * ((2 - lambda) / 2 - a * m))
end

"""
    DRLSIteration(; <keyword-arguments>)

Instantiate the Douglas-Rachford line-search algorithm (see [1]) for solving optimization problems
of the form

    minimize f(x) + g(x),

where `f` is smooth.

# Arguments
- `x0`: initial point.
- `f=Zero()`: smooth objective term.
- `g=Zero()`: proximable objective term.
- `muf=nothing`: convexity modulus of f.
- `Lf=nothing`: Lipschitz constant of the gradient of f.
- `gamma`: stepsize to use, chosen appropriately based on Lf and muf by defaults.
- `max_backtracks=20`: maximum number of line-search backtracks.
- `directions=LBFGS(5)`: strategy to use to compute line-search directions.

# References
1. Themelis, Stella, Patrinos, "Douglas-Rachford splitting and ADMM for nonconvex optimization: Accelerated and Newton-type linesearch algorithms", arXiv:2005.10230, 2020.
"""
Base.@kwdef struct DRLSIteration{R,Tx,Tf,Tg,Tmuf,TLf,D}
    f::Tf = Zero()
    g::Tg = Zero()
    x0::Tx
    alpha::R = real(eltype(x0))(0.95)
    beta::R = real(eltype(x0))(0.5)
    lambda::R = real(eltype(x0))(1)
    muf::Tmuf = nothing
    Lf::TLf = nothing
    gamma::R = drls_default_gamma(f, muf, Lf, alpha, lambda)
    c::R = beta * drls_C(f, muf, Lf, gamma, lambda)
    dre_sign::Int = muf === nothing || muf <= 0 ? 1 : -1
    max_backtracks::Int = 20
    directions::D = LBFGS(5)
end

Base.IteratorSize(::Type{<:DRLSIteration}) = Base.IsInfinite()

Base.@kwdef mutable struct DRLSState{R,Tx,TH}
    x::Tx
    u::Tx
    v::Tx
    w::Tx
    res::Tx
    res_prev::Tx = similar(x)
    xbar::Tx
    xbar_prev::Tx = copy(xbar)
    d::Tx = similar(x)
    x_d::Tx = similar(x)
    gamma::R
    f_u::R
    g_v::R
    H::TH
    tau::R = zero(gamma)
    u0::Tx = similar(x)
    u1::Tx = similar(x)
    temp_x1::Tx = similar(x)
    temp_x2::Tx = similar(x)
end

DRE(f_u::Number, g_v::Number, x, u, res, gamma) = f_u + g_v - real(dot(x - u, res)) / gamma + 1 / (2 * gamma) * norm(res)^2

DRE(state::DRLSState) = DRE(state.f_u, state.g_v, state.x, state.u, state.res, state.gamma)

function Base.iterate(iter::DRLSIteration)
    x = copy(iter.x0)
    u, f_u = prox(iter.f, x, iter.gamma)
    w = 2 * u - x
    v, g_v = prox(iter.g, w, iter.gamma)
    res = u - v
    xbar = x - iter.lambda * res
    state = DRLSState(
        x=x, u=u, v=v, w=w, res=res, xbar=xbar, gamma=iter.gamma, f_u=f_u,
        g_v=g_v, H=initialize(iter.directions, x),
    )
    return state, state
end

set_next_direction!(::QuasiNewtonStyle, ::DRLSIteration, state::DRLSState) = mul!(state.d, state.H, -state.res)
set_next_direction!(::NesterovStyle, ::DRLSIteration, state::DRLSState) = state.d .= iterate(state.H)[1] .* (state.xbar .- state.xbar_prev) .+ (state.xbar .- state.x)
set_next_direction!(::NoAccelerationStyle, ::DRLSIteration, state::DRLSState) = state.d .= state.xbar .- state.x
set_next_direction!(iter::DRLSIteration, state::DRLSState) = set_next_direction!(acceleration_style(typeof(iter.directions)), iter, state)

update_direction_state!(::QuasiNewtonStyle, ::DRLSIteration, state::DRLSState) = update!(state.H, state.d, state.res - state.res_prev)
update_direction_state!(::NesterovStyle, ::DRLSIteration, state::DRLSState) = return
update_direction_state!(::NoAccelerationStyle, ::DRLSIteration, state::DRLSState) = return
update_direction_state!(iter::DRLSIteration, state::DRLSState) = update_direction_state!(acceleration_style(typeof(iter.directions)), iter, state)

function Base.iterate(iter::DRLSIteration{R}, state::DRLSState) where R
    DRE_curr = DRE(state)

    set_next_direction!(iter, state)

    state.x_d .= state.x .+ state.d
    state.xbar_prev, state.xbar = state.xbar, state.xbar_prev
    state.res_prev, state.res = state.res, state.res_prev
    state.tau = R(1)

    state.x .= state.x_d
    state.f_u = prox!(state.u, iter.f, state.x, iter.gamma)
    state.w .= 2 .* state.u .- state.x
    state.g_v = prox!(state.v, iter.g, state.w, iter.gamma)
    state.res .= state.u .- state.v
    state.xbar .= state.x .- iter.lambda * state.res

    update_direction_state!(iter, state)

    a, b, c = R(0), R(0), R(0)

    for k in 1:iter.max_backtracks
        if iter.dre_sign * DRE(state) <= iter.dre_sign * DRE_curr - iter.c / iter.gamma * norm(state.res_prev)^2
            break
        end

        state.tau = k == iter.max_backtracks ? R(0) : state.tau / 2
        state.x .= state.tau .* state.x_d .+ (1 - state.tau) .* state.xbar_prev

        if k == 1 && ProximalOperators.is_generalized_quadratic(iter.f)
            copyto!(state.u1, state.u)
            c = prox!(state.u0, iter.f, state.xbar_prev, iter.gamma)
            state.temp_x1 .= state.xbar_prev .- state.x_d
            state.temp_x2 .= state.xbar_prev .- state.u0
            b = real(dot(state.temp_x1, state.temp_x2)) / iter.gamma
            a = state.f_u - b - c
        end

        if ProximalOperators.is_generalized_quadratic(iter.f)
            state.u .= state.tau .* state.u1 .+ (1 - state.tau) .* state.u0
            state.f_u = a * state.tau ^ 2 + b * state.tau + c
        else
            state.f_u = prox!(state.u, iter.f, state.x, iter.gamma)
        end

        state.w .= 2 .* state.u .- state.x
        state.g_v = prox!(state.v, iter.g, state.w, iter.gamma)
        state.res .= state.u .- state.v
        state.xbar .= state.x .- iter.lambda * state.res
    end

    return state, state
end

# Solver

struct DRLS{R, K}
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int
    kwargs::K
end

function (solver::DRLS)(x0; kwargs...)
    stop(state::DRLSState) = norm(state.res, Inf) / state.gamma <= solver.tol
    disp((it, state)) = @printf(
        "%5d | %.3e | %.3e | %.3e\n",
        it,
        state.gamma / state.gamma,
        norm(state.res, Inf),
        state.tau,
    )
    iter = DRLSIteration(; x0=x0, solver.kwargs..., kwargs...)
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose
        iter = tee(sample(iter, solver.freq), disp)
    end
    num_iters, state_final = loop(iter)
    return state_final.u, state_final.v, num_iters
end

DRLS(; maxit=1_000, tol=1e-8, verbose=false, freq=10, kwargs...) =
    DRLS(maxit, tol, verbose, freq, kwargs)


text:
using Clustering
using Test
using CodecZlib

@testset "hclust() (hierarchical clustering)" begin

@testset "param checks" begin
    D = rand(5, 5)
    Dsym = D + D'
    Dnan = copy(Dsym)
    Dnan[1, 3] = Dnan[3, 1] = NaN
    @testset "hclust()" begin
        @test_throws ArgumentError hclust(Dsym, linkage=:typo)
        @test_throws ArgumentError hclust(D, linkage=:single)
        @test_throws ArgumentError hclust(Dnan, linkage=:single)
    end
    hclu = @inferred(hclust(Dsym, linkage=:single))
    @test hclu isa Clustering.Hclust
    @testset "cutree()" begin
        @test_throws ArgumentError cutree(hclu)
        @test_throws ArgumentError cutree(hclu, k=-1)
        @test_throws ArgumentError cutree(hclu, k=0)
        @test cutree(hclust(Dsym), k=10) isa Vector{Int}
        @test cutree(hclust(fill(0.0, 0, 0)), k=0) == Int[]
        @test cutree(hclust(fill(0.0, 0, 0)), k=1) == Int[]
    end
end

@testset "R hclust() generated examples" begin
# load the examples
hclu_examples_filename = joinpath(@__DIR__, "data", "hclust_generated_examples.jl.gz")
Base.include_string(@__MODULE__,
                    open(io -> read(io, String), GzipDecompressorStream, hclu_examples_filename),
                    hclu_examples_filename)

# test to make sure many random examples match R's implementation
@testset "example #$i (linkage=:$(example["linkage"]), n=$(size(example["D"], 1)))" for
        (i, example) in enumerate(examples)

    linkage = example["linkage"]
    hclu = @inferred(hclust(example["D"], linkage=linkage))
    @test hclu isa Clustering.Hclust
    @test Clustering.nnodes(hclu) == size(example["D"], 1)
    @test Clustering.nmerges(hclu) == Clustering.nnodes(hclu)-1
    @test Clustering.height(hclu) ≈ maximum(example["height"]) atol=1e-5
    @test hclu.merges == example["merge"]
    @test hclu.heights ≈ example["height"] atol=1e-5
    @test hclu.order == example["order"]

    # compare hclust_nn_lw() (the default) and hclust_nn() (slower) methods
    if linkage ∈ [:complete, :average]
        @testset "hclust_nn()" begin
            hclu2 = Hclust(Clustering.hclust_nn(example["D"], linkage == :complete ? Clustering.slicemaximum : Clustering.slicemean),
                           linkage)
            @test hclu2.merges == hclu.merges
            @test hclu2.heights ≈ hclu.heights atol=1e-5
            @test hclu2.order == hclu.order
        end
    end

    local cut_k = example["cut_k"]
    local cut_h = example["cut_h"]
    if cut_h !== nothing
        # due to small arithmetic differences between R and Julia heights might be slightly different
        # find the matching height
        cut_h_r = cut_h
        cut_h_ix = findmin(abs.(hclu.heights .- cut_h_r))[2]
        cut_h = hclu.heights[cut_h_ix]
        @assert isapprox(cut_h, cut_h_r, atol=1e-6) "h=$cut_h ≈ h_R=$cut_h_r"
    end
    @testset "cutree(hclu, k=$(repr(cut_k)), h=$(repr(cut_h)))" begin
        cutt = @inferred(cutree(hclu, k=cut_k, h=cut_h))
        @test cutt isa Vector{Int}
        @test length(cutt) == Clustering.nnodes(hclu)
        @test cutt == example["cutree"]
    end
end
end

@testset "hclust_n3()" begin
    # no thorough testing (it's O(N³)), just test one example
    example_n3 = examples[10]
    hclu_n3 = @inferred(Clustering.hclust_n3(example_n3["D"], Clustering.slicemaximum))
    @test hclu_n3.mleft == example_n3["merge"][:, 1]
    @test hclu_n3.mright == example_n3["merge"][:, 2]
    @test hclu_n3.heights ≈ example_n3["height"] atol=1e-5
end

local hclust_linkages = [:single, :average, :complete]

@testset "hclust(0×0 matrix, linkage=$linkage)" for linkage in hclust_linkages
    hclu = hclust(fill(0.0, 0, 0), linkage=linkage)
    @test Clustering.nnodes(hclu) == 0
    @test Clustering.nmerges(hclu) == 0
    @test Clustering.height(hclu) == -Inf
    cut1 = @inferred(cutree(hclu, h=1))
    @test cut1 == Int[]
end

@testset "hclust([$dist] 1×1 matrix, linkage=$linkage)" for
    linkage in hclust_linkages, dist in [-Inf, 0, 1, 2, Inf]
    hclu = hclust(fill(dist, 1, 1), linkage=linkage)
    @test Clustering.nnodes(hclu) == 1
    @test Clustering.nmerges(hclu) == 0
    @test Clustering.height(hclu) == -Inf
    cut1 = @inferred(cutree(hclu, h=1))
    @test cut1 == [1]
end

@testset "hclust(linkage=$linkage) when data contains an isolated point (#109)" for linkage in hclust_linkages
    # point #2 is isolated: distances to all the other points are Inf
    mdist = [
        0.0 Inf 0.1  Inf  Inf  Inf  Inf  Inf  Inf  Inf;
        Inf 0.0 Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf;
        0.1 Inf 0.0  0.11 Inf  Inf  Inf  Inf  Inf  Inf;
        Inf Inf 0.11 0.0  0.86 0.86 Inf  Inf  Inf  Inf;
        Inf Inf Inf  0.86 0.0  0.67 0.72 1.93 Inf  Inf;
        Inf Inf Inf  0.86 0.67 0.0  0.72 Inf  Inf  Inf;
        Inf Inf Inf  Inf  0.72 0.72 0.0  Inf  0.72 Inf;
        Inf Inf Inf  Inf  1.93 Inf  Inf  0.0  3.91 Inf;
        Inf Inf Inf  Inf  Inf  Inf  0.72 3.91 0.0  0.52;
        Inf Inf Inf  Inf  Inf  Inf  Inf  Inf  0.52 0.0]

    hclu = hclust(mdist, linkage=linkage)
    @test Clustering.nnodes(hclu) == 10
    @test Clustering.nmerges(hclu) == 9
    @test Clustering.height(hclu) == Inf
    cut1 = @inferred(cutree(hclu, h=1))
    @test cut1 isa Vector{Int}
    @test length(cut1) == Clustering.nnodes(hclu)
    if linkage == :single
        @test cut1 == [1, 2, 1, 1, 1, 1, 1, 3, 1, 1]
    end
end

@testset "cutree(hclust, h=h) when the height of all subtrees greater than h (#141)" begin
    A = [0.0 0.7; 0.7 0.0]
    hA = hclust(A, linkage=:average)
    @test cutree(hA, h=0.5) == [1, 2]
    @test cutree(hA, h=0.7) == [1, 1]
    @test cutree(hA, h=0.9) == [1, 1]
end

end


text:
# dependencies
using LFAToolkit
using LinearAlgebra
using Plots
# setup
p = 2
dimension = 1
mesh = Mesh1D(1.0)

# operator
diffusion = GalleryOperator("diffusion", p+1, p+1, mesh)

# NoPreco smoother
noPreco = NoPreco(diffusion)
# full operator symbols
numbersteps = 250
maxeigenvalue = 0
θ_min = -π/2
θ_max = 3π/2
θ_step = 2π/(numbersteps-1)
θ_range = θ_min:θ_step:θ_max

# compute and plot smoothing factor
# setup
ω = [1.00]
eigenvalues = zeros(numbersteps, p)

# compute
for i in 1:numbersteps
    θ = [θ_range[i]]
    if abs(θ[1]) >  π/512
        A = computesymbols(noPreco, ω, θ)
        currenteigenvalues = [real(val) for val in eigvals(I-A)]
        eigenvalues[i, :] = currenteigenvalues
    end
end

# plot
xrange = θ_range/π
plot(
    xrange,
    xlabel="θ/π",
    xtickfont=font(12, "Courier"),
    eigenvalues,
    ytickfont=font(12, "Courier"),
    ylabel="λ",
    linewidth=3,
    legend=:none,
    title="Spectrum of NoPreco Symbol",
    palette=palette(:tab10)
)
ylims!(min(0.0, eigenvalues...) * 1.1, max(eigenvalues...) * 1.1)
savefig("noPreco_spectrum_2")

## p = 4

# setup
p = 4
dimension = 1
mesh = Mesh1D(1.0)

# operator
diffusion = GalleryOperator("diffusion", p+1, p+1, mesh)

# NoPreco smoother
noPreco = NoPreco(diffusion)
# full operator symbols
numbersteps = 250
maxeigenvalue = 0
θ_min = -π/2
θ_max = 3π/2
θ_step = 2π/(numbersteps-1)
θ_range = θ_min:θ_step:θ_max

# compute and plot smoothing factor
# setup
ω = [1.00]
eigenvalues = zeros(numbersteps, p)

# compute
for i in 1:numbersteps
    θ = [θ_range[i]]
    if abs(θ[1]) >  π/512
        A = computesymbols(noPreco, ω, θ)
        currenteigenvalues = [real(val) for val in eigvals(I-A)]
        eigenvalues[i, :] = currenteigenvalues
    end
end

# plot
xrange = θ_range/π
plot(
    xrange,
    xlabel="θ/π",
    xtickfont=font(12, "Courier"),
    eigenvalues,
    ytickfont=font(12, "Courier"),
    ylabel="λ",
    linewidth=3,
    legend=:none,
    title="Spectrum of NoPreco Symbol",
    palette=palette(:tab10)
)
ylims!(min(0.0, eigenvalues...) * 1.1, max(eigenvalues...) * 1.1)
savefig("noPreco_spectrum_4")

## p = 7

# setup
p = 7
dimension = 1
mesh = Mesh1D(1.0)

# operator
diffusion = GalleryOperator("diffusion", p+1, p+1, mesh)

# NoPreco smoother
noPreco = NoPreco(diffusion)
# full operator symbols
numbersteps = 250
maxeigenvalue = 0
θ_min = -π/2
θ_max = 3π/2
θ_step = 2π/(numbersteps-1)
θ_range = θ_min:θ_step:θ_max

# compute and plot smoothing factor
# setup
ω = [1.00]
eigenvalues = zeros(numbersteps, p)

# compute
for i in 1:numbersteps
    θ = [θ_range[i]]
    if abs(θ[1]) >  π/512
        A = computesymbols(noPreco, ω, θ)
        currenteigenvalues = [real(val) for val in eigvals(I-A)]
        eigenvalues[i, :] = currenteigenvalues
    end
end

# plot
xrange = θ_range/π
plot(
    xrange,
    xlabel="θ/π",
    xtickfont=font(12, "Courier"),
    eigenvalues,
    ytickfont=font(12, "Courier"),
    ylabel="λ",
    linewidth=3,
    legend=:none,
    title="Spectrum of NoPreco Symbol",
    palette=palette(:tab10)
)
ylims!(min(0.0, eigenvalues...) * 1.1, max(eigenvalues...) * 1.1)
savefig("noPreco_spectrum_7")


text:
using Typeclass
import ..Types.Atom
import ..Abstract

@instance MolecularContainer Atom begin
    clone(x::Atom, new_props::Dict=Dict()) = Abstract.clone(x, new_props)
    obtain(x::Atom, k::Symbol) = Abstract.obtain(x, k)
    features(x::Atom) = Abstract.features(x)
    ==(x::Atom, y::Atom) = Abstract.==(x, y)
end


text:
const orig_stdin  = Ref{IO}()
const orig_stdout = Ref{IO}()
const orig_stderr = Ref{IO}()

# the following constants need to be initialized in init().
const read_stdout = Ref{Base.PipeEndpoint}()
const read_stderr = Ref{Base.PipeEndpoint}()


text:
"""
    downside_deviation(returns, threshold; method=:full)

Calculates the downside deviation / semi-standard deviation which captures the downside risk.

# Arguments
- `returns`:     Vector of asset returns.
- `threshold`:   Scalar value or vector denoting the threshold returns.
- `method`:      One of `:full` (default) or `:partial`. Indicates whether to use the number of all returns (`:full`), or only the number of returns below the threshold (`:partial`) in the denominator.
"""
function downside_deviation(returns, threshold; method::Symbol=:full)
    sqrt(lower_partial_moment(returns, threshold, 2, method))
end

"""
    upside_deviation(returns, threshold; method=:full)

Calculates the upside deviation / semi-standard deviation which captures the upside "risk".

# Arguments
- `returns`:     Vector of asset returns.
- `threshold`:   Scalar value or vector denoting the threshold returns.
- `method`:      One of `:full` (default) or `:partial`. Indicates whether to use the number of all returns (`:full`), or only the number of returns above the threshold (`:partial`) in the denominator.
"""
function upside_deviation(returns, threshold; method::Symbol=:full)
    sqrt(higher_partial_moment(returns, threshold, 2, method))
end


text:
# Julia wrapper for header: cexport.h
# Automatically generated using Clang.jl

function aiGetExportFormatCount()
    ccall((:aiGetExportFormatCount, assimp), Csize_t, ())
end

function aiGetExportFormatDescription(pIndex)
    ccall((:aiGetExportFormatDescription, assimp), Ptr{aiExportFormatDesc}, (Csize_t,), pIndex)
end

function aiReleaseExportFormatDescription(desc)
    ccall((:aiReleaseExportFormatDescription, assimp), Cvoid, (Ptr{aiExportFormatDesc},), desc)
end

function aiCopyScene(pIn, pOut)
    ccall((:aiCopyScene, assimp), Cvoid, (Ptr{aiScene}, Ptr{Ptr{aiScene}}), pIn, pOut)
end

function aiFreeScene(pIn)
    ccall((:aiFreeScene, assimp), Cvoid, (Ptr{aiScene},), pIn)
end

function aiExportScene(pScene, pFormatId, pFileName, pPreprocessing)
    ccall((:aiExportScene, assimp), aiReturn, (Ptr{aiScene}, Cstring, Cstring, UInt32), pScene, pFormatId, pFileName, pPreprocessing)
end

function aiExportSceneEx(pScene, pFormatId, pFileName, pIO, pPreprocessing)
    ccall((:aiExportSceneEx, assimp), aiReturn, (Ptr{aiScene}, Cstring, Cstring, Ptr{aiFileIO}, UInt32), pScene, pFormatId, pFileName, pIO, pPreprocessing)
end

function aiExportSceneToBlob(pScene, pFormatId, pPreprocessing)
    ccall((:aiExportSceneToBlob, assimp), Ptr{aiExportDataBlob}, (Ptr{aiScene}, Cstring, UInt32), pScene, pFormatId, pPreprocessing)
end

function aiReleaseExportBlob(pData)
    ccall((:aiReleaseExportBlob, assimp), Cvoid, (Ptr{aiExportDataBlob},), pData)
end
# Julia wrapper for header: cimport.h
# Automatically generated using Clang.jl

function aiImportFile(pFile, pFlags)
    ccall((:aiImportFile, assimp), Ptr{aiScene}, (Cstring, UInt32), pFile, pFlags)
end

function aiImportFileEx(pFile, pFlags, pFS)
    ccall((:aiImportFileEx, assimp), Ptr{aiScene}, (Cstring, UInt32, Ptr{aiFileIO}), pFile, pFlags, pFS)
end

function aiImportFileExWithProperties(pFile, pFlags, pFS, pProps)
    ccall((:aiImportFileExWithProperties, assimp), Ptr{aiScene}, (Cstring, UInt32, Ptr{aiFileIO}, Ptr{aiPropertyStore}), pFile, pFlags, pFS, pProps)
end

function aiImportFileFromMemory(pBuffer, pLength, pFlags, pHint)
    ccall((:aiImportFileFromMemory, assimp), Ptr{aiScene}, (Cstring, UInt32, UInt32, Cstring), pBuffer, pLength, pFlags, pHint)
end

function aiImportFileFromMemoryWithProperties(pBuffer, pLength, pFlags, pHint, pProps)
    ccall((:aiImportFileFromMemoryWithProperties, assimp), Ptr{aiScene}, (Cstring, UInt32, UInt32, Cstring, Ptr{aiPropertyStore}), pBuffer, pLength, pFlags, pHint, pProps)
end

function aiApplyPostProcessing(pScene, pFlags)
    ccall((:aiApplyPostProcessing, assimp), Ptr{aiScene}, (Ptr{aiScene}, UInt32), pScene, pFlags)
end

function aiGetPredefinedLogStream(pStreams, file)
    ccall((:aiGetPredefinedLogStream, assimp), aiLogStream, (aiDefaultLogStream, Cstring), pStreams, file)
end

function aiAttachLogStream(stream)
    ccall((:aiAttachLogStream, assimp), Cvoid, (Ptr{aiLogStream},), stream)
end

function aiEnableVerboseLogging(d)
    ccall((:aiEnableVerboseLogging, assimp), Cvoid, (aiBool,), d)
end

function aiDetachLogStream(stream)
    ccall((:aiDetachLogStream, assimp), aiReturn, (Ptr{aiLogStream},), stream)
end

function aiDetachAllLogStreams()
    ccall((:aiDetachAllLogStreams, assimp), Cvoid, ())
end

function aiReleaseImport(pScene)
    ccall((:aiReleaseImport, assimp), Cvoid, (Ptr{aiScene},), pScene)
end

function aiGetErrorString()
    ccall((:aiGetErrorString, assimp), Cstring, ())
end

function aiIsExtensionSupported(szExtension)
    ccall((:aiIsExtensionSupported, assimp), aiBool, (Cstring,), szExtension)
end

function aiGetExtensionList(szOut)
    ccall((:aiGetExtensionList, assimp), Cvoid, (Ptr{aiString},), szOut)
end

function aiGetMemoryRequirements(pIn, in)
    ccall((:aiGetMemoryRequirements, assimp), Cvoid, (Ptr{aiScene}, Ptr{aiMemoryInfo}), pIn, in)
end

function aiCreatePropertyStore()
    ccall((:aiCreatePropertyStore, assimp), Ptr{aiPropertyStore}, ())
end

function aiReleasePropertyStore(p)
    ccall((:aiReleasePropertyStore, assimp), Cvoid, (Ptr{aiPropertyStore},), p)
end

function aiSetImportPropertyInteger(store, szName, value)
    ccall((:aiSetImportPropertyInteger, assimp), Cvoid, (Ptr{aiPropertyStore}, Cstring, Cint), store, szName, value)
end

function aiSetImportPropertyFloat(store, szName, value)
    ccall((:aiSetImportPropertyFloat, assimp), Cvoid, (Ptr{aiPropertyStore}, Cstring, ai_real), store, szName, value)
end

function aiSetImportPropertyString(store, szName, st)
    ccall((:aiSetImportPropertyString, assimp), Cvoid, (Ptr{aiPropertyStore}, Cstring, Ptr{aiString}), store, szName, st)
end

function aiSetImportPropertyMatrix(store, szName, mat)
    ccall((:aiSetImportPropertyMatrix, assimp), Cvoid, (Ptr{aiPropertyStore}, Cstring, Ptr{aiMatrix4x4}), store, szName, mat)
end

function aiCreateQuaternionFromMatrix(quat, mat)
    ccall((:aiCreateQuaternionFromMatrix, assimp), Cvoid, (Ptr{aiQuaternion}, Ptr{aiMatrix3x3}), quat, mat)
end

function aiDecomposeMatrix(mat, scaling, rotation, position)
    ccall((:aiDecomposeMatrix, assimp), Cvoid, (Ptr{aiMatrix4x4}, Ptr{aiVector3D}, Ptr{aiQuaternion}, Ptr{aiVector3D}), mat, scaling, rotation, position)
end

function aiTransposeMatrix4(mat)
    ccall((:aiTransposeMatrix4, assimp), Cvoid, (Ptr{aiMatrix4x4},), mat)
end

function aiTransposeMatrix3(mat)
    ccall((:aiTransposeMatrix3, assimp), Cvoid, (Ptr{aiMatrix3x3},), mat)
end

function aiTransformVecByMatrix3(vec, mat)
    ccall((:aiTransformVecByMatrix3, assimp), Cvoid, (Ptr{aiVector3D}, Ptr{aiMatrix3x3}), vec, mat)
end

function aiTransformVecByMatrix4(vec, mat)
    ccall((:aiTransformVecByMatrix4, assimp), Cvoid, (Ptr{aiVector3D}, Ptr{aiMatrix4x4}), vec, mat)
end

function aiMultiplyMatrix4(dst, src)
    ccall((:aiMultiplyMatrix4, assimp), Cvoid, (Ptr{aiMatrix4x4}, Ptr{aiMatrix4x4}), dst, src)
end

function aiMultiplyMatrix3(dst, src)
    ccall((:aiMultiplyMatrix3, assimp), Cvoid, (Ptr{aiMatrix3x3}, Ptr{aiMatrix3x3}), dst, src)
end

function aiIdentityMatrix3(mat)
    ccall((:aiIdentityMatrix3, assimp), Cvoid, (Ptr{aiMatrix3x3},), mat)
end

function aiIdentityMatrix4(mat)
    ccall((:aiIdentityMatrix4, assimp), Cvoid, (Ptr{aiMatrix4x4},), mat)
end

function aiGetImportFormatCount()
    ccall((:aiGetImportFormatCount, assimp), Csize_t, ())
end

function aiGetImportFormatDescription(pIndex)
    ccall((:aiGetImportFormatDescription, assimp), Ptr{aiImporterDesc}, (Csize_t,), pIndex)
end

function aiVector2AreEqual(a, b)
    ccall((:aiVector2AreEqual, assimp), Cint, (Ptr{aiVector2D}, Ptr{aiVector2D}), a, b)
end

function aiVector2AreEqualEpsilon(a, b, epsilon)
    ccall((:aiVector2AreEqualEpsilon, assimp), Cint, (Ptr{aiVector2D}, Ptr{aiVector2D}, Cfloat), a, b, epsilon)
end

function aiVector2Add(dst, src)
    ccall((:aiVector2Add, assimp), Cvoid, (Ptr{aiVector2D}, Ptr{aiVector2D}), dst, src)
end

function aiVector2Subtract(dst, src)
    ccall((:aiVector2Subtract, assimp), Cvoid, (Ptr{aiVector2D}, Ptr{aiVector2D}), dst, src)
end

function aiVector2Scale(dst, s)
    ccall((:aiVector2Scale, assimp), Cvoid, (Ptr{aiVector2D}, Cfloat), dst, s)
end

function aiVector2SymMul(dst, other)
    ccall((:aiVector2SymMul, assimp), Cvoid, (Ptr{aiVector2D}, Ptr{aiVector2D}), dst, other)
end

function aiVector2DivideByScalar(dst, s)
    ccall((:aiVector2DivideByScalar, assimp), Cvoid, (Ptr{aiVector2D}, Cfloat), dst, s)
end

function aiVector2DivideByVector(dst, v)
    ccall((:aiVector2DivideByVector, assimp), Cvoid, (Ptr{aiVector2D}, Ptr{aiVector2D}), dst, v)
end

function aiVector2Length(v)
    ccall((:aiVector2Length, assimp), Cfloat, (Ptr{aiVector2D},), v)
end

function aiVector2SquareLength(v)
    ccall((:aiVector2SquareLength, assimp), Cfloat, (Ptr{aiVector2D},), v)
end

function aiVector2Negate(dst)
    ccall((:aiVector2Negate, assimp), Cvoid, (Ptr{aiVector2D},), dst)
end

function aiVector2DotProduct(a, b)
    ccall((:aiVector2DotProduct, assimp), Cfloat, (Ptr{aiVector2D}, Ptr{aiVector2D}), a, b)
end

function aiVector2Normalize(v)
    ccall((:aiVector2Normalize, assimp), Cvoid, (Ptr{aiVector2D},), v)
end

function aiVector3AreEqual(a, b)
    ccall((:aiVector3AreEqual, assimp), Cint, (Ptr{aiVector3D}, Ptr{aiVector3D}), a, b)
end

function aiVector3AreEqualEpsilon(a, b, epsilon)
    ccall((:aiVector3AreEqualEpsilon, assimp), Cint, (Ptr{aiVector3D}, Ptr{aiVector3D}, Cfloat), a, b, epsilon)
end

function aiVector3LessThan(a, b)
    ccall((:aiVector3LessThan, assimp), Cint, (Ptr{aiVector3D}, Ptr{aiVector3D}), a, b)
end

function aiVector3Add(dst, src)
    ccall((:aiVector3Add, assimp), Cvoid, (Ptr{aiVector3D}, Ptr{aiVector3D}), dst, src)
end

function aiVector3Subtract(dst, src)
    ccall((:aiVector3Subtract, assimp), Cvoid, (Ptr{aiVector3D}, Ptr{aiVector3D}), dst, src)
end

function aiVector3Scale(dst, s)
    ccall((:aiVector3Scale, assimp), Cvoid, (Ptr{aiVector3D}, Cfloat), dst, s)
end

function aiVector3SymMul(dst, other)
    ccall((:aiVector3SymMul, assimp), Cvoid, (Ptr{aiVector3D}, Ptr{aiVector3D}), dst, other)
end

function aiVector3DivideByScalar(dst, s)
    ccall((:aiVector3DivideByScalar, assimp), Cvoid, (Ptr{aiVector3D}, Cfloat), dst, s)
end

function aiVector3DivideByVector(dst, v)
    ccall((:aiVector3DivideByVector, assimp), Cvoid, (Ptr{aiVector3D}, Ptr{aiVector3D}), dst, v)
end

function aiVector3Length(v)
    ccall((:aiVector3Length, assimp), Cfloat, (Ptr{aiVector3D},), v)
end

function aiVector3SquareLength(v)
    ccall((:aiVector3SquareLength, assimp), Cfloat, (Ptr{aiVector3D},), v)
end

function aiVector3Negate(dst)
    ccall((:aiVector3Negate, assimp), Cvoid, (Ptr{aiVector3D},), dst)
end

function aiVector3DotProduct(a, b)
    ccall((:aiVector3DotProduct, assimp), Cfloat, (Ptr{aiVector3D}, Ptr{aiVector3D}), a, b)
end

function aiVector3CrossProduct(dst, a, b)
    ccall((:aiVector3CrossProduct, assimp), Cvoid, (Ptr{aiVector3D}, Ptr{aiVector3D}, Ptr{aiVector3D}), dst, a, b)
end

function aiVector3Normalize(v)
    ccall((:aiVector3Normalize, assimp), Cvoid, (Ptr{aiVector3D},), v)
end

function aiVector3NormalizeSafe(v)
    ccall((:aiVector3NormalizeSafe, assimp), Cvoid, (Ptr{aiVector3D},), v)
end

function aiVector3RotateByQuaternion(v, q)
    ccall((:aiVector3RotateByQuaternion, assimp), Cvoid, (Ptr{aiVector3D}, Ptr{aiQuaternion}), v, q)
end

function aiMatrix3FromMatrix4(dst, mat)
    ccall((:aiMatrix3FromMatrix4, assimp), Cvoid, (Ptr{aiMatrix3x3}, Ptr{aiMatrix4x4}), dst, mat)
end

function aiMatrix3FromQuaternion(mat, q)
    ccall((:aiMatrix3FromQuaternion, assimp), Cvoid, (Ptr{aiMatrix3x3}, Ptr{aiQuaternion}), mat, q)
end

function aiMatrix3AreEqual(a, b)
    ccall((:aiMatrix3AreEqual, assimp), Cint, (Ptr{aiMatrix3x3}, Ptr{aiMatrix3x3}), a, b)
end

function aiMatrix3AreEqualEpsilon(a, b, epsilon)
    ccall((:aiMatrix3AreEqualEpsilon, assimp), Cint, (Ptr{aiMatrix3x3}, Ptr{aiMatrix3x3}, Cfloat), a, b, epsilon)
end

function aiMatrix3Inverse(mat)
    ccall((:aiMatrix3Inverse, assimp), Cvoid, (Ptr{aiMatrix3x3},), mat)
end

function aiMatrix3Determinant(mat)
    ccall((:aiMatrix3Determinant, assimp), Cfloat, (Ptr{aiMatrix3x3},), mat)
end

function aiMatrix3RotationZ(mat, angle)
    ccall((:aiMatrix3RotationZ, assimp), Cvoid, (Ptr{aiMatrix3x3}, Cfloat), mat, angle)
end

function aiMatrix3FromRotationAroundAxis(mat, axis, angle)
    ccall((:aiMatrix3FromRotationAroundAxis, assimp), Cvoid, (Ptr{aiMatrix3x3}, Ptr{aiVector3D}, Cfloat), mat, axis, angle)
end

function aiMatrix3Translation(mat, translation)
    ccall((:aiMatrix3Translation, assimp), Cvoid, (Ptr{aiMatrix3x3}, Ptr{aiVector2D}), mat, translation)
end

function aiMatrix3FromTo(mat, from, to)
    ccall((:aiMatrix3FromTo, assimp), Cvoid, (Ptr{aiMatrix3x3}, Ptr{aiVector3D}, Ptr{aiVector3D}), mat, from, to)
end

function aiMatrix4FromMatrix3(dst, mat)
    ccall((:aiMatrix4FromMatrix3, assimp), Cvoid, (Ptr{aiMatrix4x4}, Ptr{aiMatrix3x3}), dst, mat)
end

function aiMatrix4FromScalingQuaternionPosition(mat, scaling, rotation, position)
    ccall((:aiMatrix4FromScalingQuaternionPosition, assimp), Cvoid, (Ptr{aiMatrix4x4}, Ptr{aiVector3D}, Ptr{aiQuaternion}, Ptr{aiVector3D}), mat, scaling, rotation, position)
end

function aiMatrix4Add(dst, src)
    ccall((:aiMatrix4Add, assimp), Cvoid, (Ptr{aiMatrix4x4}, Ptr{aiMatrix4x4}), dst, src)
end

function aiMatrix4AreEqual(a, b)
    ccall((:aiMatrix4AreEqual, assimp), Cint, (Ptr{aiMatrix4x4}, Ptr{aiMatrix4x4}), a, b)
end

function aiMatrix4AreEqualEpsilon(a, b, epsilon)
    ccall((:aiMatrix4AreEqualEpsilon, assimp), Cint, (Ptr{aiMatrix4x4}, Ptr{aiMatrix4x4}, Cfloat), a, b, epsilon)
end

function aiMatrix4Inverse(mat)
    ccall((:aiMatrix4Inverse, assimp), Cvoid, (Ptr{aiMatrix4x4},), mat)
end

function aiMatrix4Determinant(mat)
    ccall((:aiMatrix4Determinant, assimp), Cfloat, (Ptr{aiMatrix4x4},), mat)
end

function aiMatrix4IsIdentity(mat)
    ccall((:aiMatrix4IsIdentity, assimp), Cint, (Ptr{aiMatrix4x4},), mat)
end

function aiMatrix4DecomposeIntoScalingEulerAnglesPosition(mat, scaling, rotation, position)
    ccall((:aiMatrix4DecomposeIntoScalingEulerAnglesPosition, assimp), Cvoid, (Ptr{aiMatrix4x4}, Ptr{aiVector3D}, Ptr{aiVector3D}, Ptr{aiVector3D}), mat, scaling, rotation, position)
end

function aiMatrix4DecomposeIntoScalingAxisAnglePosition(mat, scaling, axis, angle, position)
    ccall((:aiMatrix4DecomposeIntoScalingAxisAnglePosition, assimp), Cvoid, (Ptr{aiMatrix4x4}, Ptr{aiVector3D}, Ptr{aiVector3D}, Ptr{ai_real}, Ptr{aiVector3D}), mat, scaling, axis, angle, position)
end

function aiMatrix4DecomposeNoScaling(mat, rotation, position)
    ccall((:aiMatrix4DecomposeNoScaling, assimp), Cvoid, (Ptr{aiMatrix4x4}, Ptr{aiQuaternion}, Ptr{aiVector3D}), mat, rotation, position)
end

function aiMatrix4FromEulerAngles(mat, x, y, z)
    ccall((:aiMatrix4FromEulerAngles, assimp), Cvoid, (Ptr{aiMatrix4x4}, Cfloat, Cfloat, Cfloat), mat, x, y, z)
end

function aiMatrix4RotationX(mat, angle)
    ccall((:aiMatrix4RotationX, assimp), Cvoid, (Ptr{aiMatrix4x4}, Cfloat), mat, angle)
end

function aiMatrix4RotationY(mat, angle)
    ccall((:aiMatrix4RotationY, assimp), Cvoid, (Ptr{aiMatrix4x4}, Cfloat), mat, angle)
end

function aiMatrix4RotationZ(mat, angle)
    ccall((:aiMatrix4RotationZ, assimp), Cvoid, (Ptr{aiMatrix4x4}, Cfloat), mat, angle)
end

function aiMatrix4FromRotationAroundAxis(mat, axis, angle)
    ccall((:aiMatrix4FromRotationAroundAxis, assimp), Cvoid, (Ptr{aiMatrix4x4}, Ptr{aiVector3D}, Cfloat), mat, axis, angle)
end

function aiMatrix4Translation(mat, translation)
    ccall((:aiMatrix4Translation, assimp), Cvoid, (Ptr{aiMatrix4x4}, Ptr{aiVector3D}), mat, translation)
end

function aiMatrix4Scaling(mat, scaling)
    ccall((:aiMatrix4Scaling, assimp), Cvoid, (Ptr{aiMatrix4x4}, Ptr{aiVector3D}), mat, scaling)
end

function aiMatrix4FromTo(mat, from, to)
    ccall((:aiMatrix4FromTo, assimp), Cvoid, (Ptr{aiMatrix4x4}, Ptr{aiVector3D}, Ptr{aiVector3D}), mat, from, to)
end

function aiQuaternionFromEulerAngles(q, x, y, z)
    ccall((:aiQuaternionFromEulerAngles, assimp), Cvoid, (Ptr{aiQuaternion}, Cfloat, Cfloat, Cfloat), q, x, y, z)
end

function aiQuaternionFromAxisAngle(q, axis, angle)
    ccall((:aiQuaternionFromAxisAngle, assimp), Cvoid, (Ptr{aiQuaternion}, Ptr{aiVector3D}, Cfloat), q, axis, angle)
end

function aiQuaternionFromNormalizedQuaternion(q, normalized)
    ccall((:aiQuaternionFromNormalizedQuaternion, assimp), Cvoid, (Ptr{aiQuaternion}, Ptr{aiVector3D}), q, normalized)
end

function aiQuaternionAreEqual(a, b)
    ccall((:aiQuaternionAreEqual, assimp), Cint, (Ptr{aiQuaternion}, Ptr{aiQuaternion}), a, b)
end

function aiQuaternionAreEqualEpsilon(a, b, epsilon)
    ccall((:aiQuaternionAreEqualEpsilon, assimp), Cint, (Ptr{aiQuaternion}, Ptr{aiQuaternion}, Cfloat), a, b, epsilon)
end

function aiQuaternionNormalize(q)
    ccall((:aiQuaternionNormalize, assimp), Cvoid, (Ptr{aiQuaternion},), q)
end

function aiQuaternionConjugate(q)
    ccall((:aiQuaternionConjugate, assimp), Cvoid, (Ptr{aiQuaternion},), q)
end

function aiQuaternionMultiply(dst, q)
    ccall((:aiQuaternionMultiply, assimp), Cvoid, (Ptr{aiQuaternion}, Ptr{aiQuaternion}), dst, q)
end

function aiQuaternionInterpolate(dst, start, _end, factor)
    ccall((:aiQuaternionInterpolate, assimp), Cvoid, (Ptr{aiQuaternion}, Ptr{aiQuaternion}, Ptr{aiQuaternion}, Cfloat), dst, start, _end, factor)
end
# Julia wrapper for header: material.h
# Automatically generated using Clang.jl

function TextureTypeToString(in)
    ccall((:TextureTypeToString, assimp), Cstring, (aiTextureType,), in)
end

function aiGetMaterialProperty(pMat, pKey, type, index, pPropOut)
    ccall((:aiGetMaterialProperty, assimp), aiReturn, (Ptr{aiMaterial}, Cstring, UInt32, UInt32, Ptr{Ptr{aiMaterialProperty}}), pMat, pKey, type, index, pPropOut)
end

function aiGetMaterialFloatArray(pMat, pKey, type, index, pOut, pMax)
    ccall((:aiGetMaterialFloatArray, assimp), aiReturn, (Ptr{aiMaterial}, Cstring, UInt32, UInt32, Ptr{ai_real}, Ptr{UInt32}), pMat, pKey, type, index, pOut, pMax)
end

function aiGetMaterialIntegerArray(pMat, pKey, type, index, pOut, pMax)
    ccall((:aiGetMaterialIntegerArray, assimp), aiReturn, (Ptr{aiMaterial}, Cstring, UInt32, UInt32, Ptr{Cint}, Ptr{UInt32}), pMat, pKey, type, index, pOut, pMax)
end

function aiGetMaterialColor(pMat, pKey, type, index, pOut)
    ccall((:aiGetMaterialColor, assimp), aiReturn, (Ptr{aiMaterial}, Cstring, UInt32, UInt32, Ptr{aiColor4D}), pMat, pKey, type, index, pOut)
end

function aiGetMaterialUVTransform(pMat, pKey, type, index, pOut)
    ccall((:aiGetMaterialUVTransform, assimp), aiReturn, (Ptr{aiMaterial}, Cstring, UInt32, UInt32, Ptr{aiUVTransform}), pMat, pKey, type, index, pOut)
end

function aiGetMaterialString(pMat, pKey, type, index, pOut)
    ccall((:aiGetMaterialString, assimp), aiReturn, (Ptr{aiMaterial}, Cstring, UInt32, UInt32, Ptr{aiString}), pMat, pKey, type, index, pOut)
end

function aiGetMaterialTextureCount(pMat, type)
    ccall((:aiGetMaterialTextureCount, assimp), UInt32, (Ptr{aiMaterial}, aiTextureType), pMat, type)
end

function aiGetMaterialTexture(mat, type, index, path, mapping, uvindex, blend, op, mapmode, flags)
    ccall((:aiGetMaterialTexture, assimp), aiReturn, (Ptr{aiMaterial}, aiTextureType, UInt32, Ptr{aiString}, Ptr{aiTextureMapping}, Ptr{UInt32}, Ptr{ai_real}, Ptr{aiTextureOp}, Ptr{aiTextureMapMode}, Ptr{UInt32}), mat, type, index, path, mapping, uvindex, blend, op, mapmode, flags)
end
# Julia wrapper for header: qnan.h
# Automatically generated using Clang.jl

function is_qnan()
    ccall((:is_qnan, assimp), Cint, ())
end

function is_special_float()
    ccall((:is_special_float, assimp), Cint, ())
end

function get_qnan()
    ccall((:get_qnan, assimp), ai_real, ())
end
# Julia wrapper for header: version.h
# Automatically generated using Clang.jl

function aiGetLegalString()
    ccall((:aiGetLegalString, assimp), Cstring, ())
end

function aiGetVersionPatch()
    ccall((:aiGetVersionPatch, assimp), UInt32, ())
end

function aiGetVersionMinor()
    ccall((:aiGetVersionMinor, assimp), UInt32, ())
end

function aiGetVersionMajor()
    ccall((:aiGetVersionMajor, assimp), UInt32, ())
end

function aiGetVersionRevision()
    ccall((:aiGetVersionRevision, assimp), UInt32, ())
end

function aiGetBranchName()
    ccall((:aiGetBranchName, assimp), Cstring, ())
end

function aiGetCompileFlags()
    ccall((:aiGetCompileFlags, assimp), UInt32, ())
end


text:
# CIM.jl: different algorithms to simulate coherent Ising machines.

export
       OpticalOscillators,
       OPODynamics,
       evolve_optical_oscillators,
       noisy_mean_field_annealing

@inline zeros_like(x::AbstractArray) = zeros(eltype(x), size(x))

# Optical Parametric Oscillator (OPO)
# Based on https://arxiv.org/pdf/1901.08927.pdf
"""
$(TYPEDEF)

$(TYPEDFIELDS)

"""
struct OpticalOscillators{T <: Real}
    ig::IsingGraph
    scale::T
    noise::Distributions.Sampleable
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)

"""
struct OPODynamics{T <: Real}
    initial_state::Vector{T}
    saturation::T
    pump::Vector{T}
    momentum::T
end

"""
$(TYPEDSIGNATURES)

"""
function evolve_optical_oscillators(
    opo::OpticalOscillators{T},
    dyn::OPODynamics{T}
)  where T <: Real
    J, h = couplings(opo.ig), biases(opo.ig)
    x = dyn.initial_state
    L = length(x)
    Δm = zeros_like(x)
    for p ∈ dyn.pump
        Δx = p .* x .+ opo.scale .* (J * x .+ h) .+ rand(opo.noise, L)
        m = (1.0 - dyn.momentum) .* Δx + dyn.momentum .* Δm
        x .+= m .* (abs.(x .+ m) .< dyn.saturation)
        Δm = m
    end
    Int.(sign.(x))
end

# Noisy mean-field annealing (NMFA)
# Based on https://arxiv.org/pdf/1806.08422.pdf
"""
$(TYPEDSIGNATURES)

"""
function noisy_mean_field_annealing(
    opo::OpticalOscillators{T},
    dyn::OPODynamics{T}
)  where T <: Real
    J, h = couplings(opo.ig), biases(opo.ig)
    x = dyn.initial_state
    L = length(x)
    nmr = sqrt.(h .^ 2 .+ dropdims(sum(J .^ 2, dims=2), dims=2))
    for p ∈ dyn.pump
        ϕ = (J * x .+ h) ./ nmr .+ rand(opo.noise, L)
        x = (1.0 - dyn.momentum) .* x .- dyn.momentum .* tanh.(ϕ ./ p)
    end
    Int.(sign.(x))
end

# Degenerate Optical Parametric Oscillator (DOPO)
# Based on https://www.nature.com/articles/s41467-018-07328-1
"""
$(TYPEDEF)

$(TYPEDFIELDS)

"""
struct DegenerateOscillators{T <: Real}
    ig::IsingGraph
    pump::T
    saturation::T
end

"""
$(TYPEDSIGNATURES)

"""
function evolve_degenerate_oscillators()
    # To be written
end