text:
handle("gotosymbol") do data
  @destruct [
    word,
    path || nothing,
    # local context
    column || 1,
    row || 1,
    startRow || 0,
    context || "",
    onlyGlobal || true,
    # module context
    mod || "Main",
    text || "",
  ] = data
  gotosymbol(
    word, path,
    column, row, startRow, context, onlyGlobal,
    mod, text,
  )
end

function gotosymbol(
  word, path = nothing,
  # local context
  column = 1, row = 1, startrow = 0, context = "", onlyglobal = false,
  # module context
  mod = "Main", text = ""
)
  try
    # local goto
    if !onlyglobal
      localitems = localgotoitem(word, path, column, row, startrow, context)
      isempty(localitems) || return Dict(
        :error => false,
        :items => todict.(localitems)
      )
    end

    # global goto
    m = getmodule(mod)
    globalitems = if m == Main && mod ≠ "Main" # unloaded module
      globalgotoitems_unloaded(word, mod)
    else
      globalgotoitems(word, m, path, text)
    end
    isempty(globalitems) || return Dict(
      :error => false,
      :items => todict.(globalitems),
    )
  catch err
    return Dict(:error => true)
  end

  return Dict(:error => true) # nothing hits
end

struct GotoItem
  name::String
  text::String
  file::String
  line::Int

  GotoItem(name::String, text::String, file::String, line::Int = 0) =
    new(name, text, normpath(file), line)
end
GotoItem(name::String, file::String, line::Int = 0) = GotoItem(name, name, file, line)

# for messaging over julia ⟷ Atom
todict(gotoitem::GotoItem) = Dict(
  :text      => gotoitem.text,
  :file      => gotoitem.file,
  :line      => gotoitem.line,
  :secondary => string(gotoitem.file, ':', gotoitem.line + 1),
)

### local goto

function localgotoitem(word, path, column, row, startrow, context)
  word = first(split(word, '.')) # always ignore dot accessors
  position = row - startrow
  ls = locals(context, position, column)
  filter!(l -> l.name == word && l.line < position, ls)
  return map(ls) do l # there should be zero or one element in `ls`
    name = l.name
    line = startrow + l.line - 1
    GotoItem(name, path, line)
  end
end
localgotoitem(word, ::Nothing, column, row, startrow, context) = GotoItem[] # when called from docpane/workspace

### global goto - bundles toplevel gotos & method gotos

# goto for loaded module -- using runtime information in need
function globalgotoitems(word, mod, path, text)
  # strip a dot-accessed module if exists
  identifiers = split(word, '.')
  head = string(identifiers[1])
  if head ≠ word && (nextmod = getfield′(mod, head)) isa Module
    # if `head` is a module, update `word` and `mod`
    nextword = join(identifiers[2:end], '.')
    return globalgotoitems(nextword, nextmod, text, path)
  end

  val = getfield′(mod, word)
  val isa Module && return [GotoItem(val)] # module goto

  items = toplevelgotoitems(word, mod, path, text)

  # append method gotos that are not caught by `toplevelgotoitems`
  ml = methods(val)
  files = map(item -> item.file, items)
  methoditems = filter!(item -> item.file ∉ files, methodgotoitems(ml))
  append!(items, methoditems)
  return items
end
# goto for unloaded package
function globalgotoitems_unloaded(word, pkg::String)
  pathitemsmap = if haskey(SYMBOLSCACHE, pkg)
    SYMBOLSCACHE[pkg]
  else
    (path = Base.find_package(pkg)) === nothing && return GotoItem[]
    SYMBOLSCACHE[pkg] = _collecttoplevelitems_static(pkg, path)
  end

  ret = GotoItem[]
  for (_, items) in pathitemsmap
    @>> filter(let name = word
      item -> name == item.name
    end, items) append!(ret)
  end
  return ret
end

## module goto

function GotoItem(mod::Module)
  name = string(mod)
  file, line = mod == Main ? MAIN_MODULE_LOCATION[] : moduledefinition(mod)
  GotoItem(name, file, line - 1)
end

## toplevel goto

const PathItemsMap = Dict{String, Vector{GotoItem}}

# NOTE:
# this data structure should only keep "simple" data so that it doesn't eat
# too much memory when the cache has got fat (like via `regeneratesymbols`)
# TODO?:
# maybe we will end up creating a more generic struct that represents a "symbol",
# and wanting `SYMBOLSCACHE` to store them instead of `GotoItem`
# (should be kept "simple" atst though, and ideally it handles both toplevel and local symbols)
# so that we can use `SYMBOLSCACHE` as a general cache interface for various static features
# like goto, find references, and so on.

"""
    Atom.SYMBOLSCACHE

"module" (`String`) ⟶ "path" (`String`) ⟶ "symbols" (`Vector{GotoItem}`) map.

!!! note
    "module" should be canonical, i.e.: should be identical to names that are
      constructed from `string(mod::Module)`.
"""
const SYMBOLSCACHE = Dict{String, PathItemsMap}()

function toplevelgotoitems(word, mod, path, text)
  key = string(mod)
  pathitemsmap = if haskey(SYMBOLSCACHE, key)
    SYMBOLSCACHE[key]
  else
    SYMBOLSCACHE[key] = collecttoplevelitems(mod, path, text) # caching
  end

  ret = GotoItem[]
  for (_, items) in pathitemsmap
    @>> filter(let name = word
      item -> name == item.name
    end, items) append!(ret)
  end
  return ret
end

# entry methods
function collecttoplevelitems(mod::Module, path::String, text::String)
  return if mod == Main || isuntitled(path)
    # for `Main` module and unsaved editors, always use CSTPraser-based approach
    # with a given buffer text, and don't check module validity
    _collecttoplevelitems_static(nothing, path, text)
  else
    _collecttoplevelitems(mod)
  end
end
# when `path === nothing`, e.g.: called from docpane/workspace
collecttoplevelitems(mod::Module, path::Nothing, text::String) = _collecttoplevelitems(mod)

function _collecttoplevelitems(mod::Module)
  entrypath, paths = modulefiles(mod)
  return if entrypath !== nothing # Revise-like approach
    _collecttoplevelitems_loaded(stripdotprefixes(string(mod)), [entrypath; paths])
  else # if Revise-like approach fails, fallback to CSTParser-based approach
    entrypath, line = moduledefinition(mod)
    _collecttoplevelitems_static(stripdotprefixes(string(mod)), entrypath)
  end
end

# module-traverse for given files, which are collected by Revise-like approach:
# NOTE: only works for loaded precompiled modules
function _collecttoplevelitems_loaded(mod::Union{Nothing, String}, paths::Vector{String})
  pathitemsmap = PathItemsMap()

  entrypath, paths = paths[1], paths[2:end]

  # ignore toplevel items outside of `mod`
  if isfile′(entrypath)
    items = toplevelitems(read(entrypath, String); mod = mod)
    pathitemsmap[entrypath] = GotoItem.(entrypath, items)
  end

  # collect symbols in included files (always in `mod`)
  for path in paths
    isfile′(path) || continue
    items = toplevelitems(read(path, String); mod = mod, inmod = true)
    pathitemsmap[path] = GotoItem.(path, items)
  end

  return pathitemsmap
end

# module-traverse based on CSTParser, looking for toplevel `included` calls
function _collecttoplevelitems_static(mod::Union{Nothing, String}, entrypath::String, pathitemsmap::PathItemsMap = PathItemsMap(); inmod = false)
  isfile′(entrypath) || return pathitemsmap
  # escape recursive `include` loops
  entrypath in keys(pathitemsmap) && return pathitemsmap
  text = read(entrypath, String)
  _collecttoplevelitems_static(mod, entrypath, text, pathitemsmap; inmod = inmod)
end
function _collecttoplevelitems_static(mod::Union{Nothing, String}, entrypath::String, text::String, pathitemsmap::PathItemsMap = PathItemsMap(); inmod = false)
  items = toplevelitems(text; mod = mod, inmod = inmod)
  pathitemsmap[entrypath] = GotoItem.(entrypath, items)

  # looking for toplevel `include` calls
  for item in items
    if item isa ToplevelCall
      expr = item.expr
      if isinclude(expr)
        nextfile = expr.args[3].val
        nextentrypath = joinpath(dirname(entrypath), nextfile)
        isfile′(nextentrypath) || continue
        # `nextentrypath` is always in `mod`
        _collecttoplevelitems_static(mod, nextentrypath, pathitemsmap; inmod = true)
      end
    end
  end

  return pathitemsmap
end

GotoItem(path::String, item::ToplevelItem) = GotoItem("", path) # fallback case
function GotoItem(path::String, binding::ToplevelBinding)
  expr = binding.expr
  bind = binding.bind
  name = CSTParser.defines_macro(expr) ? string('@', bind.name) : bind.name
  text = CSTParser.has_sig(expr) ? str_value(CSTParser.get_sig(expr)) : name
  line = binding.lines.start - 1
  GotoItem(name, text, path, line)
end

## update toplevel symbols cache
# NOTE:
# - only updates toplevel symbols in the current text buffer
# - handled by the `updateeditor` handler in outline.jl

function updatesymbols(mod, path::Nothing, text) end # fallback case
function updatesymbols(mod, path::String, text)
  m = getmodule(mod)
  mod = string(m)

  # initialize the cache if there is no previous one
  if !haskey(SYMBOLSCACHE, mod)
    SYMBOLSCACHE[mod] = collecttoplevelitems(m, path, text)
  end

  # ignore toplevel items outside of `mod` when `path` is an entry file
  entrypath, _ = moduledefinition(m)
  inmod = path ≠ entrypath
  items = toplevelitems(text; mod = stripdotprefixes(mod), inmod = inmod)
  push!(SYMBOLSCACHE[mod], path => GotoItem.(path, items))
end

## generate toplevel symbols cache

handle("regeneratesymbols") do
  with_logger(JunoProgressLogger()) do
    regeneratesymbols()
  end
  nothing
end

function regeneratesymbols()
  id = "regenerate_symbols_progress"
  @info "Generating symbols cache in loaded modules" progress=0 _id=id

  loaded = Set(string.(Base.loaded_modules_array()))
  pkgs = if isdefined(Pkg, :dependencies)
    getfield.(values(Pkg.dependencies()), :name)
  else
    collect(keys(Pkg.installed()))
  end
  unloaded = filter(pkg -> pkg ∉ loaded, pkgs)
  loadedlen = length(loaded)
  unloadedlen = length(unloaded)
  total = loadedlen + unloadedlen

  for (i, mod) in enumerate(Base.loaded_modules_array())
    try
      key = string(mod)
      key == "__PackagePrecompilationStatementModule" && continue # will cause error

      @logmsg -1 "Symbols: $key ($i / $total)" progress=i/total _id=id
      SYMBOLSCACHE[key] = _collecttoplevelitems(mod)
    catch err
      @error err
    end
  end

  for (i, pkg) in enumerate(unloaded)
    try
      @logmsg -1 "Symbols: $pkg ($(i + loadedlen) / $total)" progress=(i+loadedlen)/total _id=id
      (path = Base.find_package(pkg)) === nothing && continue
      SYMBOLSCACHE[pkg] = _collecttoplevelitems_static(pkg, path)
    catch err
      @error err
    end
  end

  @info "Finished generating the symbols cache" progress=1 _id=id
end

## clear toplevel symbols cache

handle("clearsymbols") do
  clearsymbols()
  nothing
end

function clearsymbols()
  for key in keys(SYMBOLSCACHE)
    delete!(SYMBOLSCACHE, key)
  end
end

## method goto

methodgotoitems(ml) = GotoItem.(aggregatemethods(ml))

# aggregate methods with default arguments to the ones with full arguments
function aggregatemethods(ml)
  ms = collect(ml)
  sort!(ms, by = m -> m.nargs, rev = true)
  return @static VERSION ≥ v"1.1" ? unique!(m -> (m.file, m.line), ms) : unique(m -> (m.file, m.line), ms)
end

function GotoItem(m::Method)
  text = replace(sprint(show, m), methodloc_regex => s"\g<sig>")
  _, link = view(m)
  file = link.file
  line = link.line - 1
  GotoItem(text, file, line)
end


text:
#AbstractTensor
TO.checked_similar_from_indices(C, T::Type, p1, p2, A::AbstractTensor, CA::Symbol = :N) =
    TO.checked_similar_from_indices(C, T, (p1..., p2...), A, CA)

TO.checked_similar_from_indices(C, T::Type, poA, poB, p1, p2, A::AbstractTensor,
        B::AbstractTensor, CA::Symbol = :N, CB::Symbol = :N) =
    TO.checked_similar_from_indices(C, T, poA, poB, (p1..., p2...), A, B, CA, CB)

TO.add!(α, A::AbstractTensor, conjA, β, C::AbstractTensor, p1, p2)=
    TO.add!(α, A, conjA, β, C, (p1..., p2...))

TO.trace!(α, A::AbstractTensor, CA, β, C::AbstractTensor, indleft, indright, cind1, cind2) =
    TO.trace!(α, A, CA, β, C, (indleft..., indright...), cind1, cind2)

TO.contract!(α, A::AbstractTensor, CA, B::AbstractTensor, CB, β, C::AbstractTensor, oindA, cindA, oindB,
        cindB, indleft, indright, syms = nothing) =
    TO.contract!(α, A, CA, B, CB, β, C, oindA, cindA, oindB, cindB, (indleft..., indright...), syms)

#DTensor
TO.scalar(A::DTensor)  = TO.scalar(A.array)

function TO.checked_similar_from_indices(C, ::Type{T}, ind, A::DTensor, CA::Symbol) where T
    sz = map(n->size(A, n), ind)
    if C !== nothing && C isa DTensor && sz == size(C) && T == eltype(C)
        CT = DTensor{T,length(sz)}
        return C::CT
    else
        return similar(A, T, sz)
    end
