text:
#For each dataset (mock communities and mvco timeseries):
#import tsv taxbyseq file, create unique sequence identifier, create fasta file for blast
#run blast from shell
#import results, screen and parse, write out fasta file of Syn sequences for next steps

#Kristen R. Hunter-Cevera, Marine Biological Laboratory, Spring 2020

using DataFrames, CSV
using BioSequences
using Gadfly
using JLD2

##
mock_df = DataFrame(CSV.File("/home/kristen/Documents/V6V8_analysis/data/vamps_downloads/taxbyseq_tables/mock_community_taxbyseq.tsv"; delim = '\t',header = 2))

#for now, remove dock samples:
delete!(mock_df,Symbol("syn_mock_communities--DockA_run7"))
select!(mock_df, Not( Symbol("syn_mock_communities--DockB_run7")))

rename!(mock_df,Symbol("syn_mock_communities--SynCommunityA_run7") => :SynA)
rename!(mock_df,Symbol("syn_mock_communities--SynCommunityB_run7") => :SynB)
rename!(mock_df,Symbol("syn_mock_communities--SynCommunityC_run7") => :SynC)
rename!(mock_df,Symbol("syn_mock_communities--SynCommunityD_run7") => :SynD)
rename!(mock_df,Symbol("syn_mock_communities--M1_run3")  => :M1)
rename!(mock_df,Symbol("syn_mock_communities--MVM2_run2") => :M2)

#remove rows that have totals of zeros after dock samples removed:
[sum(col) for col = eachcol(mock_df[!,[:SynA, :SynB]])]
zeroind=findall([sum(row) for row = eachrow(mock_df[!,Between(:M1,:SynD)])] .== 0)
deleterows!(mock_df, zeroind)
synind = findall(map(x -> occursin(r"Synecho",x), mock_df[:,:Taxonomy]))
#colwise(sum,mock_df[[ Symbol("syn_mock_communities--SynCommunityA_run7"),  Symbol("syn_mock_communities--SynCommunityB_run7")]])
mock_df.uniq_ID=collect(1:size(mock_df,1))

## write FASTA file of unqiue sequences with unique identifiers:
fastafile="/home/kristen/Documents/V6V8_analysis/analysis_products/blast/mock_communities.fasta"
w = open(FASTA.Writer, fastafile)
for q=1:size(mock_df,1)
    rec = FASTA.Record(string(mock_df[q,:uniq_ID]),mock_df[q,:Sequence])
    write(w, rec)
end
close(w)

## run blast:

outfile="/home/kristen/Documents/V6V8_analysis/analysis_products/blast/mock_community_blast_results"
run(`bash /home/kristen/Documents/V6V8_analysis/scripts/taxa_processing_scripts/run_blast.sh $fastafile $outfile`)

## Import the blast output file and merge with database:

#convert output file from blastn search to DataFrame
blast_df = DataFrame(CSV.File("/home/kristen/Documents/V6V8_analysis/analysis_products/blast/mock_community_blast_results"; delim = '\t', header = false))

#put headers on output file
blast_col = ["uniq_ID", "sseqid" ,"pident" ,"length", "mismatch","gapopen", "qstart", "qend", "sstart", "send", "evalue", "bitscore"];
names!(blast_df, Symbol.(blast_col))

#merge dataframes:
test_df=join(mock_df,blast_df,on = :uniq_ID, kind=:left)
test_df.query_length=abs.(test_df[!,:qstart] .- test_df[!,:qend])
test_df.subject_length=abs.(test_df[!,:sstart] .- test_df[!,:send])

synind=findall(occursin.(r"Synecho",test_df[!,:Taxonomy]))
noblast_hitsyn = findall(ismissing.(test_df[synind,:pident]))
#good that this is empty - all gast syn had blast syn hits
ii=setdiff(1:length(synind),noblast_hitsyn)
synind=synind[ii]

## some logic tests to see how pident, bitscore, etc relate....plus some plots

# 2 QUESTIONS:
#Of the sequences classified as Syn by GAST - how many are not?
#Of the sequences not classified as Syn by GAST - how many might be?

#🌵 Question 1:

#histograms to show distribution of blast results:
Gadfly.set_default_plot_size(20cm,20cm)
p1=plot(x=test_df[synind,:pident], Geom.histogram,Guide.title("Percent ID")) #most should be > 95
p2=plot(x=test_df[synind,:mismatch], Geom.histogram,Guide.title("Mismatch")) #most should be < 10
p3=plot(x=test_df[synind,:gapopen], Geom.histogram,Guide.title("Gaps Open")) #most should be < 3
p4=plot(x=test_df[synind,:bitscore], Geom.histogram,Guide.title("BitScore")) #should be above 800

gridstack([p1 p2; p3 p4])

##
tt=findall(test_df[synind,:pident] .<= 90) #so not too many...and usually low abundance
[sum(col) for col = eachcol(mock_df[synind[tt],Between(:M1,:SynD)])]

## if needed, write these out to blast against full nucleotide database:
w = open(FASTA.Writer, "/home/kristen/Documents/V6V8_analysis/mock_syngast_lowPID.fasta")
for q=1:size(tt,1)
    rec = FASTA.Record(string(mock_df[synind[tt[q]],:uniq_ID]),mock_df[synind[tt[q]],:Sequence])
    write(w, rec)
end
close(w)

# so, all in all for the mock syn - no worries for question 1! ✅

#🌵 Question 2:
noblast = findall(ismissing.(test_df[!,:pident])) #no blast results to syn
synblast = setdiff(1:size(test_df,1),noblast) #total blast results that had match to syn
nonsyn = setdiff(synblast,synind) #remove gast matches to syn, what came back as syn?

plot(x=test_df[nonsyn,:pident], Geom.histogram) #most should be > 95
#most of these sequences are below 90...which is good, but what about the ones that arent't?

#look at the sequences that really show a good chance of being a coastal syn:
#tt=findall((test_df[nonsyn,:pident] .>= 90) .& (test_df[nonsyn,:query_length] .>= 440) .& (test_df[nonsyn,:subject_length] .>= 440)
#    .& (test_df[nonsyn,:gapopen] .<= 3) .& (test_df[nonsyn,:mismatch] .<= 10) .& (test_df[nonsyn,:bitscore] .>= 800)) #only consider full length queries

#-or-:

#because these are mock communities - what are the other cyanos?
ss=findall(map(x -> occursin(r"Cyano",x),test_df[nonsyn,:Taxonomy]))
[sum(col) for col = eachcol(test_df[nonsyn[ss],Between(:M1,:SynD)])] #still very few

tt=findall(test_df[nonsyn,:bitscore] .>= 700) #only consider full length queries
plot(x=test_df[nonsyn[ss],:mismatch], Geom.histogram)

taxa_to_check=unique(test_df[nonsyn[tt],:Taxonomy])

[sum(col) for col = eachcol(test_df[nonsyn[tt],Between(:M1,:SynD)])] #still very few

rowsum=[sum(row) for row = eachrow(test_df[nonsyn[tt],Between(:M1,:SynD)])] #still very few
sum(rowsum) #this is also a low amount...

# ii=findall(map(x -> occursin(r"Prochlorococcus",x),test_df[nonsyn[tt],:Taxonomy]))

## hmmm...so some cyanos that are matching to Pro, Cyanobium, etc. Parse these for blast nucleotide check!
wcyano = open(FASTA.Writer, "/home/kristen/Documents/V6V8_analysis/mock_nonsyngast_highPID_cyano.fasta")
wother = open(FASTA.Writer, "/home/kristen/Documents/V6V8_analysis/mock_nonsyngast_highPID_other.fasta")

for q=1:size(tt,1)
    if occursin(r"Cyano",mock_df[nonsyn[tt[q]],:Taxonomy])
        rec = FASTA.Record(string(mock_df[nonsyn[tt[q]],:uniq_ID]),mock_df[nonsyn[tt[q]],:Sequence])
        write(wcyano, rec)
    else
        rec = FASTA.Record(string(mock_df[nonsyn[tt[q]],:uniq_ID]),mock_df[nonsyn[tt[q]],:Sequence])
        write(wother, rec)
    end
end
close(wcyano)
close(wother)
## hmmmm...so it does seem that these sequences could be considered Syn..despite being labeled as Pro or a Microcytis...probably should include them!
# Write out sequences in oliogytping format!

#index to use:
sind=union(synind,nonsyn[tt])
df=test_df[sind,:]
# Compile for oligotyping! woot-woot!

#wA = open(FASTA.Writer, "/home/kristen/Documents/V6V8_analysis/analysis_products/syn_sequences/mock_synseqs_2016exp.fasta")
w = open(FASTA.Writer, "/home/kristen/Documents/V6V8_analysis/analysis_products/syn_sequences/mock_synseqs_2018exp.fasta")

##

for n in names(df[!,Between(:SynA,:SynD)]) #names(df[!,Between(:M1,:M2)])#names(df[!,Between(:SynA,:SynD)]) #names(df[!,Between(:MV308_rd2,:MV400Ster)])  #names(df[!,Between(:MV263_rd2,:MV375)])

        ind=findall(map(x -> x != 0, df[:,n])) #returns the indexes that are not zero
        #df[ind,[n, 68]]
        @show n
        global test=sum(df[ind,n])
        m=match(r"(?<sample>(M[1,2]{1}|Syn[A,B,C,D]{1}))",string(n))

        #m=match(r"(?<sample>MV\d{3}[A,B]{0,1}\w*)",string(n)) #for MVCO time series
        #@show m[:smaple]

        global count=0 #important for oligotyping file - adds number to each sequence
        for j=1:length(ind) #different Syn seqs in this sample
            for i=1:df[ind[j],n] #for each syn seq, duplicate to create fasta file for oligotyping
                #@show
                global count += 1
                rec = FASTA.Record(string("Sample-",m[:sample],"_Read",count),df[ind[j],:Sequence])
                write(w, rec)
            end
        end
        @show test
        @show count
        if count != test
            print("Nooooo - sum does not match sequences printed?")
        end
    #end
end

close(w)

max_stars_repo_name: hsosik/NES-LTER

max_issues_repo_name: hsosik/NES-LTER

max_forks_repo_name: hsosik/NES-LTER


text:
abstract type AbstractBranchForm end

abstract type AbstractLineForm <: AbstractBranchForm end

abstract type AbstractDCLineForm <: AbstractBranchForm end

abstract type AbstractTransformerForm <: AbstractBranchForm end

abstract type PiLine <: AbstractLineForm end

abstract type SeriesLine <: AbstractLineForm end

abstract type SimpleHVDC <: AbstractDCLineForm end

include("branches/flow_variables.jl")
include("branches/bus_variables.jl")
include("branches/flow_constraints.jl")
include("branches/ptdf_model.jl")

ext: jl


text:
#
# Phase lag index
#
# See Stam, C. J., Nolte, G., & Daffertshofer, A. (2007).
# Phase lag index: Assessment of functional connectivity from multi
# channel EEG and MEG with diminished bias from common sources.
# Human Brain Mapping, 28(11), 1178–1193. doi:10.1002/hbm.20346
struct PLI <: PairwiseStatistic; end

finish(::Type{PLI}, x::Complex, n::Int) = abs(x)/n

#
# Unbiased squared phase lag index, weighted phase
#
# See Vinck, M., Oostenveld, R., van Wingerden, M., Battaglia, F., &
# Pennartz, C. M. A. (2011). An improved index of phase-synchronization
# for electrophysiological data in the presence of volume-conduction,
# noise and sample-size bias. NeuroImage, 55(4), 1548–1565.
# doi:10.1016/j.neuroimage.2011.01.055
struct PLI2Unbiased <: PairwiseStatistic; end