end

function TO.checked_similar_from_indices(C, ::Type{T}, poA, poB, ind, A::DTensor,
        B::DTensor, CA::Symbol, CB::Symbol) where {T,N}
    oszA = map(n->size(A,n), poA)
    oszB = map(n->size(B,n), poB)
    sz = let osz = (oszA..., oszB...)
        map(n->osz[n], ind)
    end
    if C !== nothing && C isa DTensor && sz == size(C) && T == eltype(C)
        CT = DTensor{T,length(sz)}
        return C::CT
    else
        return similar(A, T, sz)
    end
end

TO.add!(α, A::DTensor, conjA, β, C::DTensor, indCinA) =
    DTensor(TO.add!(α, A.array, conjA, β, C.array, indCinA))

TO.trace!(α, A::DTensor, CA, β, C::DTensor, indCinA, cindA1, cindA2) =
    DTensor(TO.trace!(α, A.array, CA, β, C.array, indCinA, cindA1, cindA2))

TO.contract!(α, A::DTensor, CA::Symbol, B::DTensor, CB::Symbol, β, C::DTensor,
        oindA, cindA, oindB, cindB, indCinoAB, syms::Union{Nothing, NTuple{3,Symbol}} = nothing) =
        DTensor(TO.contract!(α, A.array, CA, B.array, CB, β, C.array,
        oindA, cindA, oindB, cindB, indCinoAB, syms))

#DASTensor
TO.scalar(a::DASTensor) = TO.scalar(first(values(a)))

function TO.checked_similar_from_indices(C, ::Type{T}, ind,
        A::DASTensor{TA,NA,SYM,CHS,SS,CH}, CA::Symbol) where {T,TA,NA,SYM,CHS,SS,CH}
    if CA == :N
        dims = sizes(A,ind)
        ios  = in_out(A,ind)
    else
        dims = reverse.(sizes(A,ind))
        ios  = inv(in_out(A,ind))
    end
    chs = charges(A,ind)

    if C !== nothing && C isa DASTensor && symmetry(C) == SYM && dims == sizes(C) &&
            ios == in_out(C) && chs == charges(C) && T == eltype(C) &&
            charge(C) == charge(A)
        CT = DASTensor{T,length(ind),SYM,CHS,SS,CH}
        return C::CT
    else
        return DASTensor{T,length(ind)}(SYM, chs, deepcopy(dims), ios, charge(A))
    end
end

function TO.checked_similar_from_indices(C, ::Type{T}, poA, poB, ind,
        A::DASTensor{TA,NA,SYM,CHS,SS,CH},
        B::DASTensor{TB,NB,SYM,CHS,SS,CH},
        CA::Symbol, CB::Symbol) where {T,TA,NA,TB,NB,SYM,CHS,SS,CH}
    chs = TT.getindices(TT.vcat(charges(A,poA), charges(B,poB)), ind)
    dims = TT.getindices(TT.vcat(
                TT.getindices(CA == :N ? sizes(A) : reverse.(sizes(A)), poA),
                TT.getindices(CB == :N ? sizes(B) : reverse.(sizes(B)), poB)),
                ind)
    dimsA =  CA == :N ? sizes(A,poA) : map(reverse, sizes(A,poA))
    dimsB =  CB == :N ? sizes(B,poB) : map(reverse, sizes(B,poB))
    dims = TT.getindices(TT.vcat(dimsA,dimsB), ind)
    ios = vcat(ifelse(CA == :N, in_out(A,poA), inv(in_out(A,poA))),
                ifelse(CB == :N, in_out(B,poB), inv(in_out(B,poB))))[ind]
    if C !== nothing && C isa DASTensor && symmetry(C) == SYM && dims == sizes(C) &&
            ios == in_out(C) && chs == charges(C) && T == eltype(C) &&
            charge(C) == charge(A) ⊕ charge(B) CT = DASTensor{T,length(ind),SYM,CHS,SS,CH}
        return C::CT
    else
        return DASTensor{T,length(ind)}(SYM,chs, deepcopy(dims), ios, charge(A) ⊕ charge(B))
    end
end

function _errorsadd(A::DASTensor{T,N,SYM}, C, indCinA::NTuple{M}) where {T,N,SYM,M}
    mask = in_out(A,indCinA).v .== in_out(C).v
    for (m, iA, iC) in zip(mask, indCinA, 1:M)
        sizes(A, iA) == ifelse(m,sizes(C, iC), reverse(sizes(C,iC))) ||
            throw(DimensionMismatch())
        charges(A, iA) == ifelse(m, charges(C, iC), inv(charges(C, iC))) ||
            throw(ArgumentError("charges don't agree"))
    end
    modio = InOut(ntuple(i -> ifelse(mask[i],1,-1),Val(N))...)
    return ⊗(modio)
end

function TO.add!(α::Number, A::DASTensor{T,N}, CA, β::Number, C::DASTensor{S,N},
            indCinA) where {T,S,N}
    maskfun = _errorsadd(A, C, indCinA)
    for (sector, degeneracy) in tensor(A)
        permsector = maskfun(sector)[indCinA]
        if haskey(tensor(C), permsector)
            TO.add!(α, degeneracy, CA, β, C[permsector], indCinA)
        else
            C[permsector] = TO.checked_similar_from_indices(nothing, T, indCinA, degeneracy, CA)
            TO.add!(α, degeneracy, CA, 0, C[permsector], indCinA)
        end
    end
    return C
end

function _errorstrace(A::DASTensor{T,N,<:Any},
                         cindA1::NTuple{M,Int},
                         cindA2::NTuple{M,Int},
                         C::DASTensor{TC,NC},
                         indCinA::NTuple{NC}) where {M,T,N,TC,NC}
    maskA = in_out(A, cindA1).v .== (inv(in_out(A, cindA2))).v
    for (m, iA1, iA2) in zip(maskA, cindA1, cindA2) #trace in A
        sizes(A, iA1) == ifelse(m, sizes(A, iA2), reverse(sizes(A,iA2))) ||
            throw(DimensionMismatch())
        charges(A, iA1) == ifelse(m, charges(A, iA2), inv(charges(A,iA2))) ||
            throw(ArgumentError("charges don't agree"))
    end

    maskC = in_out(A, indCinA).v .== in_out(C).v
    for (m, iA, iC) in zip(maskC, indCinA, 1:NC) #trace in A
        sizes(A, iA) == ifelse(m, sizes(C, iC), reverse(sizes(C,iC))) ||
            throw(DimensionMismatch())
        charges(A, iA) == ifelse(m, charges(C, iC), inv(charges(C,iC))) ||
            throw(ArgumentError("charges don't agree"))
    end
    maskAio = InOut(ntuple(i -> ifelse(maskA[i],1,-1), Val(M))...)
    maskCio = InOut(ntuple(i -> ifelse(maskC[i],1,-1), Val(NC))...)
    return (⊗(maskAio), ⊗(maskCio))
end

function TO.trace!(α, A::DASTensor{T,N}, CA, β, C::DASTensor{S,M},
                    indCinA, cindA1, cindA2) where {T,N,S,M}
    #conditions
    maskAfun, maskCfun = _errorstrace(A, cindA1, cindA2, C, indCinA)

    sectors = filter(x -> isequal(maskAfun(x[cindA1]), x[cindA2]), keys(A))
    cinds = TT.vcat(cindA1, cindA2)

    passedset = Set{DASSector{M,chargetype(symmetry(A))}}()
    for sector in sectors
        secC = maskCfun(sector[indCinA])
        if haskey(C, secC)
            if !in(secC, passedset)
                TO.trace!(α, A[sector], CA, β, C[secC], indCinA, cindA1, cindA2)
                push!(passedset, secC)
            else
                TO.trace!(α, A[sector], CA, 1, C[secC], indCinA, cindA1, cindA2)
            end
        else
            C[secC] = TO.checked_similar_from_indices(nothing, T, indCinA, A[sector], CA)
            TO.trace!(α, A[sector], CA, 0, C[secC], indCinA, cindA1, cindA2)
            push!(passedset, secC)
        end
    end
    return C
end

function _errorscontract(A::DASTensor{TA,NA,<:Any},
                        (oindA, cindA)::Tuple{NTuple{NoA,Int}, NTuple{CoA,Int}},
                        B::DASTensor{TB,NB,<:Any},
                        (oindB, cindB)::Tuple{NTuple{NoB,Int}, NTuple{CoB,Int}},
                        C::DASTensor{TC,NC,<:Any},
                        indCinoAB::NTuple{NoC,Int}) where
                            {NA, NB, NC, TA, TB, TC, NoA, CoA, NoB, CoB, NoC}
    CoA == CoB || throw(ArgumentError("indices to contract don't pair up"))
    NoA + NoB == NoC || throw(ArgumentError("indices to contract don't pair up"))

    maskB = in_out(A, cindA).v .== inv(in_out(B, cindB)).v
    for (iA, iB, m)  in zip(cindA, cindB, maskB)
        charges(A,iA) == ifelse(m, charges(B,iB), inv(charges(B,iB)))||
            throw(ArgumentError("charges don't agree"))
        sizes(A,iA) == ifelse(m, sizes(B,iB), reverse(sizes(B,iB))) ||
            throw(DimensionMismatch())
    end
    ioAB  = vcat(in_out(A, oindA), in_out(B, oindB))
    chsAB = TT.vcat(charges(A, oindA), charges(B, oindB))
    dsAB  = TT.vcat(sizes(A, oindA), sizes(B, oindB))

    maskAB = ioAB[indCinoAB].v .== in_out(C).v
    for (iC, iAB, m) in zip(1:NoC, indCinoAB, maskAB)
        chsAB[iAB] == ifelse(m, charges(C, iC), inv(charges(C,iC)))||
            throw(ArgumentError("charges don't agree"))
        dsAB[iAB] == ifelse(m, sizes(C,iC), reverse(sizes(C,iC))) ||
             throw(DimensionMismatch())
    end

    maskBio = InOut(tuple(ifelse.(maskB,1,-1)...)...)
    maskABio = InOut(tuple(ifelse.(maskAB,1,-1)...)...)
    return (⊗(maskBio), ⊗(maskABio))
end

function TO.contract!(α, A::DASTensor{TA,NA,SYM}, CA,
                      B::DASTensor{TB,NB,SYM}, CB, β,
                      C::DASTensor{TC,NC,SYM}, oindA, cindA, oindB, cindB,
                      indCinoAB, syms=nothing) where
                      {TA,NA,TB,NB,TC,NC,M,SYM}
    #conditions
    maskBfun, maskABfun = _errorscontract(A, (oindA, cindA), B, (oindB, cindB), C, indCinoAB)

    oinAB = TT.vcat(oindA, oindB .+ NA)
    indCinAB = TT.getindices(oinAB, indCinoAB)
    secsA = groupby(x -> x[cindA], keys(A))
    secsB = groupby(x -> maskBfun(x[cindB]), keys(B))
    # collect sectors that contract with each other
    secsAB = intersect(keys(secsA), keys(secsB))
    passedset = Set{DASSector{NC,chargetype(SYM)}}()
    for sector in secsAB
        for secA in secsA[sector], secB in secsB[sector]
            secC = maskABfun(permute(vcat(secA[oindA], secB[oindB]), indCinoAB))
            if haskey(C, secC)
                if in(secC, passedset) #firstpass
                    TO.contract!(α, A[secA], CA, B[secB], CB, 1, C[secC],
                              oindA, cindA, oindB, cindB, indCinoAB,syms)
                else
                    push!(passedset, secC)
                    TO.contract!(α, A[secA], CA, B[secB], CB, β, C[secC],
                              oindA, cindA, oindB, cindB, indCinoAB,syms)
                 end
             else
                C[secC] = TO.checked_similar_from_indices(nothing, TC,
                                    oindA, oindB, indCinoAB, (),
                                    A[secA], B[secB], :N, :N)
                TO.contract!(α, A[secA], CA, B[secB], CB, 0, C[secC],
                          oindA, cindA, oindB, cindB, indCinoAB,syms)
                push!(passedset, secC)
             end
         end
     end
     return C
 end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ ed9ab030-e61d-11eb-32f3-5be220868ae7
using Revise

# ╔═╡ d51f978f-a2c8-410a-8836-f57fc362eed1
using CSV

# ╔═╡ ffa81468-a1cd-44f3-aa5a-c970e56b7f27
using DataFrames

# ╔═╡ 4c98c769-0aed-4b8a-bf92-fb06268100a9
using Plots

# ╔═╡ c819bc8f-cc85-4edb-8289-5884f7fbb8f7
using LaTeXStrings

# ╔═╡ e2d229a0-eaae-481c-9ee1-94d641fc3014
using ColorSchemes

# ╔═╡ 88fe20b4-303a-494e-95cc-e8183b2f474e
Ω = real.(parse.(ComplexF64, (CSV.File("multi_conductivity_data.csv") |> Tables.matrix)[2:end, 1]))

# ╔═╡ a76c9858-0fe7-45c2-9b5c-a8cf7f800c08
T = real.(parse.(ComplexF64, (CSV.File("multi_conductivity_data.csv") |> Tables.matrix)[1, 2:end]))

# ╔═╡ 6d08ec24-51d5-4105-bddf-443bc8496456
A_freq = (CSV.File("B_data.csv") |> Tables.matrix)[:, 2]

# ╔═╡ 5fa3101f-d24b-41f1-9a77-f5b7782ec4a3
B_freq = (CSV.File("B_data.csv") |> Tables.matrix)[:, 2]

# ╔═╡ 78156817-9c43-44c0-90e9-55e1100538a7
multi_freqs = real.((CSV.File("multi_data.csv") |> Tables.matrix)[1:end, 2])

# ╔═╡ 4bf7a882-1306-415a-8760-6191e5284306
multi_α = real.((CSV.File("multi_data.csv") |> Tables.matrix)[1:end, 1])

# ╔═╡ 8b91053e-28cc-47ce-8551-72e3e3c8ee93
A_α = real.((CSV.File("A_data.csv") |> Tables.matrix)[1:end, 1])

# ╔═╡ fd24b1c2-a6c2-4328-854b-ada51d40c56a
B_α = real.((CSV.File("B_data.csv") |> Tables.matrix)[1:end, 1])

# ╔═╡ db273a6b-bd14-4935-92dd-1271def823ba
ir_activity = real.((CSV.File("multi_data.csv") |> Tables.matrix)[1:end, end])

# ╔═╡ 830ce1e4-80f4-4dfa-a214-d7d0c84cfbbb
begin
	scatter(multi_freqs, multi_α, yaxis = :log, xticks = (multi_freqs, hcat(["$(round(i, digits = 2))" for i in multi_freqs]...)), label = L"\textbf{\alpha_j}", size = (500, 500), xlabel = "Phonon Frequencies (THz)", tickfontsize = 9, xrotation = 90)
	scatter!(multi_freqs, ir_activity, markershape = :diamond, label = "IR")
	plot!(multi_freqs, multi_α ./ ir_activity, linewidth = 2, linestyle = :dash, label = "ratio")
end

# ╔═╡ 4b5dc310-b4f6-4bd3-88e5-f082ea8ee47e
F_multi = real.((CSV.File("multi_vw.csv") |> Tables.matrix)[1:end, 5])

# ╔═╡ f9e47e6f-fdc5-42b6-9ea2-46c801f8cf6b
F_multi_avg = [F_multi[i] for i in 1:length(T)]

# ╔═╡ e2a502b1-cca9-40c2-a107-2b577f0662dc
F_A = real.((CSV.File("A_data.csv") |> Tables.matrix)[1:end, end])

# ╔═╡ 92f361d6-2052-45bd-a2a1-0d19d04eec1b
F_B = real.((CSV.File("B_data.csv") |> Tables.matrix)[1:end, end])

# ╔═╡ b1473f00-5aa9-4ace-a0f9-9e5cb9747af5
multi_energy_temp = plot(T, F_multi, label = hcat(["$(round(i, digits = 3)) THz" for i in multi_freqs]...), legend = :outerright, linewidth = 2, size = (650, 500), minorgrid = true, ylabel = "Free energy (meV)", xlabel = "T (K)", tickfontsize = 9, legendfontsize = 9)

# ╔═╡ 4c4d1b61-6710-478c-a8f5-067c219be342
savefig(multi_energy_temp, "multi_energy_temp.pdf")

# ╔═╡ 9493c1f8-d2b4-42b3-9af5-2d69ee7499bf
begin
	free_energy_temp = plot(T, F_A, linewidth = 2, linestyle = :dashdot, label = "H-A", tickfontsize = 12, labelfontsize = 12, legendfontsize = 12, xlabel = "T(K)", ylabel = "Free energy (meV)", minorgrid = true, size = (550, 550), legend = :topright)
	plot!(T, F_B, linewidth = 2, linestyle = :dash, label = "H-B")
	plot!(T, F_multi_avg, linewidth = 2, linestyle = :solid, label = "multi")
end

# ╔═╡ d4a13a12-e1ca-4ded-9c49-3c33b1cb69a0
savefig(free_energy_temp, "free_energy_temp.pdf")

# ╔═╡ 40fc4b41-c835-41a3-b468-765d4bbd0042
begin
	v_A = real.((CSV.File("A_data.csv") |> Tables.matrix)[1:end, end-2])
	w_A = real.((CSV.File("A_data.csv") |> Tables.matrix)[1:end, end-1])
	v_B = real.((CSV.File("B_data.csv") |> Tables.matrix)[1:end, end-2])
	w_B = real.((CSV.File("B_data.csv") |> Tables.matrix)[1:end, end-1])
end

# ╔═╡ aa1b7720-91ee-48e5-a97e-f8e3a499ad2f
v_multi = real.((CSV.File("multi_vw.csv") |> Tables.matrix)[1:end, 3])

# ╔═╡ 10b31119-9d2f-445d-9353-9e6788a464d3
w_multi = real.((CSV.File("multi_vw.csv") |> Tables.matrix)[1:end, 4])