accumulator(::Union{Type{PLI}, Type{PLI2Unbiased}}, ::Type{T}) where {T<:Real} = zero(Complex{T})
accumulate(::Union{Type{PLI}, Type{PLI2Unbiased}}, x::Complex{T},
           v1::Complex{T}, v2::Complex{T}) where {T<:Real} = (x + sign(imag(conj(v1)*v2)))
accumulate(::Union{Type{PLI}, Type{PLI2Unbiased}}, x::Complex{T},
           v1::Complex{T}, v2::Complex{T}, weight::Real) where {T<:Real} = (x + sign(imag(conj(v1)*v2))*weight)
finish(::Type{PLI2Unbiased}, x::Complex, n::Int) = abs2(x)/(n*(n-1)) - 1/(n-1)

#
# Weighted phase lag index
#
# See Vinck et al. (2011) as above.
struct WPLI <: PairwiseStatistic; end
struct WPLIAccumulator{T}
    si::Complex{T}   # Sum of imag(conj(v1)*(v2))
    sa::T            # Sum of abs(imag(conj(v1)*(v2)))
end

accumulator(::Type{WPLI}, ::Type{T}) where {T<:Real} =
    WPLIAccumulator{T}(zero(Complex{T}), zero(T))
function accumulate(::Type{WPLI}, x::WPLIAccumulator{T},
                    v1::Complex{T}, v2::Complex{T}) where T<:Real
    z = imag(conj(v1)*v2)
    WPLIAccumulator(x.si + z, x.sa + abs(z))
end
function accumulate(::Type{WPLI}, x::WPLIAccumulator{T},
                    v1::Complex{T}, v2::Complex{T}, weight::Real) where T<:Real
    z = imag(conj(v1)*v2)
    WPLIAccumulator(x.si + z*weight, x.sa + abs(z*weight))
end
finish(::Type{WPLI}, x::WPLIAccumulator, n::Int) = abs(x.si)/x.sa

#
# Debiased (i.e. still somewhat biased) WPLI^2
#
# See Vinck et al. (2011) as above.
struct WPLI2Debiased <: PairwiseStatistic; end
struct WPLI2DebiasedAccumulator{T}
    si::Complex{T}   # Sum of imag(conj(v1)*(v2))
    sa::T            # Sum of abs(imag(conj(v1)*(v2)))
    sa2::T           # Sum of abs2(imag(conj(v1)*(v2)))
end

accumulator(::Type{WPLI2Debiased}, ::Type{T}) where {T<:Real} =
    WPLI2DebiasedAccumulator{T}(zero(Complex{T}), zero(T), zero(T))
function accumulate(::Type{WPLI2Debiased}, x::WPLI2DebiasedAccumulator{T},
                    v1::Complex{T}, v2::Complex{T}) where T<:Real
    z = imag(conj(v1)*v2)
    WPLI2DebiasedAccumulator(x.si + z, x.sa + abs(z), x.sa2 + abs2(z))
end
function accumulate(::Type{WPLI2Debiased}, x::WPLI2DebiasedAccumulator{T},
                    v1::Complex{T}, v2::Complex{T}, weight::Real) where T<:Real
    z = imag(conj(v1)*v2)
    WPLI2DebiasedAccumulator(x.si + z*weight, x.sa + abs(z*weight), x.sa2 + abs2(z)*weight)
end
finish(::Type{WPLI2Debiased}, x::WPLI2DebiasedAccumulator,
       n::Int) = (abs2(x.si) - x.sa2)/(abs2(x.sa) - x.sa2)

#
# Functions applicable to all phase lag-style metrics
#
const PLStat = Union{PLI, PLI2Unbiased, WPLI, WPLI2Debiased}
Base.eltype(::PLStat, X::AbstractArray{Complex{T}}) where {T<:Real} = T
diagval(::Type{T}) where {T<:PLStat} = 0

allocwork(::PLStat, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}=X) where {T<:Complex} =
    nothing

# Single input matrix
function computestat!(t::S, out::AbstractMatrix{T}, work::Nothing,
                      X::AbstractVecOrMat{Complex{T}}) where {S<:PLStat, T<:Real}
    chkinput(out, X)
    for k = 1:size(X, 2), j = 1:k
        v = accumulator(S, T)
        @simd for i = 1:size(X, 1)
            @inbounds v = accumulate(S, v, X[i, j], X[i, k])
        end
        out[j, k] = finish(S, v, size(X, 1))
    end
    out
end

# Two input matrices
function computestat!(t::S, out::AbstractMatrix{T}, work::Nothing,
                      X::AbstractVecOrMat{Complex{T}},
                      Y::AbstractVecOrMat{Complex{T}}) where {S<:PLStat, T<:Real}
    chkinput(out, X, Y)
    for k = 1:size(Y, 2), j = 1:size(X, 2)
        v = accumulator(S, T)
        @simd for i = 1:size(X, 1)
            @inbounds v = accumulate(S, v, X[i, j], Y[i, k])
        end
        out[j, k] = finish(S, v, size(X, 1))
    end
    out
end


text:
# Demonstrate simulateDSM, (simulateSNR and predictSNR) => calcSNRInfo
using RSDeltaSigmaPort
using RSDeltaSigmaPort.EasyPlot #set, cons
import RSDeltaSigmaPort: BoundingBox
import Printf: @sprintf
j=im

#==Baseband modulator
===============================================================================#
println("\n*** 5th order, 2-level, baseband modulator")
OSR = 32
N = 8192
dsm = RealDSM(order=5, OSR=OSR, Hinf=1.5, opt=1)

@info("Performing ΔΣ simulation..."); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
ftest = dsm.f0 + 1/(3*OSR); fband = default_fband(dsm)
fband[1] += 2/N #Ignore first 2 low-frequency points
(u, iftest) = genTestTone(dsm, dbv(0.5), ftest, N=N)
NTF = synthesizeNTF(dsm)
simresult = simulateDSM(u, NTF)
println("\tdone.")

@info("Plotting modulator signals")
#-------------------------------------------------------------------------------
plot = plotModTransient(u, simresult.v)
	set(plot, xyaxes=set(xmin=0, xmax=300, ymin=-1.2, ymax=1.2))
saveimage(:png, "dsdemo2_o5_sig.png", plot, AR=2/1, width=900)
displaygui(plot)

@info("Plotting output spectrum (simulated vs theory)")
#-------------------------------------------------------------------------------
specinfo = calcSpecInfo(simresult.v, NTF, fband, ftest)
plot = plotModSpectrum(specinfo)
	plot.title="Modulator Output Spectrum @ OSR = $OSR."
saveimage(:png, "dsdemo2_o5_spec.png", plot, AR=2/1, width=900)
displaygui(plot)

@info("Plotting SNR vs input power")
#-------------------------------------------------------------------------------
snrinfo = calcSNRInfo(dsm, NTF=NTF)
plot = plotSNR(snrinfo, dsm)
	plot.title="SNR curve- theory and simulation"
saveimage(:png, "dsdemo2_o5_snr.png", plot, AR=2/1, width=900)
displaygui(plot)

#==Bandpass modulator
===============================================================================#
println("\n*** 8th order, 2-level, bandpass modulator")
OSR=64
f0 = 1/8
N = 8192
dsm = RealDSM(order=8, OSR=OSR, f0=f0, Hinf=1.5, opt=1)
display(dsm)

@info("Performing ΔΣ simulation..."); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
ftest = dsm.f0 + 1/(6*OSR); fband = default_fband(dsm)
(u, iftest) = genTestTone(dsm, dbv(0.5), ftest, N=N) #Half-scale sine-wave input
NTF = synthesizeNTF(dsm)
simresult = simulateDSM(u, NTF)
println("\tdone.")

@info("Plotting modulator signals")
#-------------------------------------------------------------------------------
plot = plotModTransient(u, simresult.v)
	set(plot, xyaxes=set(xmin=0, xmax=300, ymin=-1.2, ymax=1.2))
saveimage(:png, "dsdemo2_bp_o5_sig.png", plot, AR=2/1, width=900)
displaygui(plot)

@info("Plotting output spectrum (simulated vs theory)")
#-------------------------------------------------------------------------------
specinfo = calcSpecInfo(simresult.v, NTF, fband, ftest)
plot = plotModSpectrum(specinfo)
	plot.title="Modulator Output Spectrum @ OSR = $OSR."
saveimage(:png, "dsdemo2_bp_o5_spec.png", plot, AR=2/1, width=900)
displaygui(plot)

@info("Plotting SNR vs input power")
#-------------------------------------------------------------------------------
snrinfo = calcSNRInfo(dsm, NTF=NTF)
plot = plotSNR(snrinfo, dsm)
	plot.title="SNR curve- theory and simulation"
saveimage(:png, "dsdemo2_bp_o5_snr.png", plot, AR=2/1, width=900)
displaygui(plot)

#==15-step baseband modulator
===============================================================================#
println("\n*** 7th order, 15-step, baseband modulator")
OSR = 8
M = 16 #Shouldn'nt this be 15 (nlev=16?)
N = 8192
Hinf_list = [2.0, 8.0]
dsm_list = [RealDSM(order=7, OSR=OSR, M=M, opt=1, Hinf=Hinf) for Hinf in Hinf_list]
color_list = [:blue, :green]

@info("Performing ΔΣ simulation for H(∞)=$(Hinf_list)..."); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
id_list = [@sprintf("H(∞)=%.1f", Hinf) for Hinf in Hinf_list]
ftest = 1/(7*OSR); fband = default_fband(OSR)
fband[1] += 2/N #Ignore first 2 low-frequency points
(u, iftest) = genTestTone(dsm, dbv(0.5*M), ftest, N=N) #Half-scale sine-wave input
NTF_list = [synthesizeNTF(dsm) for dsm in dsm_list]
v_list = [simulateDSM(u, NTF, nlev=M+1).v for NTF in NTF_list] #simulateDSM(..).v is mod output, v.
println("\tdone.")

@info("Plotting input/output characteristics of ΔΣ simulations")
#-------------------------------------------------------------------------------
ioplotc = cons(:plot_collection, title="15-step / 7th-order / H(∞)=$(Hinf_list)")

#Plot input & output transients:
for (v, c) in zip(v_list, color_list) #Each simulated output
	local plot = plotModTransient(u, v, color=c, legend=false)
		set(plot, xyaxes=set(xmin=0, xmax=100, ymin=-16, ymax=16))
	push!(ioplotc, plot)
end

#Append SNR vs input curves:
snrinfo = calcSNRInfo(dsm_list[1], NTF=NTF_list[1])
plot = plotSNR(snrinfo, dsm_list[1], color=color_list[1])
snrinfo = calcSNRInfo(dsm_list[2], NTF=NTF_list[2])
plotSNR!(plot, snrinfo, dsm_list[2], color=color_list[2])
	set(plot, xyaxes=set(xmin=-100, xmax=0, ymin=0, ymax=120))
push!(ioplotc, plot)

#Specify plot locations to help readability:
ioplotc.bblist = [ #Format plot locations
	BoundingBox(0, 0.5, 0, 0.5), #1st modulator transient
	BoundingBox(0, 0.5, 0.5, 1), #2nd modulator transient
	BoundingBox(0.5, 1, 0, 1), #SNR curve
]
saveimage(:png, "dsdemo2_o7_15s_io.png", ioplotc, AR=2/1, width=900)
displaygui(ioplotc)

@info("Plotting output spectrum (simulated vs theory)")
#-------------------------------------------------------------------------------
plot = plotModSpectrum()
	plot.title="Modulator Output Spectrum @ OSR = $OSR."
	set(plot, xyaxes=set(ymin=-160))
for i in keys(NTF_list)
	local specinfo = calcSpecInfo(v_list[i], NTF_list[i], fband, ftest, M=M)
	plotModSpectrum!(plot, specinfo, id=id_list[i], color=color_list[i])