# ╔═╡ eef00978-0153-4031-b748-d9e1e5e2f5c5
multi_v_temp = plot(T, v_multi .* (multi_freqs)', label = hcat(["$(round(i, digits = 3)) THz" for i in multi_freqs]...), legend = :outerright, linewidth = 1.5, size = (650, 500), minorgrid = true, ylabel = "v (THz)", xlabel = "T (K)", tickfontsize = 9, legendfontsize = 9)

# ╔═╡ 45e2101a-010d-404d-aaf6-c68d28536677
multi_w_temp = plot(T, w_multi .* multi_freqs', label = hcat(["$(round(i, digits = 3)) THz" for i in multi_freqs]...), legend = :outerright, linewidth = 1.5, size = (650, 500), minorgrid = true, ylabel = "w (THz)", xlabel = "T (K)", tickfontsize = 9, legendfontsize = 9)

# ╔═╡ 4da848c3-23dd-4640-bab3-a8b4b4784639
plot(T, (v_multi .- w_multi) .* multi_freqs', label = hcat(["$(round(i, digits = 3)) THz" for i in multi_freqs]...), legend = :outerright, linewidth = 1.5, size = (650, 500), minorgrid = true, ylabel = "v - w (THz)", xlabel = "T (K)", tickfontsize = 9, legendfontsize = 9)

# ╔═╡ 8f23b6a9-9cb8-47b5-bec1-b7f4f2fe050b
begin
	savefig(multi_v_temp, "multi_v_temp.pdf")
	savefig(multi_w_temp, "multi_w_temp.pdf")
end

# ╔═╡ e7631e4b-dc2c-4b0c-b7a9-7a92dbe0fd1e
begin
	vw_temp = plot(T, v_A, label = "H-A v", legend = :bottomright, linewidth = 2, tickfontsize = 12, labelfontsize = 12, legendfontsize = 12, xlabel = "T(K)", ylabel = "v & w (THz)", size = (550, 550), minorgrid = true, linestyle = :dashdot, color = theme_palette(:default)[1])
	plot!(T, w_A, label = "H-A w", linewidth = 2, linestyle = :dashdot, color = theme_palette(:default)[2])
	plot!(T, v_B, label = "H-B v", linewidth = 2, linestyle = :dash, color = theme_palette(:default)[3])
	plot!(T, w_B, label = "H-B w", linewidth = 2, linestyle = :dash, color = theme_palette(:default)[4])
	plot!(T, v_multi, label = "multi v", linewidth = 2, linestyle = :solid, color = theme_palette(:default)[1])
	plot!(T, w_multi, label = "multi w", linewidth = 2, linestyle = :solid, color = theme_palette(:default)[2])
end

# ╔═╡ 32c47e89-a988-4586-840b-0ac952e1d916
savefig(vw_temp, "vw_temp.pdf")

# ╔═╡ dd0f2762-4cbb-4fc4-b79e-a70b1006905e
begin
	plot(T, v_A, label = "H-A v", legend = :bottomright, linewidth = 2, tickfontsize = 12, labelfontsize = 12, legendfontsize = 12, xlabel = "T(K)", ylabel = "v & w (THz)", size = (550, 550), minorgrid = true, linestyle = :dashdot, color = theme_palette(:default)[1])
	plot!(T, w_A, label = "H-A w", linewidth = 2, linestyle = :dashdot, color = theme_palette(:default)[2])
	plot!(T, v_B, label = "H-B v", linewidth = 2, linestyle = :dash, color = theme_palette(:default)[3])
	plot!(T, w_B, label = "H-B w", linewidth = 2, linestyle = :dash, color = theme_palette(:default)[4])
end

# ╔═╡ 9d039885-65fa-4f43-ac15-ba927a68be60
multi_conductivity_data = parse.(ComplexF64, (CSV.File("multi_conductivity_data.csv") |> Tables.matrix)[2:end, 2:end])

# ╔═╡ 88c443f2-f6ea-464e-9190-aa158a0f7f50
A_conductivity_data = parse.(ComplexF64, (CSV.File("A_conductivity_data.csv") |> Tables.matrix)[2:end, 2:end])

# ╔═╡ db784d36-b859-4022-b585-0730bf081a65
B_conductivity_data = parse.(ComplexF64, (CSV.File("B_conductivity_data.csv") |> Tables.matrix)[2:end, 2:end])

# ╔═╡ 656db2a0-91d6-447b-95a8-2a2981979645
A_contour_real = Plots.contourf(Ω, T, real.(A_conductivity_data)', xlabel = "Ω (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (0, 500))

# ╔═╡ a586ecee-fd54-44e7-8dea-fa1f2a6829df
A_contour_imag = Plots.contourf(Ω, T, imag.(A_conductivity_data)', xlabel = "Ω (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), clims = (-6, 500), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm)

# ╔═╡ 2aeb1ede-b786-494f-9f15-0a8382943c29
A_contour_abs = Plots.contourf(Ω, T, abs.(A_conductivity_data)', xlabel = "Ω (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true),  clims = (0, 500), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm)

# ╔═╡ 88c2fd6b-a8cc-436a-b9e9-5412c19b42bd
begin
	savefig(A_contour_real, "A_contour_real.pdf")
	savefig(A_contour_imag, "A_contour_imag.pdf")
	savefig(A_contour_abs, "A_contour_abs.pdf")
end

# ╔═╡ 329f1be8-967d-4c5b-8536-f6f1008f2a9c
B_contour_real = Plots.contourf(Ω, T, real.(B_conductivity_data)', xlabel = "Ω (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (0, 500))

# ╔═╡ ba30cda1-9a2e-4427-9cde-3a5e3e7e4dcf
B_contour_imag = Plots.contourf(Ω, T, imag.(B_conductivity_data)', xlabel = "Ω (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (-6, 500))

# ╔═╡ e30b9ee1-0eac-44ad-bb6e-4c243672468e
B_contour_abs = Plots.contourf(Ω, T, abs.(B_conductivity_data)', xlabel = "Ω (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (0, 500))

# ╔═╡ 26ec12fd-5c2d-458c-8282-b411ef5184a8
begin
	savefig(B_contour_real, "B_contour_real.pdf")
	savefig(B_contour_imag, "B_contour_imag.pdf")
	savefig(B_contour_abs, "B_contour_abs.pdf")
end

# ╔═╡ 1e89541e-0635-47d9-96f6-99dceef27d00
multi_contour_real = Plots.contourf(Ω, T, real.(multi_conductivity_data)', xlabel = "Ω (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (0, 500))

# ╔═╡ 6b617a87-f84e-4818-be31-5806d2e8894d
multi_contour_imag = Plots.contourf(Ω, T, imag.(multi_conductivity_data)', xlabel = "Ω (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (-6, 500))

# ╔═╡ bc190237-ddff-4abc-bf29-0e23d10fad0d
multi_contour_abs = Plots.contourf(Ω, T, abs.(multi_conductivity_data)', xlabel = "Ω (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true),  linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (0, 500))

# ╔═╡ 512e223b-249d-48ff-9eec-4282f4658c76
begin
	savefig(multi_contour_real, "multi_contour_real.pdf")
	savefig(multi_contour_imag, "multi_contour_imag.pdf")
	savefig(multi_contour_abs, "multi_contour_abs.pdf")
end

# ╔═╡ 0ec007f8-9b46-4d0c-a8d2-c12ff85b0f5a
AB_contour_real = Plots.contourf(Ω, T, log.(abs.(real.(A_conductivity_data .- B_conductivity_data)))', xlabel = "Ω (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true, scale = :log), clims = (-16, 5), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 10, right_margin = 7Plots.mm)

# ╔═╡ 7a68a5ca-2c8f-40bc-ac7b-8a996f48b41b
AB_contour_imag = Plots.contourf(Ω, T, log.(abs.(imag.(A_conductivity_data .- B_conductivity_data)))', xlabel = "Ω (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true, scale = :log), clims = (-16, 5), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 10, right_margin = 7Plots.mm)

# ╔═╡ 3a14ec9b-c726-40da-a4ed-ad22ee89e8d3
AB_contour_abs = Plots.contourf(Ω, T, log.(abs.(A_conductivity_data .- B_conductivity_data))', xlabel = "Ω (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true, scale = :log), clims = (-16, 5), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 10, right_margin = 7Plots.mm)

# ╔═╡ 6071f2a4-ec70-4a12-b40a-84062ae2bb6a
begin
	savefig(AB_contour_real, "AB_contour_real.pdf")
	savefig(AB_contour_imag, "AB_contour_imag.pdf")
	savefig(AB_contour_abs, "AB_contour_abs.pdf")
end

# ╔═╡ b80aedaf-b264-47b3-a679-1c73ce1e8a0f
real.(multi_conductivity_data)[:, 291]

# ╔═╡ c30386e2-e275-4e8c-bc46-e043ce3484b2
real.(A_conductivity_data)[:, 291]

# ╔═╡ 585d49c7-0b54-4884-922c-84e9ae387837
real.(B_conductivity_data)[:, 291]

# ╔═╡ c3281d7d-3d8e-4283-ad14-0fbe30759990
begin
	multi_plot = Plots.plot(Ω[1:250], real.(multi_conductivity_data[1:250, 1:1:30]), legend = false, label = hcat(["T = $i K" for i in T[1:1:30]]...), minorgrid = true, linewidth = 1.5, size = (550, 500), ylabel = L"\mathrm{Multi\ Mobility\ } (cm^2/Vs)", xlabel = L"\nu\ (THz)", tickfontsize = 9, legendfontsize = 9, ylims = (0, 1500), palette = palette(:thermal, 30))
end

# ╔═╡ aec2c420-5d16-4130-914b-0ba94f4cd9ca
red = hcat([(real.(multi_conductivity_data[i, 1]) .- real.(multi_conductivity_data[i-1, 10])) for i in 2:350]...)

# ╔═╡ 32773157-09b8-4feb-ae44-7f41e244b101
imd = hcat([(imag.(multi_conductivity_data[i, 10]) .- imag.(multi_conductivity_data[i-1, 10])) for i in 2:350]...)

# ╔═╡ 41b5f3a6-1fb7-4cab-b712-57816e8209c5
begin
	plot(Ω[1:250], red[1:250], size = (550, 500), minorgrid = true)
	plot!(Ω[1:250], -imd[1:250])
end

# ╔═╡ a64cfd0a-67f6-4941-8ba4-2e26b3361040
begin
	multi_plot2 = Plots.plot(Ω[1:250], imag.(multi_conductivity_data[1:250, 1:1:30]), legend = false, label = hcat(["T = $i K" for i in T[1:1:30]]...), minorgrid = true, linewidth = 1.5, size = (550, 500), ylabel = L"\mathrm{Multi\ Mobility\ } (cm^2/Vs)", xlabel = L"\nu\ (THz)", tickfontsize = 9, legendfontsize = 9, ylims = (0, 1500), palette = palette(:thermal, 30))
end

# ╔═╡ 835feade-6750-48ef-baab-be9dce52500e
Plots.plot(Ω ./ 2π, abs.(real.(A_conductivity_data))[:, 400], legend = true, label = hcat(["T = $i K" for i in T[1]]...), yaxis = :log, xlabel = L"\nu\ (THz)", ylabel = L"\mathrm{Hellwarth\ A\ Mobility\ } (cm^2/Vs)", linewidth = 1.5, minorgrid = true, ylims = (5, 10^2.8), size = (500, 500), tickfontsize = 9, legendfontsize = 9)

# ╔═╡ 4bab65b0-a561-4355-8075-8cd3d8ba09f1
Plots.plot(Ω, abs.(real.(B_conductivity_data))[:, 1:60:end], legend = true, label = hcat(["T = $i K" for i in T[1:60:end]]...), yaxis = :log, xlabel = L"\Omega\ (THz)", ylabel = L"\mathrm{Hellwarth\ B\ Mobility\ } (cm^2/Vs)", linewidth = 1.5, minorgrid = true, ylims = (5, 10^2.8), size = (500, 500), tickfontsize = 9, legendfontsize = 9)

# ╔═╡ f4f970e4-6113-4782-b6dc-9072254c5b4e
Plots.plot(T, real.(multi_conductivity_data)[1:25:end, :]', legend = :topright, label = hcat(["Ω = $i THz" for i in Ω[1:25:end]]...), xlabel = L"\textrm{T\ (K)}", ylabel = L"\textbf{Multi\ Mobility\ (cm^2/Vs)}", yaxis = :log, minorgrid = true, linewidth = 1.5, ylims = (5, 10^4.5), size = (500, 500), tickfontsize = 9, legendfontsize = 9, xaxis = :log, xlims = (10, 400))

# ╔═╡ 86ff3a78-852f-424a-af6d-32e886726d0a
plot!(T, T.^(-0.2) * 18^2, linestyle = :dash, linewidth = 3, label  = "T^(-0.2)")

# ╔═╡ efca960d-a7d6-42df-8196-7e00ae0ac024
Plots.plot(T, real.(A_conductivity_data)[1:25:end, :]', legend = :topright, label = hcat(["Ω = $i THz" for i in Ω[1:25:end]]...), xlabel = L"\textrm{T\ (K)}", ylabel = L"\textrm{Hellwarth\ A\ Mobility\ (cm^2/Vs)}", yaxis = :log, minorgrid = true, linewidth = 1.5, ylims = (5, 10^5), size = (500, 500), tickfontsize = 9, legendfontsize = 9)

# ╔═╡ da3417bb-2570-41df-ac54-bdf1a46299d9
Plots.plot(T, real.(B_conductivity_data)[1:25:end, :]', legend = true, label = hcat(["Ω = $i THz" for i in Ω[1:25:end]]...), xlabel = L"\textrm{T\ (K)}", ylabel = L"\textrm{Hellwarth\ B\ Mobility\ (cm^2/Vs)}", yaxis = :log, minorgrid = true, linewidth = 1.5, ylims = (5, 10^5), size = (500, 500), tickfontsize = 9, legendfontsize = 9, xaxis = :log, xlims = (10, 400))

# ╔═╡ 5c685b4e-8702-4b74-90e8-25ffd53e5582
plot!(T, T.^(-0.5) * 50^2, linestyle = :dash, linewidth = 3, label  = "T^(-0.2)", legend = :outerright, size = (700, 500))

# ╔═╡ 9bd4dd63-c784-4a1e-86fa-935bd8c25ad2
function ridgeline(x, y, z; shift = 1, jump = 10, ylabel = :none, xlabel = :none, palette = :twilight, linewidth = 1.5, color = :black, size = (595, 842), ymirror = false, fillalpha = 1, step = 1)
	if !ymirror
		p = plot(x, z[:, 1], fillrange = z[:, jump] .- shift, fillalpha = fillalpha, ylabel = ylabel, xlabel = xlabel, xticks = (0:maximum(x)/10:maximum(x), string.(x[1:Int((length(x)-1)/10):end])), yticks = (z[1, 1:jump * step:end] .- shift .* range(0, stop = length(y[1:jump:end]) - 1, step = step), string.(y[1:jump * step:end])), legend = false, size = size, grid = false, ymirror = ymirror, palette = palette, tickfontsize = 12, labelfontsize = 12)
	else
		p = plot(x, z[:, 1], fillrange = z[:, jump] .- shift, fillalpha = fillalpha, ylabel = ylabel, xlabel = xlabel, xticks = (0:maximum(x)/10:maximum(x), string.(x[1:Int((length(x))/10):end])), yticks = (z[end, 1:jump * step:end] .- shift .* range(0, stop = length(y[1:jump:end]) - 1, step = step), string.(y[1:jump * step:end])), legend = false, size = size, grid = false, ymirror = ymirror, palette = palette, tickfontsize = 12, labelfontsize = 12)
	end

	plot!(x, z[:, 1], color = color, linewidth = linewidth)
	plot!(x, z[:, jump] .- shift, color = color, linewidth = linewidth)

	for i in 1:Int(floor(length(y)/jump) - 1)
	plot!(x, z[:, i * jump] .- shift * i, fillrange = z[:, (i + 1) * jump] .- shift * (i + 1), fillalpha = fillalpha, palette = palette)
	plot!(x, z[:, i * jump] .- shift * i, color = color, linewidth = linewidth)
	plot!(x, z[:, (i + 1) * jump] .- shift * (i + 1), color = color, linewidth = linewidth)
	end
	return p
end

# ╔═╡ 7ae0d8a1-d2ad-4599-bce1-f71865fa7497
multi_plot_temp_real = ridgeline(round.(Ω, digits = 1), Int.(T), log.(real.(multi_conductivity_data)), jump = 10, shift = 0.6, xlabel = "Ω (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ 23c44da4-8468-4987-8d5e-1069ffb41d13
multi_plot_temp_imag = ridgeline(round.(Ω, digits = 1), Int.(T), log.(abs.(imag.(multi_conductivity_data))), jump = 10, shift = 0.6, xlabel = "Ω (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 7)

# ╔═╡ 077a2c40-9aae-4045-89da-99f1db15fae7
multi_plot_temp_abs = ridgeline(round.(Ω, digits = 1), Int.(T), log.(abs.(multi_conductivity_data)), jump = 10, shift = 0.6, xlabel = "Ω (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ fbf4e5d9-c20b-4953-8477-c0144e34e9e9
begin
	savefig(multi_plot_temp_real, "multi_plot_temp_real.pdf")
	savefig(multi_plot_temp_imag, "multi_plot_temp_imag.pdf")
	savefig(multi_plot_temp_abs, "multi_plot_temp_abs.pdf")
end

# ╔═╡ bfd794a2-abf9-499e-92fe-aa02f262171d
A_plot_temp_real = ridgeline(round.(Ω, digits = 1), Int.(T), log.(abs.(real.(A_conductivity_data))), jump = 10, shift = 0.6, xlabel = "Ω (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ cbdee07e-77e5-4113-b099-4feb30d5d74a
A_plot_temp_imag = ridgeline(round.(Ω, digits = 1), Int.(T), log.(abs.(imag.(A_conductivity_data))), jump = 10, shift = 0.6, xlabel = "Ω (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 7)

# ╔═╡ 2b4876d2-6d1e-4774-89c1-3ff8be380e11
A_plot_temp_abs = ridgeline(round.(Ω, digits = 1), Int.(T), log.(abs.(A_conductivity_data)), jump = 10, shift = 0.6, xlabel = "Ω (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ 79ddc218-4710-466a-a5ba-976431cd7165
begin
	savefig(A_plot_temp_real, "A_plot_temp_real.pdf")
	savefig(A_plot_temp_imag, "A_plot_temp_imag.pdf")
	savefig(A_plot_temp_abs, "A_plot_temp_abs.pdf")
end

# ╔═╡ 9b8183fb-3a65-463b-bae5-62a8aa9f06ac
B_plot_temp_real = ridgeline(round.(Ω, digits = 1), Int.(T), log.(abs.(real.(B_conductivity_data))), jump = 10, shift = 0.6, xlabel = "Ω (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ ed8f7ce5-c4eb-44da-a747-56ac7c152760
B_plot_temp_imag = ridgeline(round.(Ω, digits = 1), Int.(T), log.(abs.(imag.(B_conductivity_data))), jump = 10, shift = 0.6, xlabel = "Ω (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 7)

# ╔═╡ 80cde1da-7d6a-4fea-a2c2-89b2a0a6cd0d
B_plot_temp_abs = ridgeline(round.(Ω, digits = 1), Int.(T), log.(abs.(B_conductivity_data)), jump = 10, shift = 0.6, xlabel = "Ω (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ 4cfb19f7-6852-4224-9b09-54a844624cde
begin
	savefig(B_plot_temp_real, "B_plot_temp_real.pdf")
	savefig(B_plot_temp_imag, "B_plot_temp_imag.pdf")
	savefig(B_plot_temp_abs, "B_plot_temp_abs.pdf")
end

# ╔═╡ ccfa612b-9cf0-4e4c-a69c-e3b002e5fea5
AB_plot_temp_real = ridgeline(round.(Ω, digits = 1), Int.(T), log.(abs.(real.(A_conductivity_data .- B_conductivity_data))), jump = 5, shift = 0.3, xlabel = "Ω (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ 39eeac1c-5db8-4d6c-8552-d308b5199d5e
AB_plot_temp_imag = ridgeline(round.(Ω, digits = 1), Int.(T), log.(abs.(imag.(A_conductivity_data .- B_conductivity_data))), jump = 5, shift = 0.3, xlabel = "Ω (THz)", ylabel = "T (K)", size = (500, 500), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 4)

# ╔═╡ a8ab6e8a-2d89-4383-9ab6-e49990712489
AB_plot_temp_abs = ridgeline(round.(Ω, digits = 1), Int.(T), log.(abs.(A_conductivity_data .- B_conductivity_data)), jump = 5, shift = 0.3, xlabel = "Ω (THz)", ylabel = "T (K)", size = (500, 500), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 4)

# ╔═╡ ba68e7aa-963f-4cdb-a141-b4ee9ff697db
begin
	savefig(AB_plot_temp_real, "AB_plot_temp_real.pdf")
	savefig(AB_plot_temp_imag, "AB_plot_temp_imag.pdf")
	savefig(AB_plot_temp_abs, "AB_plot_temp_abs.pdf")
end

# ╔═╡ 9e0e8eea-f5e9-4d91-80cd-a1e7847933e6
multi_plot_freq_real = ridgeline(Int.(T), round.(Ω, digits = 1), log.(abs.(real.(multi_conductivity_data)))', jump = 6, shift = 0.4, ylabel = "Ω (THz)", xlabel = "T (K)", ymirror = true, size = (550, 550), linewidth = 0.75, palette = :twilight, fillalpha = 0.75, step = 7)

# ╔═╡ 3e185071-6fda-4dc6-8e84-4900a106b798
multi_plot_freq_imag = ridgeline(Int.(T), round.(Ω, digits = 1), log.(abs.(imag.(multi_conductivity_data)))', jump = 6, shift = 0.4, ylabel = "Ω (THz)", xlabel = "T (K)", ymirror = true, size = (550, 550), linewidth = 0.75, palette = :twilight, fillalpha = 0.75, step = 10)

# ╔═╡ 05516d71-0f45-4e97-97b0-8a1819aed270
multi_plot_freq_abs = ridgeline(Int.(T), round.(Ω, digits = 1), log.(abs.(multi_conductivity_data))', jump = 6, shift = 0.4, ylabel = "Ω (THz)", xlabel = "T (K)", ymirror = true, size = (500, 550), linewidth = 0.75, palette = :twilight, fillalpha = 0.75, step = 7)

# ╔═╡ c3f82311-03fa-4cb0-9dca-4cf495b5e8ca
begin
	savefig(multi_plot_freq_real, "multi_plot_freq_real.pdf")
	savefig(multi_plot_freq_imag, "multi_plot_freq_imag.pdf")
	savefig(multi_plot_freq_abs, "multi_plot_freq_abs.pdf")
end

# ╔═╡ c3bdb1f6-889b-4c84-9857-88f619d6e617
A_plot_freq_real = ridgeline(Int.(T), round.(Ω, digits = 1), log.(abs.(real.(A_conductivity_data)))', jump = 6, shift = 0.6, ylabel = "Ω (THz)", xlabel = "T (K)", ymirror = true, size = (550, 550), linewidth = 0.75, palette = :twilight, fillalpha = 0.75, step = 7)

# ╔═╡ 428b7174-9d3f-4762-8941-6feabcfbb4e0
A_plot_freq_imag = ridgeline(Int.(T), round.(Ω, digits = 1), log.(abs.(imag.(A_conductivity_data)))', jump = 6, shift = 0.6, ylabel = "Ω (THz)", xlabel = "T (K)", ymirror = true, size = (550, 550), linewidth = 0.75, palette = :twilight, fillalpha = 0.75, step = 10)

# ╔═╡ 5900719d-e804-4a7d-a142-38fddf56a085
A_plot_freq_abs = ridgeline(Int.(T), round.(Ω, digits = 1), log.(abs.(A_conductivity_data))', jump = 6, shift = 0.6, ylabel = "Ω (THz)", xlabel = "T (K)", ymirror = true, size = (550, 550), linewidth = 0.75, palette = :twilight, fillalpha = 0.75, step = 7)

# ╔═╡ ff9a0b04-b007-468e-acba-677a7f25a340
begin
	savefig(A_plot_freq_real, "A_plot_freq_real.pdf")
	savefig(A_plot_freq_imag, "A_plot_freq_imag.pdf")
	savefig(A_plot_freq_abs, "A_plot_freq_abs.pdf")
end

# ╔═╡ bab7f306-94d2-4360-8c0d-08b8a619046d
B_plot_freq_real = ridgeline(Int.(T), round.(Ω, digits = 1), log.(abs.(real.(B_conductivity_data)))', jump = 6, shift = 0.6, ylabel = "Ω (THz)", xlabel = "T (K)", ymirror = true, size = (550, 550), linewidth = 0.75, palette = :twilight, fillalpha = 0.75, step = 7)

# ╔═╡ 0879342a-c552-4bc6-8dce-21f020abeb25
B_plot_freq_imag = ridgeline(Int.(T), round.(Ω, digits = 1), log.(abs.(imag.(B_conductivity_data)))', jump = 6, shift = 0.6, ylabel = "Ω (THz)", xlabel = "T (K)", ymirror = true, size = (550, 550), linewidth = 0.75, palette = :twilight, fillalpha = 0.75, step = 10)

# ╔═╡ 8acd98e2-23ef-478b-97c2-533e0e2bef98
B_plot_freq_abs = ridgeline(Int.(T), round.(Ω, digits = 1), log.(abs.(B_conductivity_data))', jump = 6, shift = 0.6, ylabel = "Ω (THz)", xlabel = "T (K)", ymirror = true, size = (550, 550), linewidth = 0.75, palette = :twilight, fillalpha = 0.75, step = 7)

# ╔═╡ 34411f71-63bc-4c91-9cd4-486e69f5876a
begin
	savefig(B_plot_freq_real, "B_plot_freq_real.pdf")
	savefig(B_plot_freq_imag, "B_plot_freq_imag.pdf")
	savefig(B_plot_freq_abs, "B_plot_freq_abs.pdf")
end

# ╔═╡ aff41df3-144f-452e-b059-c306000ad10d
AB_plot_freq_real = ridgeline(Int.(T), round.(Ω, digits = 1), log.(abs.(real.(A_conductivity_data .- B_conductivity_data)))', jump = 3, shift = 1, ylabel = "Ω (THz)", xlabel = "T (K)", ymirror = true, size = (500, 500), linewidth = 0.75, palette = :twilight, fillalpha = 0.75, step = 5)

# ╔═╡ 9288099f-4e2a-43d6-ba2c-1debb17c4901
AB_plot_freq_imag = ridgeline(Int.(T), round.(Ω, digits = 1), log.(abs.(imag.(A_conductivity_data .- B_conductivity_data)))', jump = 3, shift = 1, ylabel = "Ω (THz)", xlabel = "T (K)", ymirror = true, size = (500, 500), linewidth = 0.75, palette = :twilight, fillalpha = 0.75, step = 5)

# ╔═╡ a844da51-68ad-442f-9517-bbcde197c38e
AB_plot_freq_abs = ridgeline(Int.(T), round.(Ω, digits = 1), log.(abs.(A_conductivity_data .- B_conductivity_data))', jump = 3, shift = 1, ylabel = "Ω (THz)", xlabel = "T (K)", ymirror = true, size = (500, 500), linewidth = 0.75, palette = :twilight, fillalpha = 0.75, step = 5)

# ╔═╡ f47f71d3-870a-47b8-9ab2-426b8ceacf72
begin
	savefig(AB_plot_freq_real, "AB_plot_freq_real.pdf")
	savefig(AB_plot_freq_imag, "AB_plot_freq_imag.pdf")
	savefig(AB_plot_freq_abs, "AB_plot_freq_abs.pdf")
end

# ╔═╡ 8b7c9fd6-2289-47c6-ad73-2e60e867979e
plot(multi_freqs, F_multi[:, 1:1:20], markershape = :circle, legend = :bottomright, label = hcat(["$(T[i]) K" for i in 1:1:20]...), xticks = (multi_freqs, hcat(["$(round(i, digits = 2))" for i in multi_freqs]...)), xrotation = 90, tickfontsize = 9, ylabel = "Free energy (meV)", xlabel = "Phonon frequencies (THz)", linestyle = :dot)

# ╔═╡ 509b518e-07d9-4b5d-8938-b4a9ec188d8f
ridgeline(round.(multi_freqs, digits = 1), T[4:400], F_multi[:, 4:400], ymirror = true)

# ╔═╡ 18edac61-b8e1-4119-b49d-56d2b9ba2251
begin
	plot(multi_freqs, v_multi[:, 1:50:400] .* multi_freqs, markershape = :circle, legend = :outerright, label = hcat(["$(T[i]) K" for i in 1:50:400]...), xticks = (multi_freqs, hcat(["$(round(i, digits = 2))" for i in multi_freqs]...)), ylims = (0, 65), xrotation = 90, tickfontsize = 9, ylabel = "v & w (THz)", xlabel = "Phonon frequencies (THz)", linestyle = :dot)
	plot!(multi_freqs, w_multi[:, 1:50:400] .* multi_freqs, markershape = :diamond, legend = :outerright, label = :none, xticks = (multi_freqs, hcat(["$(round(i, digits = 2))" for i in multi_freqs]...)), ylims = (0, 65), xrotation = 90, tickfontsize = 9, linestyle = :dash, color = hcat([theme_palette(:default)[i] for i in 1:8]...))
end

# ╔═╡ 6dae57f3-cfd4-4598-bc1b-ea25f8d39565
plot(multi_freqs, (v_multi[:, 1:50:400] .- w_multi[:, 1:50:400]) .* multi_freqs ./ ir_activity, markershape = :circle, legend = :outerright, label = hcat(["$(T[i]) K" for i in 1:50:400]...), xticks = (multi_freqs, hcat(["$(round(i, digits = 2))" for i in multi_freqs]...)), xrotation = 90, tickfontsize = 9, ylabel = "(v - w) / IR", xlabel = "Phonon frequencies (THz)", linestyle = :dot)

# ╔═╡ Cell order:
# ╠═ed9ab030-e61d-11eb-32f3-5be220868ae7
# ╠═d51f978f-a2c8-410a-8836-f57fc362eed1
# ╠═ffa81468-a1cd-44f3-aa5a-c970e56b7f27
# ╠═4c98c769-0aed-4b8a-bf92-fb06268100a9
# ╠═c819bc8f-cc85-4edb-8289-5884f7fbb8f7
# ╠═88fe20b4-303a-494e-95cc-e8183b2f474e
# ╠═a76c9858-0fe7-45c2-9b5c-a8cf7f800c08
# ╠═6d08ec24-51d5-4105-bddf-443bc8496456
# ╠═5fa3101f-d24b-41f1-9a77-f5b7782ec4a3
# ╠═78156817-9c43-44c0-90e9-55e1100538a7
# ╠═4bf7a882-1306-415a-8760-6191e5284306
# ╠═8b91053e-28cc-47ce-8551-72e3e3c8ee93
# ╠═fd24b1c2-a6c2-4328-854b-ada51d40c56a
# ╠═db273a6b-bd14-4935-92dd-1271def823ba
# ╠═830ce1e4-80f4-4dfa-a214-d7d0c84cfbbb
# ╠═4b5dc310-b4f6-4bd3-88e5-f082ea8ee47e
# ╠═f9e47e6f-fdc5-42b6-9ea2-46c801f8cf6b
# ╠═e2a502b1-cca9-40c2-a107-2b577f0662dc
# ╠═92f361d6-2052-45bd-a2a1-0d19d04eec1b
# ╠═b1473f00-5aa9-4ace-a0f9-9e5cb9747af5
# ╠═4c4d1b61-6710-478c-a8f5-067c219be342
# ╠═9493c1f8-d2b4-42b3-9af5-2d69ee7499bf
# ╠═d4a13a12-e1ca-4ded-9c49-3c33b1cb69a0
# ╠═40fc4b41-c835-41a3-b468-765d4bbd0042
# ╠═aa1b7720-91ee-48e5-a97e-f8e3a499ad2f
# ╠═10b31119-9d2f-445d-9353-9e6788a464d3
# ╠═eef00978-0153-4031-b748-d9e1e5e2f5c5
# ╠═45e2101a-010d-404d-aaf6-c68d28536677
# ╠═4da848c3-23dd-4640-bab3-a8b4b4784639
# ╠═8f23b6a9-9cb8-47b5-bec1-b7f4f2fe050b
# ╠═e7631e4b-dc2c-4b0c-b7a9-7a92dbe0fd1e
# ╠═32c47e89-a988-4586-840b-0ac952e1d916
# ╠═dd0f2762-4cbb-4fc4-b79e-a70b1006905e
# ╠═9d039885-65fa-4f43-ac15-ba927a68be60
# ╠═88c443f2-f6ea-464e-9190-aa158a0f7f50
# ╠═db784d36-b859-4022-b585-0730bf081a65
# ╠═656db2a0-91d6-447b-95a8-2a2981979645
# ╠═a586ecee-fd54-44e7-8dea-fa1f2a6829df
# ╠═2aeb1ede-b786-494f-9f15-0a8382943c29
# ╠═88c2fd6b-a8cc-436a-b9e9-5412c19b42bd
# ╠═329f1be8-967d-4c5b-8536-f6f1008f2a9c
# ╠═ba30cda1-9a2e-4427-9cde-3a5e3e7e4dcf
# ╠═e30b9ee1-0eac-44ad-bb6e-4c243672468e
# ╠═26ec12fd-5c2d-458c-8282-b411ef5184a8
# ╠═1e89541e-0635-47d9-96f6-99dceef27d00
# ╠═6b617a87-f84e-4818-be31-5806d2e8894d
# ╠═bc190237-ddff-4abc-bf29-0e23d10fad0d
# ╠═512e223b-249d-48ff-9eec-4282f4658c76
# ╠═0ec007f8-9b46-4d0c-a8d2-c12ff85b0f5a
# ╠═7a68a5ca-2c8f-40bc-ac7b-8a996f48b41b
# ╠═3a14ec9b-c726-40da-a4ed-ad22ee89e8d3
# ╠═6071f2a4-ec70-4a12-b40a-84062ae2bb6a
# ╠═b80aedaf-b264-47b3-a679-1c73ce1e8a0f
# ╠═c30386e2-e275-4e8c-bc46-e043ce3484b2
# ╠═585d49c7-0b54-4884-922c-84e9ae387837
# ╠═c3281d7d-3d8e-4283-ad14-0fbe30759990
# ╠═aec2c420-5d16-4130-914b-0ba94f4cd9ca
# ╠═32773157-09b8-4feb-ae44-7f41e244b101
# ╠═41b5f3a6-1fb7-4cab-b712-57816e8209c5
# ╠═a64cfd0a-67f6-4941-8ba4-2e26b3361040
# ╠═835feade-6750-48ef-baab-be9dce52500e
# ╠═4bab65b0-a561-4355-8075-8cd3d8ba09f1
# ╠═f4f970e4-6113-4782-b6dc-9072254c5b4e
# ╠═86ff3a78-852f-424a-af6d-32e886726d0a
# ╠═efca960d-a7d6-42df-8196-7e00ae0ac024
# ╠═da3417bb-2570-41df-ac54-bdf1a46299d9
# ╠═5c685b4e-8702-4b74-90e8-25ffd53e5582
# ╠═e2d229a0-eaae-481c-9ee1-94d641fc3014
# ╠═9bd4dd63-c784-4a1e-86fa-935bd8c25ad2
# ╠═7ae0d8a1-d2ad-4599-bce1-f71865fa7497
# ╠═23c44da4-8468-4987-8d5e-1069ffb41d13
# ╠═077a2c40-9aae-4045-89da-99f1db15fae7
# ╠═fbf4e5d9-c20b-4953-8477-c0144e34e9e9
# ╠═bfd794a2-abf9-499e-92fe-aa02f262171d
# ╠═cbdee07e-77e5-4113-b099-4feb30d5d74a
# ╠═2b4876d2-6d1e-4774-89c1-3ff8be380e11
# ╠═79ddc218-4710-466a-a5ba-976431cd7165
# ╠═9b8183fb-3a65-463b-bae5-62a8aa9f06ac
# ╠═ed8f7ce5-c4eb-44da-a747-56ac7c152760
# ╠═80cde1da-7d6a-4fea-a2c2-89b2a0a6cd0d
# ╠═4cfb19f7-6852-4224-9b09-54a844624cde
# ╠═ccfa612b-9cf0-4e4c-a69c-e3b002e5fea5
# ╠═39eeac1c-5db8-4d6c-8552-d308b5199d5e
# ╠═a8ab6e8a-2d89-4383-9ab6-e49990712489
# ╠═ba68e7aa-963f-4cdb-a141-b4ee9ff697db
# ╠═9e0e8eea-f5e9-4d91-80cd-a1e7847933e6
# ╠═3e185071-6fda-4dc6-8e84-4900a106b798
# ╠═05516d71-0f45-4e97-97b0-8a1819aed270
# ╠═c3f82311-03fa-4cb0-9dca-4cf495b5e8ca
# ╠═c3bdb1f6-889b-4c84-9857-88f619d6e617
# ╠═428b7174-9d3f-4762-8941-6feabcfbb4e0
# ╠═5900719d-e804-4a7d-a142-38fddf56a085
# ╠═ff9a0b04-b007-468e-acba-677a7f25a340
# ╠═bab7f306-94d2-4360-8c0d-08b8a619046d
# ╠═0879342a-c552-4bc6-8dce-21f020abeb25
# ╠═8acd98e2-23ef-478b-97c2-533e0e2bef98
# ╠═34411f71-63bc-4c91-9cd4-486e69f5876a
# ╠═aff41df3-144f-452e-b059-c306000ad10d
# ╠═9288099f-4e2a-43d6-ba2c-1debb17c4901
# ╠═a844da51-68ad-442f-9517-bbcde197c38e
# ╠═f47f71d3-870a-47b8-9ab2-426b8ceacf72
# ╠═8b7c9fd6-2289-47c6-ad73-2e60e867979e
# ╠═509b518e-07d9-4b5d-8938-b4a9ec188d8f
# ╠═18edac61-b8e1-4119-b49d-56d2b9ba2251
# ╠═6dae57f3-cfd4-4598-bc1b-ea25f8d39565


text:
@everywhere using CSV, GLMNet, Combinatorics, DataFrames, Plots, ArgParse
#@everywhere using ScikitLearn: fit!, predict, @sk_import, fit_transform!
#@everywhere @sk_import ensemble: RandomForestRegressor
include("../forc-helpers/OLSestimator.jl")
include("../forc-helpers/create_matrices.jl")
include("../forc-models/all_models.jl")

function prepare(level, h, J, oos_i)
    ## Choose the head line price to construct the inflation rates
    head_line = true
    J_cv = (J + oos_i + h_lag + h) ÷ 2
    println("h = $h, J = $J, J_cv = $J_cv")
    Y, Y_lag, X_lag, X_lag2 = create_matrices(level, h, h_lag, head_line)
    return head_line, Y, Y_lag, X_lag, X_lag2, J_cv
end

max_stars_repo_name: jungangl/forecast-machine-learning

max_issues_repo_name: jungangl/forecast-machine-learning

max_forks_repo_name: jungangl/forecast-machine-learning


text:
@info "Starting precompile exec."

using Distributions, LinearAlgebra, DensityInterface, InverseFunctions
using BAT, ValueShapes
import Cuba

prior = HierarchicalDistribution(
    NamedTupleDist(a = Normal(), b = Weibull(), c = 5)
) do x
    NamedTupleDist(y = Normal(x.a, x.b), z = MvNormal([1.3 0.5; 0.5 2.2]))
end
likelihood = logfuncdensity((x -> -norm(x)^2) ∘ inverse(varshape(prior)))
density = PosteriorDensity(likelihood, prior)

samples_mh = bat_sample(density, MCMCSampling(mcalg = MetropolisHastings(), nchains = 4, nsteps = 10^4, strict = false)).result

samples_hmc = bat_sample(density, MCMCSampling(mcalg = HamiltonianMC(), nchains = 4, nsteps = 10^4, strict = false)).result
@show sd = SampledDensity(density, samples_hmc)
@show ess = bat_eff_sample_size(samples_hmc)

integral_ahmi = bat_integrate(samples_hmc, AHMIntegration()).result
integral_vegas = bat_integrate(density, VEGASIntegration(strict = false)).result

ENV["GKSwstype"] = "nul"
import Plots

for plt in (
    Plots.histogram(randn(1000)),
    Plots.stephist(randn(1000)),
    Plots.histogram2d(randn(1000), randn(1000)),
    Plots.scatter(randn(10), randn(10), randn(10)),
    Plots.scatter(randn(10), randn(10), randn(10)),
    Plots.plot(samples_mh),
    Plots.plot(samples_hmc),
)
    Plots.prepare_output(plt)
    precompile(display, (Plots.PlotsDisplay, typeof(plt)))

    show(IOBuffer(), MIME"image/png"(), plt)
    show(IOBuffer(), MIME"image/svg+xml"(), plt)
end

@info "Precompile exec finished."

max_stars_repo_name: legend-exp/legend-julia-tutorial

max_issues_repo_name: legend-exp/legend-julia-tutorial

max_forks_repo_name: legend-exp/legend-julia-tutorial


text:
class C {
  public static void main(String[] args) {
      int x = ({H})5;

      condition c;// = false;

      int y = ({L /c H}) 9;

  }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
@testset "object" begin

@testset "Object" begin
    mat = Material("material", ε=rand(), μ=rand())

    @testset "Cuboid" begin
        ba = sort(rand(3,2), dims=2)  # array
        b = (ba[:,1], ba[:,2])

        box_mat = Object(Cuboid(b), mat)  # new convenience constructor

        @test shape(box_mat) == Cuboid(b)
        @test max∆l(box_mat) == fill(Inf,3)
        @test matparam(box_mat,EE)==mat.param[nE] && matparam(box_mat,HH)==mat.param[nH]
    end  # @testset "Cuboid"
end  # @testset "Object"

# @testset "periodize" begin
#     # Square lattice
#     Si = Material("Si", ε = 12)
#     ge = PRIM
#     eSi = EncodedMaterial(ge, Si)
#     c_eSi = Object(Cylinder([0,0,0], 1, 5, [0,0,1]), eSi)
#
#     b = Cuboid([0,0,0], [10,10,5])
#     A = [1 0 0; 0 1 0; 0 0 5]'
#     obj_array = periodize(c_eSi, A, b)
#
#     @test length(obj_array) == 11^2
#     @test all(map(x->(x==max∆l(c_eSi)), max∆l.(obj_array)))
#     @test all(map(x->(x==matparam(c_eSi,PRIM)), matparam.(obj_array,PRIM)))
#     @test all(map(x->(x==matparam(c_eSi,DUAL)), matparam.(obj_array,DUAL)))
#
#     oind2obj = Object{3}[]
#     pind2matprm = (S²ComplexF3[], S²ComplexF3[])
#     add_obj!(oind2obj, pind2matprm, obj_array)
# end  # @testset "periodize"

# @testset "union" begin
#     b1 = Cuboid(((0,1),(0,1),(0,1)))
#     b2 = Cuboid(((-1,0),(-1,0),(-1,0)))
#     s = sphere((0,0,0), 1)
#     lsf_union = union(lsf(b1), lsf(s), lsf(b2))
#     @test all(contains.(lsf_union, [(1,1,1), (-1,-1,-1), (-1/√3,-1/√3,1/√3), (1/√3,1/√3,-1/√3)]))
# end
#
# @testset "intersect" begin
#     b1 = Cuboid(((-1,2),(-1,2),(-1,2)))
#     b2 = Cuboid(((-2,1),(-2,1),(-2,1)))
#     lsf_intersect = intersect(lsf(b1), lsf(b2))
#     @test all(contains.(lsf_intersect, [(x,y,z) for x = -1:1, y = -1:1, z = -1:1]))
#     @test all(!contains.(lsf_intersect, [(x,y,z) for x = (-2,2), y = (-2,2), z = (-2,2)]))
# end
#
# @testset "flip" begin
#     b = Cuboid(((0,1),(0,1),(0,1)))
#     lsf0 = lsf(b)
#     lsf_x = flip(lsf0, X̂, 0)
#     lsf_y = flip(lsf0, Ŷ, 0)
#     lsf_z = flip(lsf0, Ẑ, 0)
#
#     @test all(contains.(lsf0, [(x,y,z) for x = 0:1, y = 0:1, z = 0:1]))
#     @test all(contains.(lsf_x, [(x,y,z) for x = -1:0, y = 0:1, z = 0:1]))
#     @test all(contains.(lsf_y, [(x,y,z) for x = 0:1, y = -1:0, z = 0:1]))
#     @test all(contains.(lsf_z, [(x,y,z) for x = 0:1, y = 0:1, z = -1:0]))
# end
#
# @testset "shift" begin
#     b = Cuboid(((0,1),(0,1),(0,1)))
#     lsf0 = lsf(b)
#     lsf_x = shift(lsf0, X̂, -1)
#     lsf_y = shift(lsf0, Ŷ, -1)
#     lsf_z = shift(lsf0, Ẑ, -1)
#
#     @test all(contains.(lsf0, [(x,y,z) for x = 0:1, y = 0:1, z = 0:1]))
#     @test all(contains.(lsf_x, [(x,y,z) for x = -1:0, y = 0:1, z = 0:1]))
#     @test all(contains.(lsf_y, [(x,y,z) for x = 0:1, y = -1:0, z = 0:1]))
#     @test all(contains.(lsf_z, [(x,y,z) for x = 0:1, y = 0:1, z = -1:0]))
# end

end  # @testset "object"


text:
# function map

# functions with original name

const cuInit = :cuInit
const cuDeviceGetCount = :cuDeviceGetCount
const cuDriverGetVersion = :cuDriverGetVersion
const cuDeviceGet = :cuDeviceGet
const cuDeviceGetAttribute = :cuDeviceGetAttribute
const cuDeviceGetCount = :cuDeviceGetCount
const cuDeviceGetName = :cuDeviceGetName
const cuLaunchKernel = :cuLaunchKernel
const cuModuleGetFunction = :cuModuleGetFunction
const cuModuleLoad = :cuModuleLoad
const cuModuleUnload = :cuModuleUnload
const cuStreamSynchronize = :cuStreamSynchronize

# when api version >= 3020
const cuDeviceTotalMem = :cuDeviceTotalMem_v2
const cuCtxCreate = :cuCtxCreate_v2
const cuModuleGetGlobal = :cuModuleGetGlobal_v2
const cuMemGetInfo = :cuMemGetInfo_v2
const cuMemAlloc = :cuMemAlloc_v2
const cuMemAllocPitch = :cuMemAllocPitch_v2
const cuMemFree = :cuMemFree_v2
const cuMemGetAddressRange = :cuMemGetAddressRange_v2
const cuMemAllocHost = :cuMemAllocHost_v2
const cuMemHostGetDevicePointer = :cuMemHostGetDevicePointer_v2
const cuMemcpyHtoD = :cuMemcpyHtoD_v2
const cuMemcpyDtoH = :cuMemcpyDtoH_v2
const cuMemcpyDtoD = :cuMemcpyDtoD_v2
const cuMemcpyDtoA = :cuMemcpyDtoA_v2
const cuMemcpyAtoD = :cuMemcpyAtoD_v2
const cuMemcpyHtoA = :cuMemcpyHtoA_v2
const cuMemcpyAtoH = :cuMemcpyAtoH_v2
const cuMemcpyAtoA = :cuMemcpyAtoA_v2
const cuMemcpyHtoAAsync = :cuMemcpyHtoAAsync_v2
const cuMemcpyAtoHAsync = :cuMemcpyAtoHAsync_v2
const cuMemcpy2D = :cuMemcpy2D_v2
const cuMemcpy2DUnaligned = :cuMemcpy2DUnaligned_v2
const cuMemcpy3D = :cuMemcpy3D_v2
const cuMemcpyHtoDAsync = :cuMemcpyHtoDAsync_v2
const cuMemcpyDtoHAsync = :cuMemcpyDtoHAsync_v2
const cuMemcpyDtoDAsync = :cuMemcpyDtoDAsync_v2
const cuMemcpy2DAsync = :cuMemcpy2DAsync_v2
const cuMemcpy3DAsync = :cuMemcpy3DAsync_v2
const cuMemsetD8 = :cuMemsetD8_v2
const cuMemsetD16 = :cuMemsetD16_v2
const cuMemsetD32 = :cuMemsetD32_v2
const cuMemsetD2D8 = :cuMemsetD2D8_v2
const cuMemsetD2D16 = :cuMemsetD2D16_v2
const cuMemsetD2D32 = :cuMemsetD2D32_v2
const cuArrayCreate = :cuArrayCreate_v2
const cuArrayGetDescriptor = :cuArrayGetDescriptor_v2
const cuArray3DCreate = :cuArray3DCreate_v2
const cuArray3DGetDescriptor = :cuArray3DGetDescriptor_v2
const cuTexRefSetAddress = :cuTexRefSetAddress_v2
const cuTexRefGetAddress = :cuTexRefGetAddress_v2
const cuGraphicsResourceGetMappedPointer = :cuGraphicsResourceGetMappedPointer_v2

# when api version >= 4000
const cuCtxDestroy = :cuCtxDestroy_v2
const cuCtxPopCurrent = :cuCtxPopCurrent_v2
const cuCtxPushCurrent = :cuCtxPushCurrent_v2
const cuStreamDestroy = :cuStreamDestroy_v2
const cuEventDestroy = :cuEventDestroy_v2


text:
function peerconnection()
  """
      window.moz = !!navigator.mozGetUserMedia;
    var chromeVersion = !!navigator.mozGetUserMedia ? 0 : parseInt(navigator.userAgent.match( /Chrom(e|ium)\/([0-9]+)\./ )[2]);

    function RTCPeerConnection(options) {
        var w = window,
            PeerConnection = w.mozRTCPeerConnection || w.webkitRTCPeerConnection,
            SessionDescription = w.mozRTCSessionDescription || w.RTCSessionDescription,
            IceCandidate = w.mozRTCIceCandidate || w.RTCIceCandidate;

        var iceServers = [];

        if (moz) {
            iceServers.push({
                url: 'stun:23.21.150.121'
            });

            iceServers.push({
                url: 'stun:stun.services.mozilla.com'
            });
        }

        if (!moz) {
            iceServers.push({
                url: 'stun:stun.l.google.com:19302'
            });

            iceServers.push({
                url: 'stun:stun.anyfirewall.com:3478'
            });
        }

        if (!moz && chromeVersion < 28) {
            iceServers.push({
                url: 'turn:homeo@turn.bistri.com:80',
                credential: 'homeo'
            });
        }

        if (!moz && chromeVersion >= 28) {
            iceServers.push({
                url: 'turn:turn.bistri.com:80',
                credential: 'homeo',
                username: 'homeo'
            });

            iceServers.push({
                url: 'turn:turn.anyfirewall.com:443?transport=tcp',
                credential: 'webrtc',
                username: 'webrtc'
            });
        }

        if (options.iceServers) iceServers = options.iceServers;

        iceServers = {
            iceServers: iceServers
        };

        console.debug('ice-servers', JSON.stringify(iceServers.iceServers, null, '\t'));

        var optional = {
            optional: []
        };

        if (!moz) {
            optional.optional = [{
                DtlsSrtpKeyAgreement: true
            }];
        }

        console.debug('optional-arguments', JSON.stringify(optional.optional, null, '\t'));

        var peer = new PeerConnection(iceServers, optional);

        openOffererChannel();

        peer.onicecandidate = function(event) {
            if (event.candidate)
                options.onICE(event.candidate);
        };

        // attachStream = MediaStream;
        if (options.attachStream) peer.addStream(options.attachStream);

        // attachStreams[0] = audio-stream;
        // attachStreams[1] = video-stream;
        // attachStreams[2] = screen-capturing-stream;
        if (options.attachStreams && options.attachStream.length) {
            var streams = options.attachStreams;
            for (var i = 0; i < streams.length; i++) {
                peer.addStream(streams[i]);
            }
        }

        peer.onaddstream = function(event) {
            var remoteMediaStream = event.stream;

            // onRemoteStreamEnded(MediaStream)
            remoteMediaStream.onended = function() {
                if (options.onRemoteStreamEnded) options.onRemoteStreamEnded(remoteMediaStream);
            };

            // onRemoteStream(MediaStream)
            if (options.onRemoteStream) options.onRemoteStream(remoteMediaStream);

            console.debug('on:add:stream', remoteMediaStream);
        };

        var constraints = options.constraints || {
            optional: [],
            mandatory: {
                OfferToReceiveAudio: true,
                OfferToReceiveVideo: true
            }
        };

        console.debug('sdp-constraints', JSON.stringify(constraints.mandatory, null, '\t'));

        // onOfferSDP(RTCSessionDescription)

        function createOffer() {
            if (!options.onOfferSDP) return;

            peer.createOffer(function(sessionDescription) {
                sessionDescription.sdp = setBandwidth(sessionDescription.sdp);
                peer.setLocalDescription(sessionDescription);
                options.onOfferSDP(sessionDescription);

                console.debug('offer-sdp', sessionDescription.sdp);
            }, onSdpError, constraints);
        }

        // onAnswerSDP(RTCSessionDescription)

        function createAnswer() {
            if (!options.onAnswerSDP) return;

            //options.offerSDP.sdp = addStereo(options.offerSDP.sdp);
            console.debug('offer-sdp', options.offerSDP.sdp);
            peer.setRemoteDescription(new SessionDescription(options.offerSDP), onSdpSuccess, onSdpError);
            peer.createAnswer(function(sessionDescription) {
                sessionDescription.sdp = setBandwidth(sessionDescription.sdp);
                peer.setLocalDescription(sessionDescription);
                options.onAnswerSDP(sessionDescription);
                console.debug('answer-sdp', sessionDescription.sdp);
            }, onSdpError, constraints);
        }

        // if Mozilla Firefox & DataChannel; offer/answer will be created later
        if ((options.onChannelMessage && !moz) || !options.onChannelMessage) {
            createOffer();
            createAnswer();
        }

        // options.bandwidth = { audio: 50, video: 256, data: 30 * 1000 * 1000 }
        var bandwidth = options.bandwidth;

        function setBandwidth(sdp) {
            if (moz || !bandwidth /* || navigator.userAgent.match( /Android|iPhone|iPad|iPod|BlackBerry|IEMobile/i ) */) return sdp;

            // remove existing bandwidth lines
            sdp = sdp.replace( /b=AS([^\r\n]+\r\n)/g , '');

            if (bandwidth.audio) {
                sdp = sdp.replace( /a=mid:audio\r\n/g , 'a=mid:audio\r\nb=AS:' + bandwidth.audio + '\r\n');
            }

            if (bandwidth.video) {
                sdp = sdp.replace( /a=mid:video\r\n/g , 'a=mid:video\r\nb=AS:' + bandwidth.video + '\r\n');
            }

            if (bandwidth.data) {
                sdp = sdp.replace( /a=mid:data\r\n/g , 'a=mid:data\r\nb=AS:' + bandwidth.data + '\r\n');
            }

            return sdp;
        }

        // DataChannel management
        var channel;

        function openOffererChannel() {
            if (!options.onChannelMessage)
                return;

            _openOffererChannel();
        }

        function _openOffererChannel() {
            // protocol: 'text/chat', preset: true, stream: 16
            // maxRetransmits:0 && ordered:false
            var dataChannelDict = { };
            channel = peer.createDataChannel(options.channel || 'sctp-channel', dataChannelDict);
            setChannelEvents();
        }

        function setChannelEvents() {
            channel.onmessage = function(event) {
                if (options.onChannelMessage) options.onChannelMessage(event);
            };

            channel.onopen = function() {
                if (options.onChannelOpened) options.onChannelOpened(channel);
            };
            channel.onclose = function(event) {
                if (options.onChannelClosed) options.onChannelClosed(event);

                console.warn('WebRTC DataChannel closed', event);
            };
            channel.onerror = function(event) {
                if (options.onChannelError) options.onChannelError(event);

                console.error('WebRTC DataChannel error', event);
            };
        }

        if (options.onAnswerSDP && options.onChannelMessage) {
            openAnswererChannel();
        }

        function openAnswererChannel() {
            peer.ondatachannel = function(event) {
                channel = event.channel;
                setChannelEvents();
            };
        }

        // fake:true is also available on chrome under a flag!

        function useless() {
            console.error('Error in fake:true');
        }

        function onSdpSuccess() {
        }

        function onSdpError(e) {
            var message = JSON.stringify(e, null, '\t');

            if (message.indexOf('RTP/SAVPF Expects at least 4 fields') != -1) {
                message = 'It seems that you are trying to interop RTP-datachannels with SCTP. It is not supported!';
            }

            console.error('onSdpError:', message);
        }

        return {
            addAnswerSDP: function(sdp) {
                console.debug('adding answer-sdp', sdp.sdp);
                peer.setRemoteDescription(new SessionDescription(sdp), onSdpSuccess, onSdpError);
            },
            addICE: function(candidate) {
                peer.addIceCandidate(new IceCandidate({
                    sdpMLineIndex: candidate.sdpMLineIndex,
                    candidate: candidate.candidate
                }));

                console.debug('adding-ice', candidate.candidate);
            },

            peer: peer,
            channel: channel,
            sendData: function(message) {
                channel && channel.send(message);
            }
        };
    }

    // getUserMedia
    var video_constraints = {
        mandatory: { },
        optional: []
    };

    function getUserMedia(options) {
        var n = navigator,
            media;
        n.getMedia = n.webkitGetUserMedia || n.mozGetUserMedia;
        n.getMedia(options.constraints || {
                audio: true,
                video: video_constraints
            }, streaming, options.onerror || function(e) {
                console.error(e);
            });

        function streaming(stream) {
            var video = options.video;
            if (video) {
                video[moz ? 'mozSrcObject' : 'src'] = moz ? stream : window.webkitURL.createObjectURL(stream);
                video.play();
            }
            options.onsuccess && options.onsuccess(stream);
            media = stream;
        }

        return media;
    }
  """
end


text:
function tcm_rankingStats(N, lambda, orderings)
	CS_size_by_ranking = (orderings .== N) * collect(0:(N-1))
	avg_CS_size = mean(CS_size_by_ranking)
	weighted_avg_CS_size = dot(lambda, CS_size_by_ranking)
	max_CS_size = maximum(CS_size_by_ranking)

	return avg_CS_size, weighted_avg_CS_size, max_CS_size
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
# Create plots and a LaTeX table for the parameters derived from the photodynamic
# model.

using Statistics
using Printf
using JLD2; using PyPlot
using MCMCDiagnostics

if ~@isdefined(CGS)
  include("../../../src/CGS.jl")
  using Main.CGS
end

# This is run from T1_photodynamics_MCMC_02.ipynb
# (200,000 steps; 50 chains took ~3 days to run):
@load "../../../data/T1_pd_MCMC_run_000_001.jld2"
chain001 = chain
flatchain001 = flatchain
@load "../../../data/T1_pd_MCMC_run_000_002.jld2"
chain002 = chain
flatchain002 = flatchain
@load "../../../data/T1_pd_MCMC_run_000_003.jld2"
chain003 = chain
flatchain003 = flatchain
# (This doesn't have variable contamination factor).

#istart = 15000; iend = 500000
istart = 1001; iend = 200000
flatchain = hcat(flatchain001[:,istart:iend],flatchain002[:,istart:iend],flatchain003[:,istart:iend])
nsamples = size(flatchain)[2]

iburn = 201; ilength = 4000
nwalker = 50
# number of effective samples:
for i=1:19
 neff_tot = 0
 for j=1:nwalker
   neff_tot += effective_sample_size(chain001[i,j,iburn:ilength])
   neff_tot += effective_sample_size(chain002[i,j,iburn:ilength])
   neff_tot += effective_sample_size(chain003[i,j,iburn:ilength])
 end
 println(i," ",neff_tot)
end

# First, make histograms of some parameters:

function histogram(param,nbin)
  p1 = minimum(param)-1e-15; p2 = maximum(param)+1e-15
  pbin_square = [p1]
  hist_square = [0.0]
  pbin = zeros(nbin)
  hist = zeros(nbin)
  psort = sort(param)
  i1 = 1; np = size(param)[1]
  #println(p1," ",psort[i1]," ",p2," ",psort[np])
  for i=1:nbin
    while psort[i1] <= p1+(p2-p1)/nbin*i
      hist[i] += 1.0
      i1 += 1
      if i1 == np
        break
      end
    end
    push!(hist_square,hist[i]); push!(hist_square,hist[i])
    pbin[i] = p1+(p2-p1)/nbin*(i-0.5)
    push!(pbin_square,p1+(p2-p1)/nbin*(i-1))
    push!(pbin_square,p1+(p2-p1)/nbin*i)
    if i1 == np
      break
    end
  end
  push!(hist_square,0.0)
  push!(pbin_square,p2)
  return pbin,hist,pbin_square,hist_square
end

# plot histogram of each parameter:
pname = ["b","c","d","e","f","g","h"]

clf()
nplanet = n-1
rp_mean = zeros(nplanet); rp_sig = zeros(nplanet)
depth_mean = zeros(nplanet); depth_sig = zeros(nplanet)
dur_mean = zeros(nplanet); dur_sig = zeros(nplanet)
tau_mean = zeros(nplanet); tau_sig = zeros(nplanet)
med_b = zeros(nplanet); b_sig = zeros(nplanet)
# Normalized velocity from the dynamical model (edge-on) computed
# at mid-transit for trappist-1 in units of AU/day:
vnorm = [0.10676978639803744, 0.09139150615416661, 0.07710600470245187, 0.06732285130838131, 0.0586281758248469,
         0.05318938396566236, 0.04623982742725911]
rpstring = ""
depthstring = ""
durstring = ""
bstring = ""
taustring = ""
for i=1:nplanet
  # Impact parameter:
  bcurr = sort(flatchain[i+nplanet,:])
  b_bin, b_hist, b_bin_square,b_hist_square = histogram(bcurr,50)
  plot(b_bin,b_hist./maximum(b_hist),label=pname[i],linewidth=2)
  rcurr= flatchain[i,:]
  rp_mean[i] = mean(rcurr)
  rp_sig[i] = std(rcurr)
  depth_mean[i] = mean(flatchain[i,:].^2)
  depth_sig[i] = std(flatchain[i,:].^2)
  med_b[i] = median(bcurr)
  sig_b_low = med_b[i] - bcurr[ceil(Int64,0.158655*nsamples)]
  sig_b_high = bcurr[ceil(Int64,0.841345*nsamples)]-med_b[i]
  # Compute the transit duration:
  fac = flatchain[15,:].^(1//3)*AU/RSUN
  dur = 2sqrt.((1 .+flatchain[i,:]).^2 .- flatchain[i+nplanet,:].^2)./(vnorm[i]*fac)*(24*60)
  dur_mean[i] = mean(dur)
  dur_sig[i] = std(dur)
  tau = (sqrt.((1 .+flatchain[i,:]).^2 .- flatchain[i+nplanet,:].^2) -
          sqrt.((1 .-flatchain[i,:]).^2 .- flatchain[i+nplanet,:].^2))./(vnorm[i]*fac)*(24*60)
  tau_mean[i] = mean(tau)
  tau_sig[i] = std(tau)
  global rpstring = string(rpstring," \$",@sprintf("%6.5f",rp_mean[i])," \\pm ",@sprintf("%6.5f",rp_sig[i]),"\$")
  global depthstring = string(depthstring," \$",@sprintf("%6.4f",depth_mean[i]*100)," \\pm ",@sprintf("%6.4f",depth_sig[i]*100),"\$")
  global durstring = string(durstring," \$",@sprintf("%6.2f",dur_mean[i])," \\pm ",@sprintf("%6.2f",dur_sig[i]),"\$")
  global taustring = string(taustring," \$",@sprintf("%6.3f",tau_mean[i])," \\pm ",@sprintf("%6.3f",tau_sig[i]),"\$")
  if i < nplanet
    global rpstring = string(rpstring," &")
    global depthstring = string(depthstring," &")
    global durstring = string(durstring," &")
    global taustring = string(taustring," &")
  else
    global rpstring = string(rpstring," \\cr")
    global depthstring = string(depthstring," \\cr")
    global durstring = string(durstring," \\cr")
    global taustring = string(taustring," \\cr")
  end
#  global bstring = string(bstring," \$",@sprintf("%6.4f",b_mean[i])," \\pm ",@sprintf("%6.4f",b_sig[i]),"\$")
  global bstring = string(bstring," \$",@sprintf("%6.3f",med_b[i]),"_{- ",@sprintf("%6.3f",sig_b_low),"}^{+",@sprintf("%6.3f",sig_b_high),"}\$")
  if i < nplanet
    global bstring = string(bstring," &")
  else
    global bstring = string(bstring," \\cr")
  end
end
legend()
axis([0.0,0.5,0,1.1])
#xlabel("Radius ratio"); ylabel("Probability")
xlabel("Impact parameter",fontsize=15); ylabel("Probability",fontsize=15)
savefig("../impact_parameter_noprior.pdf",bbox_inches="tight")

# Now, create a table of values:

# Compute 68% confidence interval and median of density:
density = sort(flatchain[15,:])
med_dens = median(density); sig_low = med_dens - density[ceil(Int64,0.158655*nsamples)]
sig_high = density[ceil(Int64,0.841345*nsamples)]-med_dens
#println("\$\\rho_*/\\rho_\\odot\$ & \$",@sprintf("%6.4f",med_dens),"\_\{- ",@sprintf("%6.4f",sig_low)"\}\^\{+",@sprintf("%6.4f",sig_high),"\$ & & & & & &\\cr")

# Also derive a/R_* and inclinations of the orbits, and add to the table. [ ]

q11 = flatchain[16,:]
q21 = flatchain[17,:]
q12 = flatchain[18,:]
q22 = flatchain[19,:]
#println("Parameter & \$\\rho_*/\\rho_\\odot\$ & \$q_{1,Ch 1}\$ & \$q_{2,Ch 1}\$ & \$ q_{1,Ch 2}\$ & \$q_{2,Ch 2}\$ & \$ & &\\cr")
#println("Value & \$",@sprintf("%6.4f",med_dens),"_{- ",@sprintf("%6.4f",sig_low),"}^{+",@sprintf("%6.4f",sig_high),"}\$ & \$",@sprintf("%6.4f",mean(q11))," \\pm ",@sprintf("%6.4f",std(q11))," \$ & \$",@sprintf("%6.4f",mean(q21))," \\pm ",@sprintf("%6.4f",std(q21))," \$ & \$",@sprintf("%6.4f",mean(q12))," \\pm ",@sprintf("%6.4f",std(q12))," \$ & \$",@sprintf("%6.4f",mean(q22))," \\pm ",@sprintf("%6.4f",std(q22))," \$ & &\\cr")
#println("\\hline")
#println(" Planet & b & c & d & e & f & g \\cr")
#println("\$R_p/R_*\$ & ",rpstring)
#println("\$b/R_*\$ & ",bstring)

#read(stdin,Char)
clf()
@load "../../../data/T1_density_incprior.jld2"

rho_bin, rho_hist, rho_bin_square,rho_hist_square = histogram(density,100)
plot(rho_bin,rho_hist./maximum(rho_hist),linewidth=2,label="No prior")
plot(rho_bin_incprior[1:end-2],rho_hist_incprior[1:end-2]./maximum(rho_hist_incprior),linestyle="--",linewidth=2,label="Inclination prior")
xlabel(L"$\rho_*/\rho_\odot$",fontsize=15)
ylabel("Probability",fontsize=15)
legend(fontsize=15)
axis([47,56,0,1.2])
savefig("../stellar_density_noprior.pdf",bbox_inches="tight")

#plot(log.(rho_bin_square./54.0),log.(rho_hist_square))

#loglog(rho_bin_square./54,(rho_bin_square/54.0).^61 .* exp(12)./(exp.(rho_bin_square/54) .-1))

#loglog(rho_bin_square./54,(rho_bin_square ./54).^3 ./(exp.(rho_bin_square ./54).^3 .-1))

#loglog(rho_bin_square./54,(rho_bin_square ./54).^60 .*exp(11))
#loglog(rho_bin_square./54,exp(-(rho_bin_square./54))*1e5)

#rho_bin_min = copy(rho_bin_square)
#rhoc = 53.7
#plot(rho_bin_square ./ rhoc,rho_hist_square)
#rho_bin_min[rho_bin_min .< rhoc] .= rhoc
#semilogy(rho_bin_square./54,exp.((rho_bin_square .-45) .*1.25))
#plot(rho_bin_square ./ 54,exp.(-((rho_bin_min ./ 54) .-1).^2 .* 0.5 .*1e4)*6e4)
# This doesn't do too bad of a job in reproducing the distribution:
#plot(rho_bin_square./rhoc,exp.((rho_bin_square .-45) .*1.25).*exp.(-((rho_bin_min ./ rhoc) .-1).^2 .* 0.5 .*1e4))

# Plot histograms of a/R_* for all planets:
#read(stdin,Char)
#clf()

P = [1.510826,2.421938,4.049218,6.101013,9.207541,12.352445,18.772863]  # Orbit periods

aonr0 = (GRAV*MSUN*(24*3600)^2/(4pi^2*RSUN^3))^(1//3)

med_dens = median(density); dens_low = density[ceil(Int64,0.158655*nsamples)]; dens_high = density[ceil(Int64,0.841345*nsamples)]
aonr_string = "\$a/R_*\$ & \$"
for i=1:7
  aonr_range = [dens_low,med_dens,dens_high].^(1//3) .* aonr0 .* P[i]^(2//3)
  println("aon_range: ",aonr_range)
  global aonr_string = string(aonr_string,@sprintf("%5.3f",aonr_range[2]),"_{- ",@sprintf("%5.3f",aonr_range[2]-aonr_range[1]),"}^{+",@sprintf("%5.3f",aonr_range[3]-aonr_range[2]),"}\$")
  if i < 7
    aonr_string = string(aonr_string," & \$")
  end
end
aonr_string = string(aonr_string," \\cr")

#  aonr = density.^(1//3) .* (aonr0 * P[i]^(2//3))
#  aonr_bin,aonr_hist,aonr_bin_sq,aonr_hist_sq = histogram(aonr,50)
#  plot(aonr_bin_sq,aonr_hist_sq,label=pname[i])
#end
#xlabel(L"$a/R_*$")
#ylabel(L"Probability")
#legend()

# Derive the inclinations of the planets:
#read(stdin,Char)
clf()
nplanet = n-1
inc_mean = zeros(nplanet); inc_sig = zeros(nplanet)
incstring = " \$ I \$(\$^\\circ\$) & "
cname = ["C0","C1","C2","C3","C4","C5","C6","C7"]
nb = [30,30,30,40,40,40,40,40]
for i=1:nplanet
  density = sort(flatchain[15,:])
  aonr = density.^(1//3) .* (aonr0 * P[i]^(2//3))
  bcurr = flatchain[i+nplanet,:]
  inc_planet = acos.(bcurr ./ aonr) .* (180/pi)
  inc_bin, inc_hist, inc_bin_square,inc_hist_square = histogram(inc_planet,nb[i])
#  plot(inc_bin_square,inc_hist_square./maximum(inc_hist_square),label=pname[i],color=cname[i])
  plot([inc_bin;reverse(180 .-inc_bin)],[inc_hist;reverse(inc_hist)]./maximum(inc_hist),label=pname[i],color=cname[i],linewidth=2)
#  plot(180 .-inc_bin_square,inc_hist_square./maximum(inc_hist_square),color=cname[i])
  inc_mean[i] = mean(inc_planet)
  inc_sig[i] = std(inc_planet)
  global incstring = string(incstring," \$",@sprintf("%6.3f",inc_mean[i])," \\pm ",@sprintf("%6.3f",inc_sig[i]),"\$")
  if i < nplanet
    global incstring = string(incstring," & ")
  else
    global incstring = string(incstring," \\cr")
  end
end
legend(fontsize=15)
axis([89.25,90.75,0,1.1])
#xlabel("Radius ratio"); ylabel("Probability")
xlabel("Inclination [deg]",fontsize=15); ylabel("Probability",fontsize=15)
tight_layout()
savefig("../inclination_noprior.pdf",bbox_inches="tight")
#read(stdin,Char)

# Limb-darkening:
clf()
u21 = sqrt.(q11).*(1 .-2q21)
u11 = 2*sqrt.(q11).*q21
u22 = sqrt.(q12).*(1 .-2q22)
u12 = 2*sqrt.(q12).*q22
#plot(u11,u21,color="r",".",alpha=0.005)
#plot(u12,u22,color="g",".",alpha=0.005)
errorbar([0.1633],[0.2549],xerr=[0.0364],yerr=[0.057],color="k",linewidth=3)
plot([0.1633],[0.2549],color="r","o",label="Ch 1")
errorbar([0.1442],[0.2173],xerr=[0.0324],yerr=[0.0482],color="k",linewidth=3)
plot([0.1442],[0.2173],color="g","o",label="Ch 2")
xlabel(L"$u_1$",fontsize=15); ylabel(L"$u_2$",fontsize=15)
legend()
include("../../../src/contour_limb.jl")
axis([u1i,u1f,u2i,u2f])
tight_layout()
savefig("../limb_darkening_nouprior.pdf",bbox_inches="tight")
rho_sun = MSUN/(4pi/3*RSUN^3)

println("Parameter & \$\\rho_*/\\rho_\\odot\$ & \$q_{1,Ch 1}\$ & \$q_{2,Ch 1}\$ & \$ q_{1,Ch 2}\$ & \$q_{2,Ch 2}\$ &  &\\cr")
println("Value & \$",@sprintf("%6.4f",med_dens),"_{- ",@sprintf("%6.4f",sig_low),"}^{+",@sprintf("%6.4f",sig_high),"}\$ & \$",@sprintf("%6.4f",mean(q11))," \\pm ",@sprintf("%6.4f",std(q11))," \$ & \$",@sprintf("%6.4f",mean(q21))," \\pm ",@sprintf("%6.4f",std(q21))," \$ & \$",@sprintf("%6.4f",mean(q12))," \\pm ",@sprintf("%6.4f",std(q12))," \$ & \$",@sprintf("%6.4f",mean(q22))," \\pm ",@sprintf("%6.4f",std(q22))," \$ & &\\cr")
#println("Parameter & \$ \\rho_*\$ \[g/cm\$\^2\$\] & \$u_\\mathrm{1,Ch 1}\$ & \$u_\\mathrm{2,Ch 1}\$ & \$u_\\mathrm{1,Ch 2}\$ & \$u_\\mathrm{2,Ch 2}\$ &  &\\cr")
println("Parameter & \$ \\rho_*\$ [g / cm\$^3\$] & \$u_\\mathrm{1,Ch 1}\$ & \$u_\\mathrm{2,Ch 1}\$ & \$u_\\mathrm{1,Ch 2}\$ & \$u_\\mathrm{2,Ch 2}\$ &  & \\cr")

println("Value & \$",@sprintf("%6.4f",med_dens*rho_sun),"_{- ",@sprintf("%6.4f",sig_low*rho_sun),"}^{+",@sprintf("%6.4f",sig_high*rho_sun),"}\$ & \$",@sprintf("%6.4f",mean(u11))," \\pm ",@sprintf("%6.4f",std(u11))," \$ & \$",@sprintf("%6.4f",mean(u21))," \\pm ",@sprintf("%6.4f",std(u21))," \$ & \$",@sprintf("%6.4f",mean(u12))," \\pm ",@sprintf("%6.4f",std(u12))," \$ & \$",@sprintf("%6.4f",mean(u22))," \\pm ",@sprintf("%6.4f",std(u22))," \$ & & \\cr")
println("\\hline")
println(" Planet & b & c & d & e & f & g \\cr")
println("\$R_p/R_*\$ & ",rpstring)
println("Depth [\\%] & ",depthstring)
println("T [min] & ",durstring)
println("\$\\tau\$ [min] & ",taustring)
println("\$b/R_*\$ & ",bstring)
println(aonr_string)
println(incstring)


text:
# domain, distribution, and ergodic manager
d = Domain([1,1,1], 100)
means = [[.2,.2,.2], [.8,.8,.2], [.5,.5,.8]]
covs = [0.01*eye(3), 0.01*eye(3), .01*eye(3)]
phi = gaussian(d, means, covs)
K = 5
em = ErgodicManagerR3(d, phi, K)

# trajectory params
x0 = [0.49, 0.01, 0.01]
dt = 0.5
N = 80
tm = TrajectoryManager(x0, dt, N, ConstantInitializer([0.0,0.0,0.0]))
dynamics!(tm, SingleIntegrator(3,dt))
tm.descender = ArmijoLineSearch(1,1e-4)

# trajectory generation and plotting
xd,ud = pto_trajectory(em, tm, dd_crit=1e-4, max_iters=1000)