end
saveimage(:png, "dsdemo2_o7_15s_spec.png", plot, AR=2/1, width=900)
displaygui(plot)

:END_OF_DEMO


text:
@Base.kwdef struct Token
    special::Bool = false
    val::String = ""
    color::Color = Color()
    align::Int = 0
end
Token(special::Bool, val::AbstractString) = Token(special, val, Color(), 0)
Token(special::Bool, val::AbstractString, color::Color) = Token(special, val, color, 0)
function Token(s::AbstractString)
    if s in ["function", "filepath", "module"]
        return Token(true, s)
    end
    v = split(s, r"(?<!\\):")
    if v[1] == "frameno"
        if length(v) == 1
            return Token(true, s)
        else
            align = v[2] == "lalign" ? -1 : 1
            return Token(true, v[1], Color(), align)
        end
    end
    if length(v) == 1
        return Token(false, v[1])
    elseif v[2] in ["lalign", "ralign"]
        align = v[2] == "lalign" ? -1 : 1
        return Token(false, v[1], Color(v[3:end]), align)
    else
        return Token(false, v[1], Color(v[2:end]))
    end
end

const FRAME_LINE = [
    Token(false, " "),
    Token(true, "frameno", Color(), 1),
    Token(false, " "),
    Token(true, "function"),
    Token(false, "\n "),
    Token(false, "@", Color(), 1),
    Token(false, " "),
    Token(true, "module"),
    Token(true, "filepath")
]
function format_frameline(s::AbstractString, level = 0, tokens = Token[])
    i = 0
    while i < ncodeunits(s)
        i = nextind(s, i)
        c = s[i]
        if c == '{' && level == 0
            i > 1 && push!(tokens, Token(s[1:prevind(s, i)]))
            format_frameline(s[nextind(s, i):end], 1, tokens)
            break
        end

        if c == '}' && level == 1
            push!(tokens, Token(s[1:prevind(s, i)]))
            format_frameline(s[nextind(s, i):end], 0, tokens)
            break
        end
    end

    isempty(tokens) && return nothing
    empty!(FRAME_LINE)
    for token in tokens
        push!(FRAME_LINE, token)
    end

    return nothing
end


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $(realpath $(dirname $0))))/examples/" "${BASH_SOURCE[0]}" "$@" -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

include(joinpath(dirname(dirname(@__FILE__)), "src", "SpelledOut.jl"))

using .SpelledOut
using Test

@testset "SpelledOut.jl" begin
### ENGLISH TESTS
    @test spelled_out(1234, lang = :en) == "one thousand, two hundred thirty-four"
    @test spelled_out(1234) == "one thousand, two hundred thirty-four" # default to english
    @test spelled_out(1234, lang = :en, british=true) == "one thousand, two hundred and thirty-four"
    @test Spelled_out(1234, lang = :en) == "One thousand, two hundred thirty-four"
    @test Spelled_Out(1234, lang = :en) == "One Thousand, Two Hundred Thirty-Four"
    @test SPELLED_OUT(1234, lang = :en) == "ONE THOUSAND, TWO HUNDRED THIRTY-FOUR"
    @test spelled_out(123456789, lang = :en) == "one hundred twenty-three million, four hundred fifty-six thousand, seven hundred eighty-nine"
    @test spelled_out(123456789, lang = :en, british=true) == "one hundred and twenty-three million, four hundred and fifty-six thousand, seven hundred and eighty-nine"
    @test spelled_out(0, lang = :en) == "zero"
    @test spelled_out(0, lang = :en, british=true) == "zero"
    @test spelled_out(100, lang = :en) == "one hundred"
    @test spelled_out(100, lang = :en, british=true) == "one hundred"
    @test spelled_out(123456789, lang = :en, british=true, dict=:european) == "one hundred and twenty-three million, four hundred and fifty-six thousand, seven hundred and eighty-nine"
    @test spelled_out(123456789, lang = :en, dict=:european) == "one hundred twenty-three million, four hundred fifty-six thousand, seven hundred eighty-nine"
    @test spelled_out(123456789, lang = :en, british=true, dict=:british) == "one hundred and twenty-three million, four hundred and fifty-six thousand, seven hundred and eighty-nine"
    @test spelled_out(1234567890123, lang = :en, british=true, dict=:british) == "one billion, two hundred and thirty-four thousand million, five hundred and sixty-seven million, eight hundred and ninety thousand, one hundred and twenty-three"
    @test spelled_out(1234567890123, lang = :en, british=true, dict=:modern) == "one trillion, two hundred and thirty-four billion, five hundred and sixty-seven million, eight hundred and ninety thousand, one hundred and twenty-three"
    @test spelled_out(3.0, lang = :en) == "three"
    @test spelled_out(3, lang = :en) == "three"
    @test spelled_out(3.141592653, lang = :en) == "three point one four one five nine two six five three"
end


text:
function apply!(penalty::ProximalPenalty, model::JuMP.Model, β::Vector{JuMP.VariableRef}, ys::Vector{JuMP.VariableRef})
    β⁰ = penalty.β
    @assert length(β) == length(β⁰)
    λ = penalty.λ
    quad_penalty = JuMP.@expression(model, sum((β⁰[i] - β[i])^2 for i in 1:length(β)))
    return 1 / λ * quad_penalty
end

function apply!(penalty::L1Regularizer, model::JuMP.Model, β::Vector{JuMP.VariableRef}, ys::Vector{JuMP.VariableRef})
    num_vars = length(β)
    w = JuMP.@variable(model, [i=1:num_vars], lower_bound=0)
    for i in 1:num_vars
        JuMP.@constraint(model, w[i] ≥ β[i])
        JuMP.@constraint(model, w[i] ≥ -β[i])
    end
    return penalty.coefficient * sum(w)
end

function apply!(penalty::L2Regularizer, model::JuMP.Model, β::Vector{JuMP.VariableRef}, ys::Vector{JuMP.VariableRef})
    quad_penalty = JuMP.@expression(model, sum((β[i])^2 for i in 1:length(β)))
    return penalty.coefficient * quad_penalty
end

struct OptModel{T <: Union{JuMP.VariableRef, Float64}}
    jump_model::JuMP.Model
    linear_model::LinearModel
    vs::Vector{T}
    ys::Vector{T}
    zs::Matrix{T}
end

function optimize!(learner::Learner, opt_model::OptModel)
    JuMP.optimize!(opt_model.jump_model)
    primal_status = JuMP.primal_status(opt_model.jump_model)
    if primal_status != MOI.FEASIBLE_POINT
        @error "MIP did not produce a feasible solution"
    end
    termination_status = JuMP.termination_status(opt_model.jump_model)
    if termination_status != MOI.OPTIMAL
        @warn "MIP did not solve to optimality"
    end
    weight_values = JuMP.value.(opt_model.linear_model.weights)
    offset_value = JuMP.value(opt_model.linear_model.offset)
    v_values = JuMP.value.(opt_model.vs)
    y_values = JuMP.value.(opt_model.ys)
    z_values = JuMP.value.(opt_model.zs)
    linear_model_value = LinearModel(weight_values, offset_value)
    _audit_weights(learner, weight_values)
    return OptModel(opt_model.jump_model,
        linear_model_value,
        v_values,
        y_values,
        z_values)
end

function formulate_big_m(learner::Learner{T}, formulation_fn::Function=formulate_data_point_graph_big_m; relax_integrality=false) where {T <: Union{MIP,LPRelaxation}}
    model = JuMP.Model(learner.algorithm.optimizer_factory)
    train_data = learner.train_data
    num_features = feature_count(train_data)
    num_samples = sample_count(train_data)
    β = JuMP.@variable(model, [i=1:num_features], lower_bound=learner.param_lb[i], upper_bound=learner.param_ub[i])
    offset = JuMP.@variable(model, lower_bound=learner.offset_lb, upper_bound=learner.offset_ub)
    vs = Vector{JuMP.VariableRef}(undef, num_samples)
    ys = Vector{JuMP.VariableRef}(undef, num_samples)
    zs = Matrix{JuMP.VariableRef}(undef, num_samples, 3)
    for i in 1:num_samples
        sample = train_data.samples[i]
        features = sample.features
        v = JuMP.@variable(model)
        y = JuMP.@variable(model, lower_bound = 0, upper_bound = sample.bid_price)
        JuMP.@constraint(model, v == sum(features[k] * β[k] for k in 1:num_features) + offset)
        z = formulation_fn(model, v, y, sample.pay_price, sample.bid_price, M⁻(learner, features), M⁺(learner, features), relax_integrality)
        vs[i] = v
        ys[i] = y
        zs[i,:] = z
    end
    # ret_val = [formulation_fn(model, β, train_data.samples[i], L, U, relax_integrality)
            # for i in 1:num_samples]
    # ys = [ret_val[i][1] for i in 1:num_samples]
    # indicators = [ret_val[i][2] for i in 1:num_samples]
    obj = 1 / num_samples * sum(ys)
    for penalty in learner.objective_penalties
        obj -= apply!(penalty, model, β, ys)
    end
    JuMP.@objective(model, Max, obj)
    return OptModel(model, LinearModel(β, offset), vs, ys, zs)
end

M⁻(features, L, U) = sum(min.(L .* features, U .* features))
M⁺(features, L, U) = sum(max.(L .* features, U .* features))

function M⁻(learner, features)
    @assert (n = length(features)) == length(learner.param_lb) == length(learner.param_ub)
    lb = learner.offset_lb
    for i in 1:n
        lb += min(learner.param_lb[i] * features[i], learner.param_ub[i] * features[i])
    end
    return lb
end

function M⁺(learner, features)
    @assert (n = length(features)) == length(learner.param_lb) == length(learner.param_ub)
    ub = learner.offset_ub
    for i in 1:n
        ub += max(learner.param_lb[i] * features[i], learner.param_ub[i] * features[i])
    end
    return ub
end

function formulate_data_point_hypograph_big_m(model::JuMP.Model,
        v::JuMP.VariableRef,
        y::JuMP.VariableRef,
        pay_price::Float64,
        bid_price::Float64,
        l::Float64,
        u::Float64,
        relax_integrality::Bool)
    if pay_price != 0.0
        error("Hypograph formulation does not currently work with pay prices.")
    end
    z = if relax_integrality
        JuMP.@variable(model, lower_bound=0.0, upper_bound=1.0)
    else
        JuMP.@variable(model, binary=true)
    end
    JuMP.@constraints(model, begin
        y <= bid_price * z
        y <= v - l * (1 - z)
        v >= u * (1 - z)
        v <= bid_price + (u - bid_price) * (1 - z)
    end)
    return z
end

function formulate_data_point_graph_big_m(model::JuMP.Model,
    v::JuMP.VariableRef,
    y::JuMP.VariableRef,
    pay_price::Float64,
    bid_price::Float64,
    l::Float64,
    u::Float64,
    relax_integrality::Bool)
    z = (if relax_integrality
        JuMP.@variable(model, [1:3], lower_bound=0.0, upper_bound=1.0)
    else
        JuMP.@variable(model, [1:3], binary=true)
    end)
    JuMP.@constraint(model, sum(z) == 1)
    JuMP.@constraints(model, begin
        # D1 = {(x,y) : y == pay_price, v <= pay_price}
        # D2 = {(x,y) : y == v, pay_price <= v <= bid_price}
        # D3 = {(x,y) : y == 0, v >= bid_price}
        y <= pay_price * z[1] + bid_price * z[2] + 0 * z[3]
        y >= pay_price * z[1] + pay_price * z[2] + 0 * z[3]
        y <= v + (pay_price - l) * z[1] - bid_price * z[3]
        y >= v - u * z[3]
    end)
    return z
end

max_stars_repo_name: joehuchette/reserve-price-optimization

max_issues_repo_name: joehuchette/reserve-price-optimization

max_forks_repo_name: joehuchette/reserve-price-optimization


text:
#Written By Johnathan Bizzano
module Utils

    export Bound, Entry

    struct Bound{T <: Number}
        lower_bound::T
        upper_bound::T
        Bound{T}(lower_bound, upper_bound) where T = new{T}(T(lower_bound), T(upper_bound))
        Bound(lower_bound::T, upper_bound::T) where T = new{T}(lower_bound, upper_bound)
    end

    struct Entry{K, V}
        key::K
        value::V
        Entry{K, V}(key, value) where {K, V} = new{K, V}(K(key), V(value))
        Entry(key::K, value::V) where {K, V} = new{K, V}(K(key), V(value))
    end

    include("Fun.jl")
    include("JuliaUtils.jl")
    include("Iterable.jl")
end


text:
using SemiclassicalOrthogonalPolynomials, ClassicalOrthogonalPolynomials, LinearAlgebra
import ClassicalOrthogonalPolynomials: associated

ρ = 0.5
T = TwoBandJacobi(ρ,-1/2,-1/2,1/2)
U = associated(T)
x = axes(T,1)
H = inv.(x .- x')

L = U \ (H * Weighted(T))
c = L[:,2:∞] \ (U \ x)
c[1] = -c[2]

u = Weighted(T) * c
u[0.5000001]

[0.1]


text:
module day23

using CircularList, ProgressMeter

ex = "389125467"
input = "925176834"

# first cup in list is current cup (cups is already shifted)
import Base: pop!
function pop!(cups::CircularList.List{T}) where T
    val = current(cups).data; delete!(cups); forward!(cups)
    return val
end

function make_lookup(cups)
    lookup = Dict{Int, CircularList.Node{Int}}()
    for i in 1:length(cups)
        node = current(cups)
        lookup[node.data] = node
        forward!(cups)
    end
    return lookup
end

pickup3!(cups) = [pop!(cups) for _ in 1:3]

function place3!(cups, items, lookup)
    for item in items
        insert!(cups, item);
        lookup[item] = current(cups)
    end
end

function play(cups, moves; verbose = false)

    lo = minimum(cups)
    hi = maximum(cups)
    lookup = make_lookup(cups)

    @showprogress for move in 1:moves

        verbose && println("-- move $move --")
        verbose && println("cups: ", cups)

        curr = head(cups).data

        # pickup three cups immediately CW of current cup
        forward!(cups)       # move pointer CW by one
        pu3 = pickup3!(cups) # pickup those three cups

        verbose && println("pick up: ", pu3)

        # get desitination cup
        dest = curr - 1
        if dest < lo
            dest = hi
        end
        while dest ∈ pu3
            dest -= 1
            if dest < lo
                dest = hi
            end
        end
        verbose && println("destination: ", dest)
        # jump!(cups, dest)
        jump!(cups, lookup[dest])

        # Place three cups immediately CW of dest
        place3!(cups, pu3, lookup)

        # Select new curr cup. immediately CW of curr cup
        # jump!(cups, curr)   # advance pointer to curr
        jump!(cups, lookup[curr])   # advance pointer to curr
        forward!(cups)              # then move one more

        verbose && println()
    end

    jump!(cups, lookup[1]) # set pointer to cup #1 and return
    return cups
end

function part1(input)
    cups = circularlist([parse(Int, c) for c in input])
    res = play(cups, 100, verbose=false)
    # jump!(res, lookup[1])
    return join(collect(res)[2:end])
end

@assert part1(ex) == "67384529"
@info "Day 23, Part 1 answer: $(part1(input))"

function part2(input)
    c1 = [parse(Int, c) for c in input]
    cups = circularlist(vcat(c1, collect(10:1_000_000)));
    res = play(cups, 10_000_000, verbose=false)
    return res.current.next.data * res.current.next.next.data
    # return res
end

@info "Day 23, Part 2 answer: 91408386135"

""" Takes ~40s to solve part 2.
julia> part2(input)
Progress: 100%|██████████████████████████████████████████████████| Time: 0:00:41
91408386135
"""

end


text:
__precompile__()
module Length

import DataStructures: counter
using SeqHax.Utils
using ArgParse

function add_args(argparser)
    add_common_args(argparser["length"])
    add_paired_args(argparser["length"])
    @add_arg_table argparser["length"] begin
        "input"
            help = "Input FASTQ read file(s) (may be gzipped)"
            required = true
            arg_type = String
    end
end

function main(args)
    r1ctr = counter(Int)
    paired = args["paired"]
    quiet = args["quiet"]
    r2ctr = ifelse(paired, counter(Int), r1ctr)

    if paired
        npair = foreach_read(args["input"], quiet=quiet) do idx, read
            push!(r1ctr, length(read.seq))
        end
    else
        npair = foreach_readpair(args["input"], quiet=quiet) do idx, r1, r2
            push!(r1ctr, length(r1.seq))
            push!(r2ctr, length(r2.seq))
        end
    end

    lengths = sort(union(keys(r1ctr), keys(r2ctr)))
    if paired
        println("read_length\tr1_count\tr2_count")
    else
        println("read_length\tcount")
    end
    for len in lengths
        if paired
            r1 = r1ctr[len]
            r2 = r2ctr[len]
            println("$len\t$r1\t$r2")
        else
            cnt = r1ctr[len]
            println("$len\t$cnt")
        end
    end
end

end # module Length


text:
function banach_rounding_order(weights::GeometricWeights{T}, bound::T) where {T<:AbstractFloat}
    isone(rate(weights)) | isinf(bound) && return typemax(Int)
    ϵ = eps(T)
    bound ≤ ϵ && return 0
    order = log(bound/ϵ)/log(rate(weights))
    isinf(order) && return typemax(Int)
    return ceil(Int, order)
end
function banach_rounding_order(weights::AlgebraicWeights{T}, bound::T) where {T<:AbstractFloat}
    iszero(rate(weights)) | isinf(bound) && return typemax(Int)
    ϵ = eps(T)
    bound ≤ ϵ && return 0
    order = exp(log(bound/ϵ)/rate(weights))-1
    isinf(order) && return typemax(Int)
    return ceil(Int, order)
end

function banach_rounding_order(weights::GeometricWeights, bound::Real)
    rate_, bound_ = promote(float(rate(weights)), float(bound))
    return banach_rounding_order(GeometricWeights(rate_), bound_)
end
function banach_rounding_order(weights::AlgebraicWeights, bound::Real)
    rate_, bound_ = promote(float(rate(weights)), float(bound))
    return banach_rounding_order(AlgebraicWeights(rate_), bound_)
end

banach_rounding_order(weights::GeometricWeights{<:Interval}, bound::Interval) =
    banach_rounding_order(GeometricWeights(sup(rate(weights))), sup(bound))
banach_rounding_order(weights::AlgebraicWeights{<:Interval}, bound::Interval) =
    banach_rounding_order(AlgebraicWeights(sup(rate(weights))), sup(bound))

banach_rounding_order(weights::GeometricWeights{<:Interval}, bound::Real) =
    banach_rounding_order(GeometricWeights(sup(rate(weights))), bound)
banach_rounding_order(weights::AlgebraicWeights{<:Interval}, bound::Real) =
    banach_rounding_order(AlgebraicWeights(sup(rate(weights))), bound)

banach_rounding_order(weights::GeometricWeights, bound::Interval) =
    banach_rounding_order(weights, sup(bound))
banach_rounding_order(weights::AlgebraicWeights, bound::Interval) =
    banach_rounding_order(weights, sup(bound))

banach_rounding_order(weights::Tuple, bound::Real) =
    map(νᵢ -> banach_rounding_order(νᵢ, bound), weights)

function banach_rounding_order(weights::Tuple, bound::Interval)
    sup_bound = sup(bound)
    return map(νᵢ -> banach_rounding_order(νᵢ, sup_bound), weights)
end

#

function banach_rounding!(a, weights, bound)
    banach_rounding!(a, weights, bound, banach_rounding_order(weights, bound))
    return a
end

function banach_rounding!(a::Sequence{TensorSpace{T},<:AbstractVector{S}}, weights::NTuple{N,Weights}, bound::Real, rounding_order::NTuple{N,Int}) where {N,T<:NTuple{N,BaseSpace},S<:Interval}
    M = typemax(Int)
    w = map(ord -> ifelse(ord == 0, 1, ord), rounding_order)
    @inbounds for α ∈ indices(space(a))
        if mapreduce((wᵢ, αᵢ) -> wᵢ == M ? 0//1 : abs(αᵢ) // wᵢ, +, w, α) ≥ 1
            μᵅ = bound / weight(weights, α)
            sup_μᵅ = sup(μᵅ)
            a[α] = Interval(-sup_μᵅ, sup_μᵅ)
        end
    end
    return a
end

function banach_rounding!(a::Sequence{TensorSpace{T},<:AbstractVector{Complex{S}}}, weights::NTuple{N,Weights}, bound::Real, rounding_order::NTuple{N,Int}) where {N,T<:NTuple{N,BaseSpace},S<:Interval}
    M = typemax(Int)
    w = map(ord -> ifelse(ord == 0, 1, ord), rounding_order)
    @inbounds for α ∈ indices(space(a))
        if mapreduce((wᵢ, αᵢ) -> wᵢ == M ? 0//1 : abs(αᵢ) // wᵢ, +, w, α) ≥ 1
            μᵅ = bound / weight(weights, α)
            sup_μᵅ = sup(μᵅ)
            x = Interval(-sup_μᵅ, sup_μᵅ)
            a[α] = complex(x, x)
        end
    end
    return a
end

# Taylor

function banach_rounding!(a::Sequence{Taylor,<:AbstractVector{T}}, weights::GeometricWeights, bound::Real, rounding_order::Int) where {T<:Interval}
    if rounding_order ≤ order(a)
        ν⁻¹ = inv(rate(weights))
        μⁱ = bound / weight(weights, rounding_order)
        @inbounds for i ∈ rounding_order:order(a)
            sup_μⁱ = sup(μⁱ)
            a[i] = Interval(-sup_μⁱ, sup_μⁱ)
            μⁱ *= ν⁻¹
        end
    end
    return a
end
function banach_rounding!(a::Sequence{Taylor,<:AbstractVector{T}}, weights::AlgebraicWeights, bound::Real, rounding_order::Int) where {T<:Interval}
    @inbounds for i ∈ rounding_order:order(a)
        sup_μⁱ = sup(bound / weight(weights, i))
        a[i] = Interval(-sup_μⁱ, sup_μⁱ)
    end
    return a
end

function banach_rounding!(a::Sequence{Taylor,<:AbstractVector{Complex{T}}}, weights::GeometricWeights, bound::Real, rounding_order::Int) where {T<:Interval}
    if rounding_order ≤ order(a)
        ν⁻¹ = inv(rate(weights))
        μⁱ = bound / weight(weights, rounding_order)
        @inbounds for i ∈ rounding_order:order(a)
            sup_μⁱ = sup(μⁱ)
            x = Interval(-sup_μⁱ, sup_μⁱ)
            a[i] = complex(x, x)
            μⁱ *= ν⁻¹
        end
    end
    return a
end
function banach_rounding!(a::Sequence{Taylor,<:AbstractVector{Complex{T}}}, weights::AlgebraicWeights, bound::Real, rounding_order::Int) where {T<:Interval}
    @inbounds for i ∈ rounding_order:order(a)
        sup_μⁱ = sup(bound / weight(weights, i))
        x = Interval(-sup_μⁱ, sup_μⁱ)
        a[i] = complex(x, x)
    end
    return a
end

# Fourier

function banach_rounding!(a::Sequence{<:Fourier,<:AbstractVector{T}}, weights::GeometricWeights, bound::Real, rounding_order::Int) where {T<:Interval}
    if rounding_order ≤ order(a)
        ν⁻¹ = inv(rate(weights))
        μⁱ = bound / weight(weights, rounding_order)
        @inbounds for i ∈ rounding_order:order(a)
            sup_μⁱ = sup(μⁱ)
            x = Interval(-sup_μⁱ, sup_μⁱ)
            a[i] = x
            a[-i] = x
            μⁱ *= ν⁻¹
        end
    end
    return a
end
function banach_rounding!(a::Sequence{<:Fourier,<:AbstractVector{T}}, weights::AlgebraicWeights, bound::Real, rounding_order::Int) where {T<:Interval}
    @inbounds for i ∈ rounding_order:order(a)
        sup_μⁱ = sup(bound / weight(weights, i))
        x = Interval(-sup_μⁱ, sup_μⁱ)
        a[i] = x
        a[-i] = x
    end
    return a
end

function banach_rounding!(a::Sequence{<:Fourier,<:AbstractVector{Complex{T}}}, weights::GeometricWeights, bound::Real, rounding_order::Int) where {T<:Interval}
    if rounding_order ≤ order(a)
        ν⁻¹ = inv(rate(weights))
        μⁱ = bound / weight(weights, rounding_order)
        @inbounds for i ∈ rounding_order:order(a)
            sup_μⁱ = sup(μⁱ)
            x = Interval(-sup_μⁱ, sup_μⁱ)
            complex_x = complex(x, x)
            a[i] = complex_x
            a[-i] = complex_x
            μⁱ *= ν⁻¹
        end
    end
    return a
end
function banach_rounding!(a::Sequence{<:Fourier,<:AbstractVector{Complex{T}}}, weights::AlgebraicWeights, bound::Real, rounding_order::Int) where {T<:Interval}
    @inbounds for i ∈ rounding_order:order(a)
        sup_μⁱ = sup(bound / weight(weights, i))
        x = Interval(-sup_μⁱ, sup_μⁱ)
        complex_x = complex(x, x)
        a[i] = complex_x
        a[-i] = complex_x
    end
    return a
end

# Chebyshev

function banach_rounding!(a::Sequence{Chebyshev,<:AbstractVector{T}}, weights::GeometricWeights, bound::Real, rounding_order::Int) where {T<:Interval}
    if rounding_order ≤ order(a)
        ν⁻¹ = inv(rate(weights))
        μⁱ = bound / weight(weights, rounding_order)
        @inbounds for i ∈ rounding_order:order(a)
            sup_μⁱ = sup(μⁱ)
            a[i] = Interval(-sup_μⁱ, sup_μⁱ)
            μⁱ *= ν⁻¹
        end
    end
    return a
end
function banach_rounding!(a::Sequence{Chebyshev,<:AbstractVector{T}}, weights::AlgebraicWeights, bound::Real, rounding_order::Int) where {T<:Interval}
    @inbounds for i ∈ rounding_order:order(a)
        sup_μⁱ = sup(bound / weight(weights, i))
        a[i] = Interval(-sup_μⁱ, sup_μⁱ)
    end
    return a
end

function banach_rounding!(a::Sequence{Chebyshev,<:AbstractVector{Complex{T}}}, weights::GeometricWeights, bound::Real, rounding_order::Int) where {T<:Interval}
    if rounding_order ≤ order(a)
        ν⁻¹ = inv(rate(weights))
        μⁱ = bound / weight(weights, rounding_order)
        @inbounds for i ∈ rounding_order:order(a)
            sup_μⁱ = sup(μⁱ)
            x = Interval(-sup_μⁱ, sup_μⁱ)
            a[i] = complex(x, x)
            μⁱ *= ν⁻¹
        end
    end
    return a
end
function banach_rounding!(a::Sequence{Chebyshev,<:AbstractVector{Complex{T}}}, weights::AlgebraicWeights, bound::Real, rounding_order::Int) where {T<:Interval}
    @inbounds for i ∈ rounding_order:order(a)
        sup_μⁱ = sup(bound / weight(weights, i))
        x = Interval(-sup_μⁱ, sup_μⁱ)
        a[i] = complex(x, x)
    end
    return a
end

#

_inf_max_order(::BaseSpace) = typemax(Int)
_inf_max_order(::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} = ntuple(i -> typemax(Int), Val(N))

function Base.:*(a::Sequence{<:SequenceSpace}, b::Sequence{<:SequenceSpace})
    new_space = image(*, space(a), space(b))
    CoefType = promote_type(eltype(a), eltype(b))
    c = Sequence(new_space, zeros(CoefType, dimension(new_space)))
    _add_mul!(c, a, b, _inf_max_order(new_space))
    return c
end

function banach_rounding_mul(a::Sequence{<:BaseSpace}, b::Sequence{<:BaseSpace}, weights::Weights)
    weighted_ℓ¹_norm = Weightedℓ¹Norm(weights)
    bound_ab = norm(a, weighted_ℓ¹_norm) * norm(b, weighted_ℓ¹_norm)
    rounding_order = banach_rounding_order(weights, bound_ab)
    new_space = image(*, space(a), space(b))
    CoefType = promote_type(eltype(a), eltype(b))
    c = Sequence(new_space, zeros(CoefType, dimension(new_space)))
    banach_rounding!(c, weights, bound_ab, rounding_order)
    _add_mul!(c, a, b, rounding_order-1)
    return c
end

function banach_rounding_mul(a::Sequence{TensorSpace{T}}, b::Sequence{TensorSpace{S}}, weights::NTuple{N,Weights}) where {N,T<:NTuple{N,BaseSpace},S<:NTuple{N,BaseSpace}}
    weighted_ℓ¹_norm = Weightedℓ¹Norm(weights)
    bound_ab = norm(a, weighted_ℓ¹_norm) * norm(b, weighted_ℓ¹_norm)
    rounding_order = banach_rounding_order(weights, bound_ab)
    new_space = image(*, space(a), space(b))
    CoefType = promote_type(eltype(a), eltype(b))
    c = Sequence(new_space, zeros(CoefType, dimension(new_space)))
    banach_rounding!(c, weights, bound_ab, rounding_order)
    _add_mul!(c, a, b, ntuple(i -> rounding_order[i]-1, Val(N)))
    return c
end

function _add_mul!(c::Sequence{<:BaseSpace}, a, b, max_order)
    C = coefficients(c)
    A = coefficients(a)
    B = coefficients(b)
    _add_mul!(C, A, B, space(c), space(a), space(b), max_order)
    return c
end

function _add_mul!(c::Sequence{<:TensorSpace}, a, b, max_order)
    space_c = space(c)
    space_a = space(a)
    space_b = space(b)
    C = _no_alloc_reshape(coefficients(c), dimensions(space_c))
    A = _no_alloc_reshape(coefficients(a), dimensions(space_a))
    B = _no_alloc_reshape(coefficients(b), dimensions(space_b))
    _add_mul!(C, A, B, space_c, space_a, space_b, max_order, 0//1)
    return c
end

function _add_mul!(C, A, B, space_c::BaseSpace, space_a, space_b, max_order)
    ord = min(order(space_c), max_order)
    for i ∈ _convolution_indices(space_c, ord)
        _convolution!(C, A, B, space_c, space_a, space_b, i)
    end
    return C
end

function _add_mul!(C, A, B, space_c::TensorSpace{<:NTuple{N,BaseSpace}}, space_a, space_b, max_order, n) where {N}
    remaining_space_c = Base.front(space_c)
    remaining_space_a = Base.front(space_a)
    remaining_space_b = Base.front(space_b)
    remaining_max_order = Base.front(max_order)
    @inbounds current_space_c = space_c[N]
    @inbounds current_space_a = space_a[N]
    @inbounds current_space_b = space_b[N]
    @inbounds current_max_order = max_order[N]
    current_order_c = order(current_space_c)
    if current_max_order == typemax(Int)
        for i ∈ _convolution_indices(current_space_c, current_order_c)
            _convolution!(C, A, B, current_space_c, current_space_a, current_space_b, remaining_space_c, remaining_space_a, remaining_space_b, remaining_max_order, n, i)
        end
    else
        ord = min(current_order_c, floor(Int, current_max_order * (1 - n)))
        w = ifelse(current_max_order == 0, 1, current_max_order)
        for i ∈ _convolution_indices(current_space_c, ord)
            n_ = n + abs(i) // w
            _convolution!(C, A, B, current_space_c, current_space_a, current_space_b, remaining_space_c, remaining_space_a, remaining_space_b, remaining_max_order, n_, i)
        end
    end
    return C
end

function _add_mul!(C, A, B, space_c::TensorSpace{<:Tuple{BaseSpace}}, space_a, space_b, max_order, n)
    @inbounds current_space_c = space_c[1]
    @inbounds current_space_a = space_a[1]
    @inbounds current_space_b = space_b[1]
    @inbounds current_max_order = max_order[1]
    current_order_c = order(current_space_c)
    ord = current_max_order == typemax(Int) ? current_order_c : min(current_order_c, floor(Int, current_max_order * (1 - n)))
    for i ∈ _convolution_indices(current_space_c, ord)
        _convolution!(C, A, B, current_space_c, current_space_a, current_space_b, i)
    end
    return C
end

# Taylor

_convolution_indices(::Taylor, order) = 0:order

function _convolution!(C, A, B, space_c::Taylor, space_a::Taylor, space_b::Taylor, i)
    order_a = order(space_a)
    order_b = order(space_b)
    @inbounds for j ∈ max(i-order_a, 0):min(i, order_b)
        C[i+1] += A[i-j+1] * B[j+1]
    end
    return C
end

function _convolution!(C::AbstractArray{T,N}, A, B, current_space_c::Taylor, current_space_a::Taylor, current_space_b::Taylor, remaining_space_c, remaining_space_a, remaining_space_b, remaining_max_order, n, i) where {T,N}
    current_order_a = order(current_space_a)
    current_order_b = order(current_space_b)
    @inbounds Cᵢ = selectdim(C, N, i+1)
    @inbounds for j ∈ max(i-current_order_a, 0):min(i, current_order_b)
        _add_mul!(Cᵢ, selectdim(A, N, i-j+1), selectdim(B, N, j+1), remaining_space_c, remaining_space_a, remaining_space_b, remaining_max_order, n)
    end
    return C
end

# Fourier

_convolution_indices(::Fourier, order) = -order:order

function _convolution!(C, A, B, space_c::Fourier, space_a::Fourier, space_b::Fourier, i)
    order_c = order(space_c)
    order_a = order(space_a)
    order_b = order(space_b)
    @inbounds for j ∈ max(i-order_a, -order_b):min(i+order_a, order_b)
        C[i+order_c+1] += A[i-j+order_a+1] * B[j+order_b+1]
    end
    return C
end

function _convolution!(C::AbstractArray{T,N}, A, B, current_space_c::Fourier, current_space_a::Fourier, current_space_b::Fourier, remaining_space_c, remaining_space_a, remaining_space_b, remaining_max_order, n, i) where {T,N}
    current_order_c = order(current_space_c)
    current_order_a = order(current_space_a)
    current_order_b = order(current_space_b)
    @inbounds Cᵢ = selectdim(C, N, i+current_order_c+1)
    @inbounds for j ∈ max(i-current_order_a, -current_order_b):min(i+current_order_a, current_order_b)
        _add_mul!(Cᵢ, selectdim(A, N, i-j+current_order_a+1), selectdim(B, N, j+current_order_b+1), remaining_space_c, remaining_space_a, remaining_space_b, remaining_max_order, n)
    end
    return C
end

# Chebyshev

_convolution_indices(::Chebyshev, order) = 0:order

function _convolution!(C, A, B, space_c::Chebyshev, space_a::Chebyshev, space_b::Chebyshev, i)
    order_a = order(space_a)
    order_b = order(space_b)
    @inbounds for j ∈ max(i-order_a, -order_b):min(i+order_a, order_b)
        C[i+1] += A[abs(i-j)+1] * B[abs(j)+1]
    end
    return C
end

function _convolution!(C::AbstractArray{T,N}, A, B, current_space_c::Chebyshev, current_space_a::Chebyshev, current_space_b::Chebyshev, remaining_space_c, remaining_space_a, remaining_space_b, remaining_max_order, n, i) where {T,N}
    current_order_a = order(current_space_a)
    current_order_b = order(current_space_b)
    @inbounds Cᵢ = selectdim(C, N, i+1)
    @inbounds for j ∈ max(i-current_order_a, -current_order_b):min(i+current_order_a, current_order_b)
        _add_mul!(Cᵢ, selectdim(A, N, abs(i-j)+1), selectdim(B, N, abs(j)+1), remaining_space_c, remaining_space_a, remaining_space_b, remaining_max_order, n)
    end
    return C
end

#

function Base.:^(a::Sequence{<:SequenceSpace}, n::Int)
    n < 0 && return throw(DomainError(n, "^ is only defined for positive integers"))
    n == 0 && return one(a)
    n == 1 && return copy(a)
    n == 2 && return _sqr(a)
    # power by squaring
    t = trailing_zeros(n) + 1
    n >>= t
    while (t -= 1) > 0
        a = _sqr(a)
    end
    c = a
    while n > 0
        t = trailing_zeros(n) + 1
        n >>= t
        while (t -= 1) ≥ 0
            a = _sqr(a)
        end
        c *= a
    end
    return c
end

function banach_rounding_pow(a::Sequence{<:SequenceSpace}, n::Int, weights)
    n < 0 && return throw(DomainError(n, "^ is only defined for positive integers"))
    n == 0 && return one(a)
    n == 1 && return copy(a)
    n == 2 && return _banach_rounding_sqr(a, weights)
    # power by squaring
    t = trailing_zeros(n) + 1
    n >>= t
    while (t -= 1) > 0
        a = _banach_rounding_sqr(a, weights)
    end
    c = a
    while n > 0
        t = trailing_zeros(n) + 1
        n >>= t
        while (t -= 1) ≥ 0
            a = _banach_rounding_sqr(a, weights)
        end
        c = banach_rounding_mul(c, a, weights)
    end
    return c
end

function _sqr(a::Sequence{<:SequenceSpace})
    new_space = image(^, space(a), 2)
    CoefType = eltype(a)
    c = Sequence(new_space, zeros(CoefType, dimension(new_space)))
    _add_sqr!(c, a, _inf_max_order(new_space))
    return c
end

function _banach_rounding_sqr(a::Sequence{<:BaseSpace}, weights::Weights)
    weighted_ℓ¹_norm = Weightedℓ¹Norm(weights)
    norm_a = norm(a, weighted_ℓ¹_norm)
    bound_a² = norm_a * norm_a
    rounding_order = banach_rounding_order(weights, bound_a²)
    new_space = image(^, space(a), 2)
    CoefType = eltype(a)
    c = Sequence(new_space, zeros(CoefType, dimension(new_space)))
    banach_rounding!(c, weights, bound_a², rounding_order)
    _add_sqr!(c, a, rounding_order-1)
    return c
end

function _banach_rounding_sqr(a::Sequence{TensorSpace{T}}, weights::NTuple{N,Weights}) where {N,T<:NTuple{N,BaseSpace}}
    weighted_ℓ¹_norm = Weightedℓ¹Norm(weights)
    norm_a = norm(a, weighted_ℓ¹_norm)
    bound_a² = norm_a * norm_a
    rounding_order = banach_rounding_order(weights, bound_a²)
    new_space = image(^, space(a), 2)
    CoefType = eltype(a)
    c = Sequence(new_space, zeros(CoefType, dimension(new_space)))
    banach_rounding!(c, weights, bound_a², rounding_order)
    _add_sqr!(c, a, ntuple(i -> rounding_order[i]-1, Val(N)))
    return c
end

function _add_sqr!(c::Sequence{<:BaseSpace}, a, max_order)
    C = coefficients(c)
    A = coefficients(a)
    _add_sqr!(C, A, space(c), space(a), max_order)
    return c
end

function _add_sqr!(c::Sequence{<:TensorSpace}, a, max_order)
    space_c = space(c)
    space_a = space(a)
    C = _no_alloc_reshape(coefficients(c), dimensions(space_c))
    A = _no_alloc_reshape(coefficients(a), dimensions(space_a))
    _add_mul!(C, A, A, space_c, space_a, space_a, max_order, 0//1)
    return c
end

# Taylor

function _add_sqr!(C, A, space_c::Taylor, space_a::Taylor, max_order)
    order_a = order(space_a)
    @inbounds A₀ = A[1]
    @inbounds C[1] = A₀ * A₀
    ord = min(order(space_c), max_order)
    @inbounds for i ∈ 1:ord
        i_odd = i%2
        i_end = (i-2+i_odd)÷2
        @inbounds for j ∈ max(i-order_a, 0):i_end
            C[i+1] += A[i-j+1] * A[j+1]
        end
        C[i+1] *= 2
        if iszero(i_odd)
            A_i½ = A[i÷2+1]
            C[i+1] += A_i½ * A_i½
        end
    end
    return C
end

# Fourier

function _add_sqr!(C, A, space_c::Fourier, space_a::Fourier, max_order)
    order_c = order(space_c)
    order_a = order(space_a)
    @inbounds for j ∈ 1:order_a
        C[order_c+1] += A[j+order_a+1] * A[-j+order_a+1]
    end
    @inbounds A₀ = A[order_a+1]
    @inbounds C[order_c+1] = 2C[order_c+1] + A₀ * A₀
    ord = min(order_c, max_order)
    @inbounds for i ∈ 1:ord
        i½, i_odd = divrem(i, 2)
        @inbounds for j ∈ i½+1:order_a
            C[i+order_c+1] += A[i-j+order_a+1] * A[j+order_a+1]
            C[-i+order_c+1] += A[j-i+order_a+1] * A[-j+order_a+1]
        end
        C[i+order_c+1] *= 2
        C[-i+order_c+1] *= 2
        if iszero(i_odd)
            A_i½ = A[i½+order_a+1]
            A_neg_i½ = A[-i½+order_a+1]
            C[i+order_c+1] += A_i½ * A_i½
            C[-i+order_c+1] += A_neg_i½ * A_neg_i½
        end
    end
    return C
end

# Chebyshev

function _add_sqr!(C, A, space_c::Chebyshev, space_a::Chebyshev, max_order)
    order_a = order(space_a)
    @inbounds for j ∈ 1:order_a
        Aⱼ = A[j+1]
        C[1] += Aⱼ * Aⱼ
    end
    @inbounds A₀ = A[1]
    @inbounds C[1] = 2C[1] + A₀ * A₀
    ord = min(order(space_c), max_order)
    @inbounds for i ∈ 1:ord
        i½, i_odd = divrem(i, 2)
        @inbounds for j ∈ i½+1:order_a
            C[i+1] += A[abs(i-j)+1] * A[j+1]
        end
        C[i+1] *= 2
        if iszero(i_odd)
            A_i½ = A[i½+1]
            C[i+1] += A_i½ * A_i½
        end
    end
    return C
end

#

function *̄(a::Sequence{<:SequenceSpace}, b::Sequence{<:SequenceSpace})
    new_space = image(*̄, space(a), space(b))
    CoefType = promote_type(eltype(a), eltype(b))
    c = Sequence(new_space, zeros(CoefType, dimension(new_space)))
    _add_mul!(c, a, b, _inf_max_order(new_space))
    return c
end

function banach_rounding_mul_bar(a::Sequence{<:BaseSpace}, b::Sequence{<:BaseSpace}, weights::Weights)
    weighted_ℓ¹_norm = Weightedℓ¹Norm(weights)
    bound_ab = norm(a, weighted_ℓ¹_norm) * norm(b, weighted_ℓ¹_norm)
    rounding_order = banach_rounding_order(weights, bound_ab)
    new_space = image(*̄, space(a), space(b))
    CoefType = promote_type(eltype(a), eltype(b))
    c = Sequence(new_space, zeros(CoefType, dimension(new_space)))
    banach_rounding!(c, weights, bound_ab, rounding_order)
    _add_mul!(c, a, b, rounding_order-1)
    return c
end

function banach_rounding_mul_bar(a::Sequence{TensorSpace{T}}, b::Sequence{TensorSpace{S}}, weights::NTuple{N,Weights}) where {N,T<:NTuple{N,BaseSpace},S<:NTuple{N,BaseSpace}}
    weighted_ℓ¹_norm = Weightedℓ¹Norm(weights)
    bound_ab = norm(a, weighted_ℓ¹_norm) * norm(b, weighted_ℓ¹_norm)
    rounding_order = banach_rounding_order(weights, bound_ab)
    new_space = image(*̄, space(a), space(b))
    CoefType = promote_type(eltype(a), eltype(b))
    c = Sequence(new_space, zeros(CoefType, dimension(new_space)))
    banach_rounding!(c, weights, bound_ab, rounding_order)
    _add_mul!(c, a, b, ntuple(i -> rounding_order[i]-1, Val(N)))
    return c
end

#

function ^̄(a::Sequence{<:SequenceSpace}, n::Int)
    n < 0 && return throw(DomainError(n, "^̄ is only defined for positive integers"))
    n == 0 && return one(a)
    n == 1 && return copy(a)
    n == 2 && return _sqr_bar(a)
    # power by squaring
    t = trailing_zeros(n) + 1
    n >>= t
    while (t -= 1) > 0
        a = _sqr_bar(a)
    end
    c = a
    while n > 0
        t = trailing_zeros(n) + 1
        n >>= t
        while (t -= 1) ≥ 0
            a = _sqr_bar(a)
        end
        c = c *̄ a
    end
    return c
end

function banach_rounding_pow_bar(a::Sequence{<:SequenceSpace}, n::Int, weights)
    n < 0 && return throw(DomainError(n, "^ is only defined for positive integers"))
    n == 0 && return one(a)
    n == 1 && return copy(a)
    n == 2 && return _banach_rounding_sqr_bar(a, weights)
    # power by squaring
    t = trailing_zeros(n) + 1
    n >>= t
    while (t -= 1) > 0
        a = _banach_rounding_sqr_bar(a, weights)
    end
    c = a
    while n > 0
        t = trailing_zeros(n) + 1
        n >>= t
        while (t -= 1) ≥ 0
            a = _banach_rounding_sqr_bar(a, weights)
        end
        c = banach_rounding_mul_bar(c, a, weights)
    end
    return c
end

function _sqr_bar(a::Sequence{<:SequenceSpace})
    new_space = image(^̄, space(a), 2)
    CoefType = eltype(a)
    c = Sequence(new_space, zeros(CoefType, dimension(new_space)))
    _add_sqr!(c, a, _inf_max_order(new_space))
    return c
end

function _banach_rounding_sqr_bar(a::Sequence{<:BaseSpace}, weights::Weights)
    weighted_ℓ¹_norm = Weightedℓ¹Norm(weights)
    norm_a = norm(a, weighted_ℓ¹_norm)
    bound_a² = norm_a * norm_a
    rounding_order = banach_rounding_order(weights, bound_a²)
    new_space = image(^̄, space(a), 2)
    CoefType = eltype(a)
    c = Sequence(new_space, zeros(CoefType, dimension(new_space)))
    banach_rounding!(c, weights, bound_a², rounding_order)
    _add_sqr!(c, a, rounding_order-1)
    return c
end

function _banach_rounding_sqr_bar(a::Sequence{TensorSpace{T}}, weights::NTuple{N,Weights}) where {N,T<:NTuple{N,BaseSpace}}
    weighted_ℓ¹_norm = Weightedℓ¹Norm(weights)
    norm_a = norm(a, weighted_ℓ¹_norm)
    bound_a² = norm_a * norm_a
    rounding_order = banach_rounding_order(weights, bound_a²)
    new_space = image(^̄, space(a), 2)
    CoefType = eltype(a)
    c = Sequence(new_space, zeros(CoefType, dimension(new_space)))
    banach_rounding!(c, weights, bound_a², rounding_order)
    _add_sqr!(c, a, ntuple(i -> rounding_order[i]-1, Val(N)))
    return c
end


text:
using Base.Test
using Lora

println("    Testing VanillaMCTune constructors and methods...")

tune = VanillaMCTune()

@test tune.accepted == 0
@test tune.proposed == 0
@test isnan(tune.rate)

tune.proposed = 10
for i in 1:4
  count!(tune)
end
rate!(tune)

@test tune.accepted == 4
@test tune.proposed == 10
@test tune.rate == 0.4

reset!(tune)

@test tune.accepted == 0
@test tune.proposed == 0
@test isnan(tune.rate)


text:
cd(@__DIR__)
Pkg.activate("..")
using Chordal
using JuMP, Hypatia, COSMO
using SparseArrays, LinearAlgebra, Random

const CD = Chordal

rand_seed = 1234
n = 200
# optimizer = Hypatia.Optimizer
optimizer = optimizer_with_attributes(COSMO.Optimizer, "decompose" => false)

# min c^Tx st ΣF_ix_i + G ∈ PSDCONE()
# D is dual variable
c, F, G, ystar, D = CD.generate_random_sdp(n, rand_seed=rand_seed)

# ------------------------------------------------------
# ------------------ Standard Form SDP -----------------
# ------------------------------------------------------
## Without decomposition
m3 = Model(optimizer)
sparse(@variable(m3, X[1:n, 1:n] in PSDCone()))
for i in 1:n
    @constraint(m3, sum(F[i] .* X) - c[i] == 0)
end
@objective(m3, Max, -sum(G .* X))
time_non_chordal = @elapsed optimize!(m3)

Xv = value.(X)
@info "termination status: $(termination_status(m3))"
@info "solution status: $(primal_status(m3))"
@info "Optimal value: $(round(objective_value(m3), digits=3))"

## Chordal Setup
m4 = Model(optimizer)
perm, iperm, Tls, cliques = get_selectors(sparsity_pattern([F..., G]); verbose=false, ret_cliques=true)
sp = CD.get_sparsity_pattern_from_cliques(cliques)
sp_inds = zip(findnz(sp)[1:2]...)
X = @variable(m4, X[i=1:n, j=1:n; (i,j) in sp_inds])

# TODO: Package this up into a function?
for i in 1:length(c)
    @constraint(m4, CD.tr(F[i][perm, perm], X) == c[i])
end

for (p, cl) in enumerate(cliques)
    len_cl = length(cl)
    Xp = @variable(m4, [1:len_cl, 1:len_cl], base_name="X$p")
    for i in 1:len_cl, j in 1:len_cl
        Xp[i,j] = X[cl[i], cl[j]]
    end
    @constraint(m4, Xp in PSDCone())
end

@objective(m4, Max, -CD.tr(G[perm, perm], X))

time_chordal = @elapsed optimize!(m4)
@info "termination status: $(termination_status(m4))"
@info "solution status: $(primal_status(m4))"
match = objective_value(m4) ≈ objective_value(m3)
@info "Optimal value: $(round(objective_value(m4), digits=3)), equals: $(match)"
@info "Time with chordal is $(round(time_chordal, digits=3))s vs $(round(time_non_chordal, digits=3))s"

X_ = CD.reconstruct_from_sparse_varref(X, n)
@assert maximum(X_ - X_') < 1e-8
X_ = 0.5*(X_ + X_')
Xcomp = maxdet_completion(X_)[iperm, iperm]
@assert sum(-G .* Xcomp) ≈ objective_value(m4)


text:
# Reading and parsing functions

#
# Reading
#

"""
    read(filename; warn=false) -> ::FDSNStationXML

Read a FDSN StationXML file with name `filename` from disk and return a
`FDSNStationXML` object.

If `warn` is `true`, then print warnings when attributes and elements are
encountered in the StationXML which are not expected.

Note that StationXML.jl reads `filename` first before then parsing the
string read in.  Do not `read` StationXML files which are larger than your
system's memory.
"""
function read(filename; warn=false)
    xml = EzXML.readxml(filename)
    _check_and_parse_xml(xml, filename, warn)
end

"""
    readstring(xml_string; filename=nothing, warn=false) -> ::FDSNStationXML

Read the FDSN StationXML contained in `xml_string` and return a `FDSNStationXML` object.

Optionally specify the `filename` from which the string was read.

If `warn` is `true`, then print warnings when attributes and elements are
encountered in the StationXML which are not expected.
"""
function readstring(xml_string; filename=nothing, warn=false)
    xml = EzXML.parsexml(xml_string)
    _check_and_parse_xml(xml, filename, warn)
end

function _check_and_parse_xml(xml, filename, warn)
    file_string = (filename === nothing || filename isa IO) ? "" : " in file $filename"
    xml_is_station_xml(xml) ||
        throw(ArgumentError("\"$file_string\" does not appear to be a StationXML file"))
    schema_version_is_okay(xml) ||
        throw(ArgumentError("StationXML$file_string does not have the correct schema version"))
    parse_node(xml.root, warn)
end

"""
    xml_is_station_xml(xml)

Return `true` if `xml` appears to be a StationXML file.
"""
xml_is_station_xml(xml) = EzXML.hasroot(xml) && xml.root.name == "FDSNStationXML"

"""
    schema_version_is_okay(xml::EzXML.Document) -> ::Bool

Return `true` if this XML document is of the correct version.

Currently fully supported is only v1.0 of the FDSN specification.
StationXML will parse v1.1 files, but ignore the additional information
provided in the new version.
"""
function schema_version_is_okay(xml::EzXML.Document)
    version = VersionNumber(xml.root["schemaVersion"])
    if version < v"1.2"
        return true
    elseif v"1.2.0" <= version < v"2"
        @warn("document is StationXML version $version; only v1.1 data will be read")
        return true
    else
        return false
    end
end

"""
    parse_node(root::EzXML.Node, warn=false) -> ::FDSNStationXML

Parse the `root` node of a StationXML document.  This can be accessed as
`EzXML.readxml(file).root`.

If `warn` is `true`, then warn about unexpected items in the StationXML.
"""
parse_node(root::EzXML.Node, warn::Bool=false) = parse_node(FDSNStationXML, root, warn)

"Types which can be directly parsed from a Node"
const ParsableTypes = Union{String, Float64, Int, DateTime}

function parse_node(::Type{T}, node::EzXML.Node, warn::Bool=false) where {T<:ParsableTypes}
    @debug("Parsing $T from \"$(node.content)\"")
    local_parse(T, node.content)
end

"""
    parse_node(::Type{Union{Missing,T}}, node::EzXML.Node, warn=false) where T -> ::T

Parse an XML node into a union of a type `T` with `Missing`, for optional
fields.
"""
function parse_node(::Type{Union{Missing,T}}, node::EzXML.Node, warn::Bool=false) where T
    @debug("Parsing a Union{Missing,$T} from \"$(node.content)\"")
    parse_node(T, node, warn)
end

"""
    parse_node(::Type{T}, node::EzXML.Node, warn::Bool=false) where T -> ::T

Parse the contents of `node` and return a type `T` containing the
information held within it.
"""
parse_node

@generated function parse_node(::Type{T}, node::EzXML.Node, warn::Bool=false) where T
    # Lists of fields
    attributes = attribute_fields(T)
    elements = collect(element_fields(T))
    fields = fieldnames(T)
    field_types = fieldtype.(Ref(T), fields)
    attribute_names = xml_attribute_name.(attributes)
    attribute_types = fieldtype.(Ref(T), attributes)
    element_names = xml_element_name.(elements)
    element_types = fieldtype.(Ref(T), elements)

    quote
        # Define function-local variables to pass in to constructor
        $([typ <: AbstractArray ?
            :($f = $(typ)()) :
            typ isa Union ?
                :(local $f = missing) :
                :(local $f)
            for (f, typ) in zip(fields, field_types)]...)

        # Attributes (always scalars)
        for att in EzXML.eachattribute(node)
            $([:(
                if att.name === $att_name
                    $var_name = local_parse($ftype, att.content)
                    continue
                end
                ) for (att_name, var_name, ftype) in
                    zip(attribute_names, attributes, attribute_types)]...)
        end

        # Elements
        for elm in EzXML.eachelement(node)
            # Scalars
            $([:(
                if elm.name === $elm_name
                    $var_name = parse_node($ftype, elm, warn)
                    continue
                end
                ) for (elm_name, var_name, ftype) in
                    zip(element_names, elements, element_types)
                    if !(ftype <: AbstractArray)]...)
            # Vectors
            $([:(
                if elm.name === $elm_name
                    push!($var_name, parse_node($(eltype(ftype)), elm, warn))
                    continue
                end
                ) for (elm_name, var_name, ftype) in
                    zip(element_names, elements, element_types)
                    if ftype <: AbstractArray]...)
        end

        # Pass variables to inner constructor
        $(T)($([:($f) for f in fields]...))
    end
end

# Version of parse which accepts String as the type.
# Don't define this for Base as this is type piracy.
local_tryparse(::Type{T}, s::AbstractString) where T<:AbstractString = s
local_tryparse(::Type{T}, s::AbstractString) where T = tryparse(T, s)
local_parse(::Type{<:AbstractString}, s::AbstractString) = s
local_parse(::Type{T}, s::AbstractString) where T = parse(T, s)
local_parse(::Type{Union{Missing, T}}, s::AbstractString) where T = local_parse(T, s)

# Handle dates with greater than millisecond precision by truncating to nearest millisecond,
# cope with UTC time zone information (ends with 'Z'), and convert non-UTC time zones to UTC
function local_parse(::Type{DateTime}, s::AbstractString)
    # Remove sub-millisecond intervals
    m = match(r"(.*T..:..:..[\.]?)([0-9]{0,3})[0-9]*([-+Z].*)*", s)
    m === nothing && throw(ArgumentError("invalid date-time string \"$(s)\""))
    dt = DateTime(m.captures[1] * m.captures[2]) # Local date to ms
    (m.captures[3] === nothing || m.captures[3] in ("", "Z", "+00:00", "-00:00")) && return dt # UTC
    pm = m.captures[3][1] # Whether ahead or behind UTC
    offset = Time(m.captures[3][2:end]) - Time("00:00")
    dt = pm == '+' ? dt + offset : dt - offset
    dt
end

#
# Writing
#

"""
    write(io, sxml::FDSNStationXML)

Write a `FDSNStationXML` structure to `io`, which may be a filename
or a `Base.IO` type.

# Example
(Note that `"example_stationxml_file.xml"` may not exist.)
```
julia> sxml = StationXML.read("example_stationxml_file.xml");

julia> write("new_file.xml", sxml)
```
"""
Base.write(io::IO, sxml::FDSNStationXML; warn=false) = EzXML.prettyprint(io, xmldoc(sxml; warn=warn))

"""
    xmldoc(sxml::FDSNStationXML; warn=false) -> xml::EzXML.XMLDocument

Create an XML document from `sxml`, a set of events of type `EventParameters`.
`xml` is an `EzXML.XMLDocument` suitable for output.

The StationXML document `xml` may be written with `write(io, xml)`
or converted to a string with `string(xml)`.

If `warn` is `true`, then warn about the presence of fields which have been
removed from the most recent StationXML specification, and hence will not be written.
"""
function xmldoc(sxml::FDSNStationXML; warn=false)
    version = "1"
    doc = EzXML.XMLDocument("1.0")
    root = EzXML.ElementNode("FDSNStationXML")
    # FIXME: Is this the only way to set a namespace in EzXML?
    namespace = EzXML.AttributeNode("xmlns", "http://www.fdsn.org/xml/station/" * version)
    EzXML.link!(root, namespace)
    EzXML.setroot!(doc, root)
    add_attributes!(root, sxml, warn)
    add_elements!(root, :fdsn_station_xml, sxml, warn)
    doc
end

"""
    add_attributes!(node, value, warn) -> node

Add the attribute fields from the structure `value` to a `node`.

If `warn` is `true`, then  warn about the presence of attributes which have
been removed from the most recent StationXML specification, and hence will not be added.
"""
function add_attributes!(node, value::T, warn) where T
    for field in attribute_fields(T)
        # Skip fields read from v1.0 but removed from v1.1
        if has_removed_fields(T) && field in removed_fields(T)
            warn && _warn_removed_field(field)
            continue
        # Write out StationXML v1.1 regardless of what was read in
        elseif field === :schema_version
            add_attribute!(node, field, DEFAULT_SCHEMA_VERSION)
            continue
        end
        content = getfield(value, field)
        content === missing && continue
        add_attribute!(node, field, content)
    end
    node
end

function add_attribute!(node, field, value)
    value === missing && return node
    name = xml_attribute_name(field)
    attr = EzXML.AttributeNode(name, string(value))
    EzXML.link!(node, attr)
    node
end

add_attribute!(node, field, value::EnumeratedStruct) = add_attribute!(node, field, value.value)

"""
    add_elements!(node, parent_field, value, warn) -> node

Add the elements to `node` contained within `value`.  `parent_field`
is the name of the field which contains `value`.  If `warn` is `true`, then
warn about the presence of elements which have been removed from the
most recent StationXML specification, and hence will not be added.
"""
function add_elements!(node, parent_field, value::T, warn) where T
    for field in element_fields(T)
        @debug("adding $parent_field: $field")
        # Skip fields read from v1.0 but removed from v1.1
        if has_removed_fields(T) && field in removed_fields(T)
            warn && _warn_removed_field(field)
            continue
        end
        content = getfield(value, field)
        if content === missing
            continue
        end
        add_element!(node, field, content, warn)
    end
    node
end

function add_elements!(node, parent_field, values::AbstractArray, warn)
    for value in values
        add_elements!(node, parent_field, value, warn)
    end
    node
end

"Union of types which can be natively written"
const WritableTypes = Union{Float64, Int, String, DateTime, Bool}

"""
    add_element!(node, field, value, warn) -> node

Add an element called `field` to `node` with content `value`.

`warn` is passed on to [`add_attributes!`](@ref) and [`add_elements!`](@ref).
"""
function add_element!(node, field, value::WritableTypes, warn)
    @debug("  adding writable type name \"$field\" with value \"$value\"")
    name = xml_element_name(field)
    elem = EzXML.ElementNode(name)
    add_text!(elem, value)
    EzXML.link!(node, elem)
    node
end

function add_element!(node, field, value::EnumeratedStruct, warn)
    @debug("  adding enumerated struct name \"$field\" with value \"$value\"")
    add_element!(node, field, value.value, warn)
end

function add_element!(node, field, values::AbstractArray, warn)
    @debug("  adding array type name \"$field\" with $(length(values)) values")
    for value in values
        add_element!(node, field, value, warn)
    end
    node
end

# Fallback for structs
function add_element!(node, field, value::T, warn) where T
    @debug("  adding compound type name \"$field\" of type \"$(typeof(value))\"")
    name = xml_element_name(field)
    elem = EzXML.ElementNode(name)
    EzXML.link!(node, elem)
    add_attributes!(elem, value, warn)
    add_elements!(elem, field, value, warn)
    if has_text_field(T)
        value_field = text_field(T)
        add_text!(elem, getfield(value, value_field))
    end
    node
end

"Add simple text to an element node"
function add_text!(node, value::WritableTypes)
    @debug("  adding simple text with value \"$value\"")
    content = EzXML.TextNode(string(value))
    EzXML.link!(node, content)
    node
end

"""
Issue a warning about an attribute or element `field` which is no longer
present in the most recent StationXML standard and hence will not be written.
"""
_warn_removed_field(field) =
    @warn("Not writing field $field: removed in StationXML v$DEFAULT_SCHEMA_VERSION")


text:
"""
`Basics.interpolate()`  ... interpolates a one- or two-dimensional function by different numerical methods, either by value or on a whole grid.

  + `("function: to new radial grid, Grasp92", (F::Array{Float64,1}, oldGrid::Radial.Grid), (G::Void, newGrid::Radial.Grid) )`
    ... to interpolate the (radial) function F from the oldgrid to newgrid, by a call to the Grasp92 Fortran procedures;
        an function G::Array{Float64,1} is returned. **Not yet implemented !**

  + `("function: to new radial grid, trapez rule", (F::Array{Float64,1}, oldGrid::Radial.Grid), (G::Void, newGrid::Radial.Grid)`
    ... to integrate the (radial) function F from the oldgrid to newgrid, by using a trapez rule from the Numerical Recipies.
        A function G::Array{Float64,1} is returned. **Not yet implemented !**
"""
function Basics.interpolate(sa::String, from::Tuple{Array{Float64,1},JAC.Radial.Grid}, to::Tuple{DataType,JAC.Radial.Grid})
    if       sa == "function: to new radial grid, Grasp92"          wa = interpolateOnGridGrasp92(from, to)
    elseif   sa == "function: to new radial grid, trapez rule"      wa = interpolateOnGridTrapezRule(from, to)
    else     error("Unsupported keystring = $sa.")
    end

end

"""
`Basics.interpolateOnGridGrasp92()`  ... interpolates ... by using Grasp92 Fortran procedures; see Basic.interpolate().
"""
function Basics.interpolateOnGridGrasp92(from::Tuple{Array{Float64,1},JAC.Radial.Grid}, to::Tuple{DataType,JAC.Radial.Grid})
    # First prepare all fields and settings to call quad_grasp92()
    error("Not yet implemented !")
end

"""
`Basics.interpolateOnGridTrapezRule()`  ... interpolates ... by using trapez rule from Numerical Recipies; see Basic.interpolate().
"""
function Basics.interpolateOnGridTrapezRule(from::Tuple{Array{Float64,1},JAC.Radial.Grid}, to::Tuple{DataType,JAC.Radial.Grid})
    error("Not yet implemented !")
end


text:
function binner(edges::AbstractRange)
    start = first(edges)
    stop = last(edges)
    dx = step(edges)
    maxi = length(edges)
    function tobin(x)::Int64
        if x < start
            1
        elseif x > stop
            maxi
        else
            trunc(Int, (x - start)/dx) + 1
        end
    end
    tobin
end

function histogram(X, Y, xedges, yedges)
    tobinx = binner(xedges)
    tobiny = binner(yedges)
    weights = zeros((length(xedges), length(yedges)))
    for (x, y) in zip(X, Y)
        i = tobinx(x)
        j = tobiny(y)
        weights[i, j] += 1
    end
    weights
end

function histogram(X::AbstractVector, edges::AbstractRange)::Vector{Int}
    tobin = binner(edges)
    weights = zeros(length(edges))
    for x in X
        i = tobin(x)
        weights[i] += 1
    end
    weights
end

function histogram_edges(X::AbstractVector; nbins=1000)
    edges = range(minimum(X), stop=maximum(X), length=nbins)
end

function histogram_edges(X, Y; nbins=1000)
    xedges = range(minimum(X), stop=maximum(X), length=nbins)
    yedges = range(minimum(Y), stop=maximum(Y), length=nbins)
    xedges, yedges
end

histogram_edges(X::AbstractMatrix; nbins=1000) = histogram_edges(X[:, 1], X[:, 2], nbins=nbins)
function histogram(X::AbstractMatrix; nbins=1000)
    edges = histogram_edges(X, nbins=nbins)
    histogram(X[:, 1], X[:, 2], edges...), edges
end

function histogram(X::AbstractVector; nbins=1000)
    edges = histogram_edges(X, nbins=nbins)
    histogram(X, edges), edges
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
module JLLPlots

import Plots;
import Unitful;

import ..JLLArrays

export plotas;

"""
    plotas(y; yunit=missing, kwargs...)
    plotas(x, y; unit=missing, xunit=missing, yunit=missing, xlabel="", ylabel="", kwargs...)

Plots quantities y or x and y, which should be arrays of Unitful.Quantities, in
either their native units if `unit`, `xunit`, and `yunit` are not given, or converted
to those units before plotting. `unit` is applied to both x and y, `xunit` and
`yunit` override it for x and y, respectively. `unit`, `xunit`, and `yunit` must
be Unitful.Units values. For example:

    using Unitful;
    using JLLUtils;

    x = (1:10)u"m";
    y = (1000:1000:10000)u"m";

    plotas(x,y) # plots x and y in meters, their native unit
    plotas(x,y; unit=u"km") # plots both x and y in kilometers
    plotas(x,y; unit=u"km", xunit=u"m") # plots x in meters, y in kilometers
    plotas(x,y; yunit=u"km") # also plots x in meters, y in kilometers

The unit will be given in parentheses in the axis label. If you pass a string to
the `xlabel` or `ylabel` keywords, that string will be printed before the unit.
For example,

    plotas(x,y; yunit=u"km", xlabel="Distance", ylabel="Altitude")

would print "Distance (m)" as the x-axis label and "Altitude (km)" as the y-axis
label.

If only y is given, then x is created as a sequential array the same size as y.
Any extra keyword arguments are passed through to Plots.plot.

NOTE: the units of x and y are not saved anywhere, so if you later use `plot!`
to add more series, it is your resposibility to make sure that the new series
are the same units.
"""
function plotas(y; yunit=missing, kwargs...)
    x = JLLArrays.seq_mat(size(y)...);
    plotas(x, y; yunit=yunit, kwargs...)
end

function plotas(x, y; unit=missing, xunit=missing, yunit=missing, xlabel="", ylabel="", kwargs...)
    x, y = _convert_xy(x, y; unit=unit, xunit=xunit, yunit=yunit);

    xlabel *= " ($(Unitful.unit(x[1])))";
    ylabel *= " ($(Unitful.unit(y[1])))";

    return Plots.plot(Unitful.ustrip(x), Unitful.ustrip(y); xlabel=xlabel, ylabel=ylabel, kwargs...);
end

function plotas!(x, y; unit=missing, xunit=missing, yunit=missing, kwargs...)
    x, y = _convert_xy(x, y; unit=unit, xunit=xunit, yunit=yunit);
    return Plots.plot!(Unitful.ustrip(x), Unitful.ustrip(y); kwargs...)
end

function _convert_xy(x, y; unit, xunit, yunit)
    xunit = ismissing(xunit) ? unit : xunit;
    yunit = ismissing(yunit) ? unit : yunit;

    x = !ismissing(xunit) ? Unitful.uconvert.(xunit, x) : x;
    y = !ismissing(yunit) ? Unitful.uconvert.(yunit, y) : y;
    return x, y
end

end #module

max_stars_repo_name: joshua-laughner/Julia-JLLUtils

max_issues_repo_name: joshua-laughner/Julia-JLLUtils

max_forks_repo_name: joshua-laughner/Julia-JLLUtils