text:
function intraperiod_segmentation(data_merged::ClustDataMerged;
                n_seg::Int=24,
                iterations::Int=300,
                norm_scope::String="full")
  #For easy access
  K=data_merged.K
  T=data_merged.T
  data=data_merged.data
  #Prepare matrices
  data_seg=zeros(n_seg*length(data_merged.data_type),K)
  deltas_seg=zeros(n_seg,K)
  #Loop over each period
  for k in 1:K
    #Take single period and reshape it to have t (1:T) as columns and types (solar-germany, wind-germany...) as rows
    period=permutedims(reshape(data[:,k],(T,length(data_merged.data_type))))
    #Run hierarchical clustering
    centers,weights,clustids,cost,iter=run_clust_segmentation(period;n_seg=n_seg,iterations=iterations,norm_scope=norm_scope)
    #Assign values back to matrices to match types*n_seg x K
    data_seg[:,k]=reshape(permutedims(centers),(size(data_seg,1),1))
    #Assign values back to matrices to match n_seg x K
    deltas_seg[:,k]=weights
  end
  return ClustDataMerged(data_merged.region,data_merged.years,K,n_seg,data_seg,data_merged.data_type,data_merged.weights,data_merged.k_ids;deltas=deltas_seg,)
end

function run_clust_segmentation(period::Array{Float64,2};
                n_seg::Int=24,
                iterations::Int=300,
                norm_scope::String="full")
  norm_period, typely_mean, typely_sdv=z_normalize(period;scope=norm_scope)
  #x,weights,clustids,x,iter= run_clust_hierarchical(norm_period,n_seg,iterations)
  data=norm_period
  clustids=run_clust_hierarchical_partitional(data::Array, n_seg::Int64)
  weights = calc_weights(clustids,n_seg)

  centers_norm = calc_centroids(norm_period,clustids)
  cost = calc_SSE(norm_period,centers_norm,clustids)
  centers = undo_z_normalize(centers_norm,typely_mean,typely_sdv;idx=clustids)
  return centers,weights,clustids,cost,1
end

function get_clustids(ends::Array{Int64,1})
  clustids=collect(1:size(data,2))
  j=1
  for i in 1:size(data,2)
    clustids[i]=j
    if i in ends
      j+=1
    end
  end
  return clustids
end

"""
function run_clust_hierarchical_partitional(data::Array, n_seg::Int64)
  Usees provided data and number of segments to aggregate them together

Beginning implementation Bahl
function run_clust_hierarchical_partitional(data::Array,
                                            n_seg::Int64)
  ends_all=LHCoptim((size(data,2)-1),(n_seg-1),1000)[1]
  distance_runs=[]
  clustids_runs=[]
  for e in 1:size(ends_all,1)
    ends=ends_all[e,:]
    for n in 1:100
      global clustids=get_clustids(ends)
      #Calculate mean of data: The number of columns is kept the same, mean is calculated for aggregated columns and the same in all with same clustid
      global means_data=get_mean_data(data,clustids)
      ends_prev=ends
      for i in 1:length(ends)
        #
        if sum(abs.(data[:,ends[i]]-means_data[:,ends[i]]))>sum(abs.(data[:,ends[i]]-means_data[:,min(ends[i]+1,size(data,2))]))
          global ends[i]+=1
          #break
        end
      end
      if ends_prev==ends
        #break
      end
    end
    distance=sum(pairwise(SqEuclidean(),means_data,data))
    push!(distance_runs,distance)
    push!(clustids_runs,clustids)
  end
  #return clustids_runs[findfirst(distance_runs.==minimum(distance_runs))]
  return clustids_runs, distance_runs
end
"""

"""
function run_clust_hierarchical_partitional(data::Array, n_seg::Int64)
  Usees provided data and number of segments to aggregate them together
"""
function run_clust_hierarchical_partitional(data::Array,
                                            n_seg::Int64)
  _dist= SqEuclidean()
  #Assign each timeperiod it's own cluster
  clustids=collect(1:size(data,2))
  #While aggregation not finished, aggregate
  #Calculate the sq distance
  d_mat=pairwise(_dist,data)
  while clustids[end]>n_seg
    #Calculate mean of data: The number of columns is kept the same, mean is calculated for aggregated columns and the same in all with same clustid
    #Initially no index is selected and distance is Inf
    NNnext=0
    NNdist=Inf
    # loop through the sq distance matrix to check:
    for i=1:(clustids[end]-1)
      # if the distance between this index [i] and it's neighbor [i+1] is lower than the minimum found so far
       #distance=sum(d_mat[findall(clustids.==i),findall(clustids.==i+1)])
      clustids_test=deepcopy(clustids)
      merge_clustids!(clustids_test,findlast(clustids.==i))
      distance=calc_SSE(data,clustids_test)
      #println(distance)
      if distance < NNdist
        #Save this index and the distance
        NNnext=findlast(clustids.==i)
        NNdist=distance
      end
    end
    # Aggregate the clustids that were closest to each other
    merge_clustids!(clustids,NNnext)
  end
  return clustids
end

"""
function merge_clustids!(clustids::Array{Int64,1},index::Int64)
Calculate the new clustids by merging the cluster of the index provided with the cluster of index+1
"""
function merge_clustids!(clustids::Array{Int64,1},index::Int64)
  clustids[index+1]=clustids[index]
  clustids[index+2:end].-=1
end

"""
function get_mean_data(data::Array, clustids::Array{Int64,1})
Calculate mean of data: The number of columns is kept the same, mean is calculated for aggregated columns and the same in all with same clustid
"""
function get_mean_data(data::Array,
                    clustids::Array{Int64,1})
  mean_data=zeros(size(data))
  for i in 1:size(data,2)
    mean_data[:,i]=mean(data[:,findall(clustids.==clustids[i])], dims=2)
  end
  return mean_data
end


text:
using BigArrays.Indexes

using Test

@testset "test index transformation" begin
    str = "2968-3480_1776-2288_16912-17424"
    range = [2969:3480, 1777:2288, 16913:17424]
    @test range == Indexes.string2unit_range( str )
    @test str == Indexes.unit_range2string( range )
end

ext: jl


text:
using Base.Test
import COBRAUtils

@testset "AV, 1" begin
    model = COBRAUtils.AVCobraModel(0.5, 0.)
    status, μ, v = COBRAUtils.optimize(model)
    @test status == :Optimal
    @test μ ≈ 17.2
    @test v ≈ [0.5, -0.1, 0.9, 17.2]

    model = COBRAUtils.reactionsIrrev(model)
    status, μ, v = COBRAUtils.optimize(model)
    @test status == :Optimal
    @test μ ≈ 17.2
    @test v ≈ [0.5, 0., 0.9, 17.2, 0.0, 0.1, 0.0, 0.0]

    model = COBRAUtils.addFluxCost(model, ones(length(model.rxns)))
    status, μ, v = COBRAUtils.optimize(model)
    @test status == :Optimal
    @test sum(v[1 : end - 1]) ≈ v[end] ≈ 1
    @test μ ≈ 0.9268292682926829
end

@testset "AV fva" begin
    model = COBRAUtils.AVCobraModel(0.5, 0.)
    fmin, fmax = COBRAUtils.fva(model)
    @test fmin ≈ [0.0, -1.0, 0.0, 0.0]
    @test fmax ≈ [0.5, 0.0, 0.9, 17.2]
end

@testset "AV find" begin
    model = COBRAUtils.AVCobraModel(0.5, 0.)
    @test COBRAUtils.metabolite_reactions(model, "atp") == [1,3,4]
    @test COBRAUtils.metabolite_reactions(model, "pyr") == [1,2,3]
end

@testset "Bounds" begin
    model = COBRAUtils.AVCobraModel(0.5, 0.)
    for i = 1 : length(model.rxns)
        @test COBRAUtils.lower_bound(model, model.rxns[i]) == model.lb[i]
        @test COBRAUtils.upper_bound(model, model.rxns[i]) == model.ub[i]
    end
end


text:
using InstrumentConfig

const EXAMPLE_FILE = "https://raw.githubusercontent.com/Orchard-Ultrasound-Innovation/TcpInstruments.jl/master/.tcp_instruments.yml"

const tcp_config = InstrumentConfig.Config(
    ".tcp_instruments.yml",
    @__MODULE__
)

function get_config()
    return InstrumentConfig.get_config(tcp_config)
end

function create_config(;dir=homedir())
    InstrumentConfig.create_config(tcp_config; dir=dir)
end

function edit_config()
    InstrumentConfig.edit_config(tcp_config)
end

function load_config()
    InstrumentConfig.load_config(tcp_config)
    create_aliases(tcp_config)
end

function create_aliases(config; ignore=[])
    for (device, data) in config.config
        if device in ignore
            continue
        end
        device_type = nothing
        try
            device_type = eval(Symbol(device))
        catch e
            error("""
            $(config.loaded_file) contains device of name:
                $device

            This is not a valid device.

            For a list of available devices use `help> Instrument`
            """)
        end
        !(data isa Dict) && continue
        alias = get(data, "alias", "")
        isempty(alias) && continue
        @eval global const $(Symbol(alias)) = $(device_type)
        @eval export $(Symbol(alias))
        alias_print("$alias = $device")
    end
end


text:
export classify_matches, set_matches, set_symbol_matches

# set equality is checked using a linear time tree isomorphism algorithm
#  https://logic.pdmi.ras.ru/~smal/files/smal_jass08.pdf

function lexicographic_sort!(strings::Vector{String})
    # TODO: implement linear time string sort from https://doc.lagout.org/science/0_Computer%20Science/2_Algorithms/The%20Design%20and%20Analysis%20of%20Computer%20Algorithms%20%5BAho,%20Hopcroft%20&%20Ullman%201974-01-11%5D.pdf page 90

    sort!(strings)
end

get_canonical_name(sym::LogicalSymbol) = "_"
function get_canonical_name(node::ExtensionalSet)
    if node == ∅
        return "10"
    end
    subnames = [get_canonical_name(child) for child ∈ elements(node)]
    lexicographic_sort!(subnames)
    string("1", join(subnames), "0")
end

# SET MATCHING
@enum MatchClassification EXACT_MATCH PARTIAL_MATCH NO_MATCHES
function classify_matches(matches::Dict{LogicalSymbol, Set{AbstractExpression}})
    match_lengths = length.(collect(values(matches)))
    if any(match_lengths .== 0)
        return NO_MATCHES
    end
    if all(match_lengths .== 1)
        return EXACT_MATCH
    end
    return PARTIAL_MATCH
end
function classify_matches(pattern::ExtensionalSet, haystack::ExtensionalSet)
    classify_matches(set_matches(pattern, haystack))
end

function reduce_matches!(unreduced_matches::Dict)
    # TODO: consider this scenario:
    # a = [1, 2]
    # b = [1, 2]
    # c = [1, 2, 3]
    # We cannot make any conclusions about a or b, but we _can_ conclude that c must not be either 1 or 2, considering
    # if it was the problem would become unsatisfiable between a and b. This kind of higher-order logical resolution
    # probably will require a neat algorithm of its own, which constitutes this problem as another issue (#21).

    continue_running = true

    completed_symbols = Set{AbstractExpression}()
    while continue_running
        continue_running = false
        for (sym, matches) ∈ unreduced_matches
            if length(matches) == 1 && sym ∉ completed_symbols
                ideal_match = first(matches)
                push!(completed_symbols, sym)
                # this could be simplified with a "graph inversion" but usually dictionary sizes are low so this doesn't cost too much
                # NOTE TO FUTURE ME: graph inversion might also be a clever way to solve #21 as well...
                for (sym2, matches2) ∈ unreduced_matches
                    if sym == sym2
                        continue
                    end
                    if ideal_match ∈ matches2
                        delete!(matches2, ideal_match)
                    end
                end

                continue_running = true
            end
        end
    end

    unreduced_matches  # but now its reduced from the while loop :)
end

# TODO: Move these elsewhere or switch fully to ExtensionalSet from Set
ES = Union{ExtensionalSet, Set}
elements(s::Set) = s
cardinality(s::Set) = length(s)
variables(s::Set) = cardinality(s) > 0 ? reduce(∪, variables.(collect(s))) : Set{LogicalSymbol}()

# TODO: fix the method signatures to be less messy
#  1. Reverse pattern and haystack to match up with other matching algorithm signatures
#  2. Make sure type specificity is as tight as possible to prevent contradictory type combinations from being valid inputs
function set_matches(pattern::ES, haystack::ES; reduce=true, strict=true)
    matching_subpatterns = Dict(sub => Set{AbstractExpression}() for sub ∈ elements(pattern))

    if strict && cardinality(pattern) != cardinality(haystack)
        return Dict()
    end

    if pattern == ∅ && haystack == ∅
        return Dict(pattern => haystack)
    end

    for subpattern ∈ elements(pattern)
        submatches = [haystack_el => set_matches(subpattern, haystack_el) for haystack_el ∈ elements(haystack)]
        matching_elements = map(x -> x.first, filter(x -> length(x.second) > 0, submatches))
        for matching_el ∈ matching_elements
            push!(matching_subpatterns[subpattern], matching_el)
        end
    end

    # PHYSICAL ROOT MATCHING
    # if ∅ has matches, remove ∅ from other potential matched symbols
    if haskey(matching_subpatterns, ∅) && length(matching_subpatterns[∅]) > 0
        for (subpatt, opts) ∈ matching_subpatterns
            if ∅ ∈ opts && subpatt != ∅
                deleteat!(opts, findall(x->x==∅, opts))
            end
        end
    end
    # END PHYSICAL ROOT MATCHING

    # we can perform naive simplification
    if any(length.(collect(values(matching_subpatterns))) .== 1)
        reduce_matches!(matching_subpatterns)
    end

    if reduce
        constrained_output_elimination!(matching_subpatterns)
    end

    matching_subpatterns
end
set_matches(sym::LogicalSymbol, haystack::Union{ES, AbstractExpression}; reduce=false) = Dict(sym => Set([haystack]))
# TODO: think about generalizing set_matches such that it integrates with `find_matches` since they fundamentally do something similar
function set_matches(pattern::LogicalExpression, haystack::AbstractExpression)
    expr_matches = Dict{LogicalSymbol, AbstractExpression}()
    if find_matches!(haystack, pattern, expr_matches)
        return Dict(k => Set([v]) for (k, v) ∈ expr_matches)
    end
    Dict()
end

function set_symbol_matches(pattern::Union{ES, LogicalSymbol}, haystack::ES; symbol_matches=nothing, reduce=true, strict=true)
    if isnothing(symbol_matches)
        symbol_matches = Dict{LogicalSymbol, Set{AbstractExpression}}(sym => Set{AbstractExpression}() for sym ∈ variables(pattern))
    end

    current_set_matches = set_matches(pattern, haystack; reduce=false)
    for (expr, matches) ∈ current_set_matches
        if expr isa LogicalSymbol
            union!(symbol_matches[expr], Set(matches))
        else
            set_symbol_matches.(flat_repeat(expr, length(matches)), matches; symbol_matches, reduce=false, strict=strict)
        end
    end

    if reduce
        constrained_output_elimination!(symbol_matches)
    end

    symbol_matches
end
function set_symbol_matches(pattern::LogicalExpression, haystack::AbstractExpression; symbol_matches=nothing, kwargs...)
    my_matches = set_matches(pattern, haystack)
    for (sym, matches) ∈ my_matches
        union!(symbol_matches[sym], matches)
    end
end


text:
cbi = function (p, i_exp)
    ode_data = ode_data_list[i_exp, :, :]
    _ts = 0:ts[end]/nplot:ts[end]
    pred = predict_neuralode(u0, p, i_exp, nothing, false)
    gold = predict_neuralode(u0, p_gold, i_exp, nothing, false)
    l_plt = []
    for i in 1:minimum([10, ns])
        plt = Plots.scatter(ts, ode_data[i,:], label="Data");
        plot!(plt, _ts, pred[i,:], label="Prediction");
        plot!(plt, _ts, gold[i,:], label="Ground truth");
        ylabel!(plt, "x$i")
        xlabel!(plt, "Time")
        push!(l_plt, plt)
    end
    plt_all = plot(l_plt..., legend=false, size=(1000, 1000));
    png(plt_all, string(fig_path, "/conditions/i_exp_", i_exp))
    return false
end

cbi! = function (p, i_exp, plt)
    ode_data = ode_data_list[i_exp, :, :]
    _ts = 0:ts[end]/nplot:ts[end]
    pred = predict_neuralode(u0, p, i_exp, nothing, false)
    gold = predict_neuralode(u0, p_gold, i_exp, nothing, false)
    ylabel!(plt, "Protein levels")
    xlabel!(plt, "Time")
    plot!(plt, _ts, gold[1,:], line=:dash, label="Ground truth", color=palette(:tab10)[1]);
    Plots.scatter!(plt, ts, ode_data[1,:], label="Data", color=palette(:tab10)[1]);
    plot!(plt, _ts, pred[1,:], label="Prediction", color=palette(:tab10)[1]);
    for i in 2:minimum([5, ns])
        plot!(plt, _ts, gold[i,:], line=:dash, label=nothing, color=palette(:tab10)[i]);
        Plots.scatter!(plt, ts, ode_data[i,:], label=nothing, color=palette(:tab10)[i]);
        plot!(plt, _ts, pred[i,:], label=nothing, color=palette(:tab10)[i]);
    end
    png(plt, string(fig_path, "/trajectories_cond_", i_exp, ".png"))
end

cbrd = function (p)
    ode_data = ode_data_list[:,:,1]
    pred_data = zeros(size(ode_data_list)[1:2])
    for i_exp in 1:n_exp
        pred = predict_neuralode(u0, p, i_exp, nothing, false)
        pred_data[i_exp,:] = pred[:,end]
    end
    plt_train = Plots.scatter([ode_data[1:n_exp_train,:]...],[pred_data[1:n_exp_train,:]...])
    xlabel!("experimental data Training set")
    ylabel!("predicted data")
    plt_val = Plots.scatter([ode_data[(n_exp_train+1):(n_exp_train+n_exp_val),:]...],[pred_data[(n_exp_train+1):(n_exp_train+n_exp_val),:]...])
    xlabel!("experimental data Validation set")
    ylabel!("predicted data")
    plt_test = Plots.scatter([ode_data[(n_exp-n_exp_test+1):end,:]...],[pred_data[(n_exp-n_exp_test+1):end,:]...])
    xlabel!("experimental data Test set")
    ylabel!("predicted data")
    plt_p = plot([plt_train, plt_val, plt_test]..., legend=false, layout = grid(1, 3), size=(1300, 400))
    png(plt_p, string(fig_path, "/predicted_response_correlation_with_real_data.png"))
end

cbp = function (p, iter)
    plt_alpha = Plots.scatter(p_gold[:,1],p[:,1])
    xlabel!("ground truth alpha")
    ylabel!("inferred alpha")
    plt_w = Plots.scatter([p_gold[:,2:end]...],[p[:,2:end]...])
    xlabel!("ground truth w")
    ylabel!("inferred w")
    plt_p = plot([plt_alpha, plt_w]..., legend=false, layout = grid(1, 2), size=(850, 400))
    png(plt_p, string(fig_path, "/p_inference_iter", iter))
end

cbloss = function ()
    plt_loss = plot(l_loss_train, yscale=:log10, label="train")
    plot!(plt_loss, l_loss_val, yscale=:log10, label="val")
    plot!(plt_loss, l_loss_test, yscale=:log10, label="test")
    plt_grad = plot(l_grad, yscale=:log10, label="grad_norm")
    plt_p = plot([l_loss_network[i][1] for i in 1:length(l_loss_network)], label="alpha")
    plot!(plt_p, [l_loss_network[i][2] for i in 1:length(l_loss_network)], label="w")
    xlabel!(plt_loss, "Epoch")
    xlabel!(plt_grad, "Epoch")
    xlabel!(plt_p, "Epoch")
    ylabel!(plt_loss, "Loss")
    ylabel!(plt_grad, "Gradient Norm")
    ylabel!(plt_p, "Parameter correlation")
    # ylims!(plt_loss, (-Inf, 1))
    plt_all = plot([plt_loss, plt_grad, plt_p]..., legend=:top, layout = grid(3, 1), size=(600, 900))
    png(plt_all, string(fig_path, "/loss_grad"))
end


text:
#=Copyright 2018 Vidrio Technologies, LLC

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

  http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

module ScanImageTiffReader

import Base: length, size
export Context, data, description, length, metadata, pxtype, size

_release = "1.3"
_libname = begin
    external = joinpath(dirname(dirname(@__FILE__)), "external")
    if Sys.iswindows()
        lp = joinpath(external, "ScanImageTiffReader-$_release-win64", "lib")
    elseif Sys.isapple()
        lp = joinpath(external, "ScanImageTiffReader-$_release-Darwin", "lib")
    elseif Sys.islinux()
        lp = joinpath(external, "ScanImageTiffReader-$_release-Linux", "lib")
    else
        error("you done goofed")
    end

    joinpath(lp, Sys.iswindows() ? "ScanImageTiffReaderAPI.dll" : "libScanImageTiffReaderAPI.so")
end
_typemap = [UInt8, UInt16, UInt32, UInt64, Int8, Int16, Int32, Int64, Float32, Float64]

"""
    Context(handle::Ptr{Cvoid}, log::Ptr{UInt8})
"""
struct Context
    handle::Ptr{Cvoid}
    log::Ptr{UInt8}
end

struct strides_array
    strides_0::Int64
    strides_1::Int64
    strides_2::Int64
    strides_3::Int64
    strides_4::Int64
    strides_5::Int64
    strides_6::Int64
    strides_7::Int64
    strides_8::Int64
    strides_9::Int64
    strides_10::Int64
end

struct dims_array
    dims_0::UInt64
    dims_1::UInt64
    dims_2::UInt64
    dims_3::UInt64
    dims_4::UInt64
    dims_5::UInt64
    dims_6::UInt64
    dims_7::UInt64
    dims_8::UInt64
    dims_9::UInt64
end

struct Size
    ndim::UInt32
    typeid::Int32
    strides::strides_array
    dims::dims_array
end

"""
    open(func::Function, filename::AbstractString, args...)

Open a ScanImage TIFF file `filename` for reading and apply `func` to it
(with optional arguments).

# Examples
```jldoctest
ScanImageTiffReader.open(mytif) do io
    length(io)
end
# output
10
```
"""
function open(func::Function, filename::AbstractString, args...)
    handle = @eval ccall((:ScanImageTiffReader_Open, $(_libname)), Context,
                         (Ptr{UInt8}, ), $(filename))

    if handle.log != C_NULL
        throw(ErrorException(unsafe_string(handle.log)))
    end

    try
        func(handle, args...)
    finally
        @eval ccall((:ScanImageTiffReader_Close, $(_libname)), Cvoid,
                    (Ptr{Context}, ), $(Ref(handle)))
    end
end

"""
    size(ctx::Context)

Return the shape of the data in the TIFF file.

# Examples
```jldoctest
ScanImageTiffReader.open(mytif) do io
    size(io)
end
# output
(512, 512, 10)
```
"""
function size(ctx::Context)
    s = @eval ccall((:ScanImageTiffReader_GetShape, $(_libname)), Size,
                    (Ptr{Context}, ), $(Ref(ctx)));

    if ctx.log != C_NULL
        throw(ErrorException(unsafe_string(ctx.log)))
    end

    Int64.(getfield.(Ref(s.dims), fieldnames(dims_array)[1:s.ndim]))
end

"""
    pxtype(ctx::Context)

Return the type of the data in the TIFF file.

# Examples
```jldoctest
ScanImageTiffReader.open(mytif) do io
    pxtype(io)
end
# output
Int16
```
"""
function pxtype(ctx::Context)
    s = @eval ccall((:ScanImageTiffReader_GetShape,$(_libname)), Size,
                    (Ptr{Context}, ), $(Ref(ctx)));

    if ctx.log!=C_NULL
        throw(ErrorException(unsafe_string(ctx.log)))
    end

    _typemap[s.typeid + 1] # julia is 1-based
end

"""
    data(ctx::Context)

Return an n-dimensional array containing the image stack.

# Examples
```jldoctest; output=false
using ScanImageTiffReader
d = ScanImageTiffReader.open(mytif) do io
    data(io)
end
d[1:10, 1:10, 1]
# output
10×10 Array{Int16,2}:
  9844   8752  10596  10153   8711   9879   8827   9004   9684   8604
 10038   9315   8696   9950   9061   8582   9013  10098  10364   9551
 10618  10610  10536   9716  10187   8813  10033  10206   8567  10299
  9535   9417   9435  10227   9336   9300   9822   9799   9511   8991
  9190  10558  10147  10129   9564   9112   9622   9074   8570   9614
  9152   9148  10318  10175   8917  10477  10200   9530   9591   8780
  9708   9749  10424   8796   8952  10358   9195   9367   9154   8618
 10136   8884  10197   9472   9141   9774   9026   9514   9945   8416
  9008  10001   9679   8854   9217   8927   9310  10315   9682   9176
  9301   9186   8765   9407   9358   8626   9943   8780   9384   9837
```
"""
function data(ctx::Context)
    out = Array{pxtype(ctx)}(undef, size(ctx)...)
    @eval ccall((:ScanImageTiffReader_GetData, $(_libname)), Int,
                (Ptr{Context}, Ptr{Cvoid}, Csize_t), $(Ref(ctx)), $(out), sizeof($(out)));

    if ctx.log != C_NULL
        throw(ErrorException(unsafe_string(ctx.log)))
    end

    out
end

"""
    length(ctx::Context)

Return the number of planes in the image stack.

# Examples
```jldoctest
ScanImageTiffReader.open(mytif) do io
    length(io)
end
# output
10
```
"""
function length(ctx::Context)
    n = @eval ccall((:ScanImageTiffReader_GetImageDescriptionCount, $(_libname)),
                    Int32, (Ptr{Context}, ), $(Ref(ctx)))

    if ctx.log != C_NULL
        throw(ErrorException(unsafe_string(ctx.log)))
    end

    n
end

"""
    description(ctx::Context, iframe::Int)

Return the contents of the image description tag for frame `iframe`.

# Examples
```jldoctest
desc = ScanImageTiffReader.open(mytif) do io
    description(io, 1)
end
print(desc)
# output
frameNumbers = 1
acquisitionNumbers = 1
frameNumberAcquisition = 1
frameTimestamps_sec = 0.000000
acqTriggerTimestamps_sec =
nextFileMarkerTimestamps_sec =
endOfAcquisition =  0
endOfAcquisitionMode = 0
dcOverVoltage = 0
epoch = [2016 6 4 13 51 7.8046]
```
"""
function description(ctx::Context, iframe::Int)
    sz = @eval ccall((:ScanImageTiffReader_GetImageDescriptionSizeBytes, $(_libname)),
                     Csize_t, (Ptr{Context}, Cint), $(Ref(ctx)), $(iframe) - 1) # convert to zero based
    sz != 0 || return ""

    str = Vector{UInt8}(undef, sz)
    @eval ccall((:ScanImageTiffReader_GetImageDescription,$(_libname)), Csize_t,
                (Ptr{Context}, Cint, Ptr{UInt8}, Csize_t),
                $(Ref(ctx)), $(iframe) - 1, $(str), sizeof($(str))) # convert to zero based
    if ctx.log != C_NULL
        throw(ErrorException(unsafe_string(ctx.log)))
    end

    unsafe_string(pointer(str))
end

"""
    metadata(ctx::Context)

Read the ScanImage metadata section from the file.

This data section is not part of the Tiff specification, so common Tiff readers
will not be able to access this data.

In ScanImage 2016 and later, this may be a JSON string.  For previous versions of
ScanImage, this is a bytestring that must be deserialized in MATLAB.

# Examples

```jldoctest
ScanImageTiffReader.open(mytif) do io
    JSON.parse(metadata(io))
end
# output
Dict{String,Any} with 2 entries:
  "SI"        => Dict{String,Any}("hConfigurationSaver"=>Dict{String,Any}("usrF…
  "RoiGroups" => Dict{String,Any}("photostimRoiGroups"=>nothing,"imagingRoiGrou…
```
"""
function metadata(ctx::Context)
    sz = @eval ccall((:ScanImageTiffReader_GetMetadataSizeBytes, $(_libname)),
                     Csize_t, (Ptr{Context}, ), $(Ref(ctx)))
    sz != 0 || return ""

    str = Vector{UInt8}(undef, sz)
    @eval ccall((:ScanImageTiffReader_GetMetadata, $(_libname)), Csize_t,
                (Ptr{Context}, Ptr{UInt8}, Csize_t),
                $(Ref(ctx)), $(str), sizeof($(str)))
    if ctx.log != C_NULL
        throw(ErrorException(unsafe_string(ctx.log)))
    end

    unsafe_string(pointer(str))
end

end


text:
using MultiScaleArrays
using OrdinaryDiffEq, DiffEqBase, Test, StochasticDiffEq, Statistics

struct Cell2{B} <: AbstractMultiScaleArrayLeaf{B}
    values::Vector{B}
end
struct Population2{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArrayHead{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end

cell1 = Cell2([1.0; 2.0; 3.0])
cell2 = Cell2([4.0; 5])
cell3 = Cell2([3.0; 2.0; 5.0])
cell4 = Cell2([4.0; 6])
pop = construct(Population2, deepcopy([cell1, cell2, cell3, cell4]))

function cell_ode4(dcell,cell,p,t)
  m = mean(cell)
  for i in eachindex(cell)
    dcell[i] = -m*cell[i]
  end
end

function f4(dembryo,embryo,p,t)
  for (cell, dcell) in LevelIter(1,embryo, dembryo)
    cell_ode4(dcell,cell,p,t)
  end
end

tstop = [0.5]

condition = function (u, t, integrator)
  t ∈ tstop
end

affect! = function (integrator)
  add_node!(integrator, integrator.u.nodes[1])
end

growing_cb = DiscreteCallback(condition, affect!)

prob = ODEProblem(f4, deepcopy(pop), (0.0, 1.0))

add_node!(pop,pop.nodes[1])

sol = solve(prob, Tsit5(), callback=growing_cb, tstops=tstop)

sol = solve(prob, Rosenbrock23(), callback=growing_cb,tstops=tstop)

@test length(sol[end]) == 13

affect_del! = function (integrator)
    remove_node!(integrator, 1)
end

shrinking_cb = DiscreteCallback(condition, affect_del!)

prob = ODEProblem(f4, deepcopy(pop), (0.0, 1.0))
sol = solve(prob, Tsit5(), callback=shrinking_cb, tstops=tstop)

prob = ODEProblem(f4, deepcopy(pop), (0.0, 1.0))
sol = solve(prob, Rosenbrock23(),callback=shrinking_cb,tstops=tstop)
@test length(sol[end]) == 10

println("Do the SDE Part")

function g4(du,u,p,t)
    for i in eachindex(u)
        du[i] = 0.1u[i]
    end
end

prob = SDEProblem(f4, g4, deepcopy(pop), (0.0, 1.0))

sol = solve(prob, SOSRI(), callback=growing_cb, tstops=tstop)

@test length(sol[end]) == 16

sol = solve(prob, SOSRI(), callback=shrinking_cb, tstops=tstop)

@test length(sol[end]) == 10


text:
@testset "CoreModelCoupled generic interface" begin
    model = load_model(CoreModelCoupled, model_paths["e_coli_core.mat"])

    @test reaction_stoichiometry(model, "EX_ac_e") == Dict("ac_e" => -1)
    @test reaction_stoichiometry(model, 44) == Dict("ac_e" => -1)

end

@testset "Conversion from and to StandardModel" begin
    cm = load_model(CoreModelCoupled, model_paths["e_coli_core.mat"])

    sm = convert(StandardModel, cm)
    cm2 = convert(CoreModelCoupled, sm)

    @test Set(reactions(cm)) == Set(reactions(sm))
    @test Set(reactions(cm)) == Set(reactions(cm2))

    @test reaction_gene_association(sm, reactions(sm)[1]) ==
          reaction_gene_association(cm, reactions(sm)[1])

    @test reaction_gene_association_vec(cm)[1:3] == [
        [["b3916"], ["b1723"]],
        [
            ["b0902", "b0903", "b2579"],
            ["b0902", "b0903"],
            ["b0902", "b3114"],
            ["b3951", "b3952"],
        ],
        [["b4025"]],
    ]
end


text:
rotatereferenceatlas = x -> reverse(permutedims(x, (1, 3, 2)), dims=(3,))
# Plot the reference volume and probe locations
function plotreferencevolume(S; dotext=true, dostructures = true, ids=:targets, resolution=(1920, 1080), kwargs...)
    channels = AN.getchannels(S)
    vol, info = AN.gettemplatevolume()
    vol = Array{Float16}(vol)
    grad = cgrad(cgrad([RGBA(1.0, 1.0, 1.0, 0.0), RGBA(0.0, 0.0, 0.0, 1.0)]), 100)
    #coords = [(1:s).*25 for s ∈ size(vol)] # 25 μm resolution
    #vol = reverse(permutedims(vol, [1 3 2]), dims=(1,3))
    #coords = [1:x for x ∈ size(vol)].*50
    s = Scene(backgroundcolor = RGBA(1.0, 1.0, 1.0, 0.0), resolution = resolution, transparency=true)
    coords = [1:s for s ∈ size(rotatereferenceatlas(vol))./2]
    coords[3] = .-coords[3]

    markercolor = :gray
    if ids isa Symbol
        if ids == :all
            # Get and plot a bunch of structures
            #ids = AN.getallstructureids()
            anns = unique(["FRP", "MO", "SS", "GU", "VISC", "AUD", "ACA", "PL", "ILA", "ORB", "AI", "RSP", "PTLp", "TEa", "PERI", "ECT", "OLF", "VISp", "VISl", "VISrl", "VISam", "VISpm", "VIS", "VISal","VISmma","VISmmp","VISli", "LGd","LD", "LP", "VPM", "TH", "MGm","MGv","MGd","PO","LGv","VL",
            "VPL","POL","Eth","PoT","PP","PIL","IntG","IGL","SGN","VPL","PF","RT", "CA1", "CA2","CA3", "DG", "SUB", "POST","PRE","ProS","HPF", "MB","SCig","SCiw","SCsg","SCzo","PPT","APN","NOT","MRN","OP","LT","RPF","CP"])
        elseif ids == :corticaltargets
            anns = unique(["VISp", "VISl", "VISrl", "VISal", "VISpm", "VISam"])
        elseif ids == :targets # As listed in the visual coding white paper
            anns =  unique(["VISp", "VISl", "VISrl", "VISal", "VISpm", "VISam", "CA1", "CA3", "DG", "SUB", "ProS", "LGd", "LP", "APN"])
        else
            @error "`ids` keyword argument not valid "
        end
        t = AN.getstructuretree()
        D = t.get_id_acronym_map()
        ids = getindex.((D,), anns)
    end

    f = volume!(s, coords..., rotatereferenceatlas(vol)[1:2:end, 1:2:end, 1:2:end]; algorithm=:mip, colorrange=extrema(vol), colormap=collect(grad), ticks=nothing, fxaa=true, transparency=true, kwargs...)
    #s.plots[1].attributes[:fxaa] = true

    chrome = FileIO.load(download("https://raw.githubusercontent.com/nidorx/matcaps/master/1024/E6BF3C_5A4719_977726_FCFC82.png"))
    for probeid in unique(channels.probe_id)
        (x, y, z) = AN.getprobecoordinates(S, probeid)./50
        if !any(length.((x, y, z)).==0)
            f = meshscatter!(s, x, z, -y; markersize=1.0, fxaa=true, color=markercolor, matcap=chrome, shading=true, transparency=true, kwargs...)

            if dotext
                _, idx = findmax(z) # The highest unit
                text!(s, string(probeid), position=Vec3f0(x[idx], z[idx], -y[idx]), space=:data, textsize=5, align=(:right, :bottom), rotation=Quaternion((-0.3390201, 0.33899, 0.6205722, -0.620517)), transparency=true)
            end
        end
    end
    #scale!(s, 1, 1, -1)
    if dostructures
        channels = subset(AN.getchannels(),    :ecephys_probe_id           => ByRow(!ismissing),
                                                :ecephys_structure_id       => ByRow(!ismissing),
                                                :ecephys_structure_acronym  => ByRow(x->x!=="grey"))
        if isnothing(ids)
            ids = unique(channels[!, :ecephys_structure_id])
            #ids = ids[AN.getstructuretreedepth.(ids) .< 9]
        end
        for id ∈ ids
            try
                try
                    mask, _ = AN.getstructuremask(id)
                catch
                    mask = AN.buildstructuremask(id)
                end
                mask = Array{Float64}(rotatereferenceatlas(mask)[1:2:end, 1:2:end, 1:2:end])
                mc_algo = NaiveSurfaceNets(iso=0, insidepositive=true)
                m = GeometryBasics.Mesh(mask, mc_algo; origin=[min(coords[i]...) for i ∈ 1:length(coords)], widths=[abs(-(extrema(coords[i])...)) for i ∈ 1:length(coords)])
                c = AN.getstructurecolor(id)
                f = mesh!(s, m; color=RGBA(c.r, c.g, c.b, 0.41), fxaa=true, shading=true, transparency=true, kwargs...)
            catch y
                @warn y
            end
        end
    end
    return s
end

function exportreferencevolume(S, file::String="plot.html"; ids=:targets, kwargs...)
    open(file, "w") do io
            println(io, """
            <html>
                <head>
                </head>
                <body>
            """)
            show(io, MIME"text/html"(), Page(exportable=true, offline=true))
            fig = plotreferencevolume(S; dostructures=true, ids, kwargs...)
            rotate_cam!(fig, Vec3f0(0, 2.35, 0))
            show(io, MIME"text/html"(), fig)
            println(io, """
                </body>
            </html>
            """)
    end
end

function formattedreferencevolume(S, file::String="plot.html")
    exportreferencevolume(S, file;
                                    dostructures=true,
                                    ids=:targets,
                                    show_axis=false,
                                    shading=true, resolution=(1080, 1080));
end


text:
function build(contractType::Type{T}, options::Dict{String,Any})::AbstractAssetModel where {T<:AbstractAssetModel}

    # initialize -
    model = eval(Meta.parse("$(contractType)()")) # empty contract model -

    # if we have options, add them to the contract model -
    if (isempty(options) == false)

        # for the result of the fields, let's lookup in the dictionary.
        # error state: if the dictionary is missing a value -
        for field_name_symbol ∈ fieldnames(contractType)

            # convert the field_name_symbol to a string -
            field_name_string = string(field_name_symbol)

            # check the for the key -
            if (haskey(options, field_name_string) == false)
                throw(ArgumentError("dictionary is missing: $(field_name_string)"))
            end

            # get the value -
            value = options[field_name_string]

            # set -
            setproperty!(model, field_name_symbol, value)
        end
    end

    # return -
    return model
end

function build(modelType::Type{CRRContractPremiumLatticeModel}, options::Dict{String,Any})::CRRContractPremiumLatticeModel

     # initialize -
     model = eval(Meta.parse("$(modelType)()")) # empty contract model -

     # if we have options, add them to the contract model -
     if (isempty(options) == false)

         # for the result of the fields, let's lookup in the dictionary.
         # error state: if the dictionary is missing a value -
         for field_name_symbol ∈ fieldnames(modelType)

             # convert the field_name_symbol to a string -
             field_name_string = string(field_name_symbol)

             # check the for the key -
             if (haskey(options, field_name_string) == false)
                 throw(ArgumentError("dictionary is missing: $(field_name_string)"))
             end

             # get the value -
             value = options[field_name_string]

             # set -
             setproperty!(model, field_name_symbol, value)
         end
     end

     # return -
     return model
end

function build(type::Type{PQContractPremiumLatticePoint};
    i::Int64=1, j::Int64=1, s::Int64=1, d::Int64 = 1)::PQContractPremiumLatticePoint

    point_model = PQContractPremiumLatticePoint()
    point_model.i = i
    point_model.j = j
    point_model.s = s
    point_model.d = d

    # return -
    return point_model
end

function build(type::Type{CRRLatticeModel}; number_of_levels::Int64 = 2, T::Float64 = (1 / 365), σ::Float64 = 0.15,
    Sₒ::Float64 = 1.0, μ::Float64 = 0.0015)::Union{ArgumentError,CRRLatticeModel}

    # check: number_of_levels ≥ 2
    if (number_of_levels < 2)
        return ArgumentError("Check failure: number_of_levels must be ≥ 2")
    end

    # initialize -
    lattice_model = CRRLatticeModel()

    # compute the movement up, down and probability -
    ΔT = T / number_of_levels
    u = exp(σ * √ΔT)
    d = exp(-σ * √ΔT)
    p = (exp(µ * ΔT) - d) / (u - d)

    # compute connectivity -
    number_items_per_level = [i for i = 1:number_of_levels]
    tmp_array = Array{Int64,1}()
    theta = 0
    for value in number_items_per_level
        for _ = 1:value
            push!(tmp_array, theta)
        end
        theta = theta + 1
    end

    N = sum(number_items_per_level[1:(number_of_levels-1)])
    connectivity_index_array = Array{Int64,2}(undef, N, 3)
    for row_index = 1:N

        # index_array[row_index,1] = tmp_array[row_index]
        connectivity_index_array[row_index, 1] = row_index
        connectivity_index_array[row_index, 2] = row_index + 1 + tmp_array[row_index]
        connectivity_index_array[row_index, 3] = row_index + 2 + tmp_array[row_index]
    end
    lattice_model.connectivity = connectivity_index_array

    # init the tree -
    total_number_of_lattice_nodes = connectivity_index_array[end, end]
    number_of_nodes_to_evaluate = connectivity_index_array[end, 1]
    tree_value_array = Array{Float64,2}(undef, total_number_of_lattice_nodes, 3) # nodes x 3 = col1: underlying price, col2: intrinsic value, col3: option price

    # First: let's compute the underlying price on the lattice -
    tree_value_array[1, 1] = Sₒ
    for node_index ∈ 1:number_of_nodes_to_evaluate

        # get index -
        parent_node_index = connectivity_index_array[node_index, 1]
        up_node_index = connectivity_index_array[node_index, 2]
        down_node_index = connectivity_index_array[node_index, 3]

        # compute prices -
        parent_price = tree_value_array[parent_node_index, 1]
        up_price = parent_price * u
        down_price = parent_price * d

        # store prices -
        tree_value_array[up_node_index, 1] = up_price
        tree_value_array[down_node_index, 1] = down_price
    end

    # Second: let's add the data to the lattice model -
    lattice_model.data = tree_value_array

    # add other stuff -
    lattice_model.p = p
    lattice_model.μ = μ
    lattice_model.u = u
    lattice_model.d = d
    lattice_model.ΔT = ΔT
    lattice_model.σ = σ
    lattice_model.T = T
    lattice_model.Sₒ = Sₒ

    # return the model -
    return lattice_model
end

function build(latticeType::Type{CRRJITContractPremiumLatticeModel},contractType::Type{Y}, underlying::Array{Float64,1},
    strike::Array{Float64,1}, iv::Array{Float64,1}, dte::Array{Float64,1}; number_of_levels::Int64 = 100, μ::Float64 = 0.0045)::CRRJITContractPremiumLatticeModel  where {Y<:AbstractDerivativeContractModel}

    # show -
    @info "Building a $(latticeType) lattice with a $(contractType) contract"

    # data fields -
    # underlying::Array{Float64,1}
    # strike::Union{Array{Float64,1}, Array{Float64,2}}
    # iv::Array{Float64,1}
    # dte::Array{Float64,1}
    # number_of_levels::Int64
    # risk_free_rate::Float64
    # contractType::Type{T} where {T <: AbstractDerivativeContractModel}

    # build an empty CRRContractPremiumLatticeModel -
    model = CRRJITContractPremiumLatticeModel()
    model.underlying = underlying
    model.strike = strike
    model.iv = iv
    model.dte = dte
    model.number_of_levels = number_of_levels
    model.risk_free_rate = μ
    model.contractType = contractType

    # return -
    return model
end

function build(latticeType::Type{CRRContractPremiumLatticeModel},contractType::Type{Y}, underlying::Array{Float64,1},
    strike::Array{Float64,1}, iv::Array{Float64,1}, dte::Array{Float64,1}; number_of_levels::Int64 = 100, μ::Float64 = 0.0045)::CRRContractPremiumLatticeModel  where {Y<:AbstractDerivativeContractModel}

    # show -
    @info "Building $(latticeType). This may take several minutes."

    # initialize -
    simulation_archive = Dict{NamedTuple, Float64}()

    # build an empty contract -
    contract = build(contractType, Dict{String,Any}())
    contract.number_of_contracts = 1
    contract.direction = 1

    # main loop to compute look ahead archive -
    for (d, T) ∈ enumerate(dte)
        for (s, Sₒ) ∈ enumerate(underlying)
            for (i, K) ∈ enumerate(strike)
                for (j, σ) ∈ enumerate(iv)

                    # build a binary tree with these market conditions -
                    model = build(CRRLatticeModel; Sₒ = Sₒ, number_of_levels = number_of_levels, σ = σ, T = (T / 365), μ = μ)

                    # set the parameters on the contract / for this version of the method we have a single contract -
                    contract.strike_price = K

                    # compute the premimum -
                    p = premium(contract, model)

                    # build the key -
                    key_tuple = (
                        s = s,
                        d = d,
                        i = i,
                        j = j
                    );

                    # cache -
                    simulation_archive[key_tuple] = p

                    @info (s,d,i,j)
                end
            end
        end
    end

    # build model -
    lattice_model = CRRContractPremiumLatticeModel()
    lattice_model.grid = simulation_archive
    lattice_model.underlying = underlying
    lattice_model.strike = strike
    lattice_model.iv = iv
    lattice_model.dte = dte
    lattice_model.number_of_levels = number_of_levels
    lattice_model.risk_free_rate = μ
    lattice_model.contractType = contractType

    # return -
    return lattice_model
end


text:
# Quadratic function

# A generic quadratic equation of the form:
# f = (x) -> ax^2 + bx + c
# f1 = (x) -> 2ax + b

include("derivative.jl")

function quadratic(f)
    f1 = derivative(f)

    c = f(0.0)
    b = f1(0.0)
    a = f(1.0) - b - c

    return (-b + sqrt(b^2 -4a*c + 0im)) / 2a, (-b - sqrt(b^2 -4a*c + 0im)) / 2a

end

max_stars_repo_name: neilpanchal/Julia-Algorithms

max_issues_repo_name: neilpanchal/Julia-Algorithms

max_forks_repo_name: neilpanchal/Julia-Algorithms


text:
"""
```
homotopy!(m, xₙ₋₁; step = .1, pnorm = Inf, verbose = :none, kwargs...)
```

solves the system of equations characterizing a risk-adjusted linearization by a homotopy method with
embedding parameter ``q``, which steps from 0 to 1, with ``q = 1`` obtaining the true solution.

Currently, the only algorithm for choosing ``q`` is a simple uniform step search. Given a step size
``\\Delta``, we solve the homotopy starting from ``q = \\Delta`` and increase ``q`` by ``\\Delta``
until ``q`` reaches 1 or passes 1 (in which case, we force ``q = 1``).

### Types:
- `S1 <: Number`

### Inputs
- `m::RiskAdjustedLinearization`: object holding functions needed to calculate
    the risk-adjusted linearization
- `xₙ₋₁::AbstractVector{S1}`: initial guess for ``(z, y, \\Psi)``

### Keywords
- `step::Float64`: size of the uniform step from `step` to 1.
- `pnorm::Float64`: norm under which to evaluate the errors after homotopy succeeds.
- `sparse_jacobian::Bool = false`: if true, exploit sparsity in the Jacobian in calls to `nlsolve` using SparseDiffTools.jl.
    If `jac_cache` and `sparsity` are `nothing`, then `homotopy!` will attempt to determine the sparsity pattern.
- `sparsity::Union{AbstractArray, Nothing} = nothing`: sparsity pattern for the Jacobian in calls to `nlsolve`
- `colorvec = nothing`: matrix coloring vector for sparse Jacobian in calls to `nlsolve`
- `jac_cache = nothing`: pre-allocated Jacobian cache for calls to `nlsolve` during the numerical algorithms
- `sparsity_detection::Bool = false`: If true, use SparsityDetection.jl to detect sparsity pattern (only relevant if
    both `jac_cache` and `sparsity` are `nothing`). If false,  then the sparsity pattern is
    determined by using finite differences to calculate a Jacobian and assuming any zeros will always be zero.
    Currently, SparsityDetection.jl fails to work.
- `verbose::Symbol`: verbosity of information printed out during solution.
    a) `:low` -> statement when homotopy continuation succeeds
    b) `:high` -> statement when homotopy continuation succeeds and for each successful iteration
"""
function homotopy!(m::RiskAdjustedLinearization, xₙ₋₁::AbstractVector{S1};
                   step::Float64 = .1, pnorm::Float64 = Inf,
                   autodiff::Symbol = :central,
                   sparse_jacobian::Bool = false,
                   sparsity::Union{AbstractArray, Nothing} = nothing, colorvec = nothing,
                   jac_cache = nothing, sparsity_detection::Bool = true,
                   verbose::Symbol = :none,
                   kwargs...) where {S1 <: Number}
    # Set up
    nl = nonlinear_system(m)
    li = linearized_system(m)
    _my_eqn = if Λ_eltype(nl) <: RALF1 && Σ_eltype(nl) <: RALF1 # only difference in this block and the next block
        (F, x, q) -> _homotopy_equations1(F, x, m, q)                # is the number of args to retrieve 𝒱_sss and JV
    else
        (F, x, q) -> _homotopy_equations2(F, x, m, q)
    end

    qguesses = step:step:1.
    if qguesses[end] != 1.
        qguesses = vcat(qguesses, 1.)
    end
    for (i, q) in enumerate(qguesses)
        solve_steadystate!(m, getvecvalues(m), _my_eqn, q; autodiff = autodiff,
                           sparse_jacobian = sparse_jacobian,
                           sparsity = sparsity, colorvec = colorvec,
                           jac_cache = jac_cache,
                           sparsity_detection = sparsity_detection,
                           verbose = verbose, kwargs...)

        if verbose == :high
            println("Success at iteration $(i) of $(length(qguesses))")
        end
    end

    if verbose in [:low, :high]
        errvec = steady_state_errors(m)

        println("Homotopy succeeded!")
        println("Error under norm = $(pnorm) is $(norm(errvec, pnorm)).")
    end

    update!(m)

    return m
end

function solve_steadystate!(m::RiskAdjustedLinearization, x0::AbstractVector{S1}, f::Function, q::Float64;
                            sparse_jacobian::Bool = false,
                            sparsity::Union{AbstractArray, Nothing} = nothing, colorvec = nothing,
                            jac_cache = nothing, sparsity_detection::Bool = true, autodiff::Symbol = :central,
                            verbose::Symbol = :none, kwargs...) where {S1 <: Real}

    if sparse_jacobian # Exploit sparsity?
        nlsolve_jacobian!, jac =
            construct_sparse_jacobian_function(m, (F, x) -> f(F, x, q), :homotopy, autodiff;
                                               sparsity = sparsity, colorvec = colorvec,
                                               jac_cache = jac_cache,
                                               sparsity_detection = sparsity_detection)
        out = nlsolve(OnceDifferentiable((F, x) -> f(F, x, q), nlsolve_jacobian!, x0, copy(x0), jac), x0; kwargs...)
    else
        # Need to declare chunk size to ensure no problems with reinterpreting the cache
        # Potential solution: https://discourse.julialang.org/t/issue-with-pdmp-and-forwardiff-differentialequation/17925/22
        # Essentially, it may require developing another cache.
        out = nlsolve(OnceDifferentiable((F, x) -> f(F, x, q), x0, copy(x0), autodiff,
                                         ForwardDiff.Chunk(ForwardDiff.pickchunksize(min(m.Nz, m.Ny)))), x0; kwargs...)
    end

    if out.f_converged
        N_zy = m.Nz + m.Ny
        m.z .= out.zero[1:m.Nz]
        m.y .= out.zero[(m.Nz + 1):N_zy]
        m.Ψ .= reshape(out.zero[(N_zy + 1):end], m.Ny, m.Nz)
    else
        if verbose == :high
            println(out)
        end
        throw(RALHomotopyError("A solution for (z, y, Ψ) to the state transition, expectational, " *
                               "and linearization equations could not be found when the embedding " *
                               "parameter q equals $(q)"))
    end
end

function _homotopy_equations1(F::AbstractArray, x::AbstractArray, m::RiskAdjustedLinearization, q::Number)
    # Unpack
    z = @view x[1:m.Nz]
    y = @view x[(m.Nz + 1):(m.Nz + m.Ny)]
    Ψ = @view x[(m.Nz + m.Ny + 1):end]
    Ψ = reshape(Ψ, m.Ny, m.Nz)

    # Given coefficients, update the model
    update!(m.nonlinear, z, y, Ψ)
    update!(m.linearization, z, y, Ψ)

    # Calculate residuals
    μ_sss              = get_tmp(m.nonlinear.μ.cache, z, y, (1, 1)) # select the first DiffCache b/c that one
    ξ_sss              = get_tmp(m.nonlinear.ξ.cache, z, y, (1, 1)) # corresponds to autodiffing both z and y
    𝒱_sss              = get_tmp(m.nonlinear.𝒱.cache, z, Ψ, (1, 1)) # This line is different than in _homotopy_equations2
    Γ₁                 = get_tmp(m.linearization.μz.cache, z, y, (1, 1))
    Γ₂                 = get_tmp(m.linearization.μy.cache, z, y, (1, 1))
    Γ₃                 = get_tmp(m.linearization.ξz.cache, z, y, (1, 1))
    Γ₄                 = get_tmp(m.linearization.ξy.cache, z, y, (1, 1))
    JV                 = get_tmp(m.linearization.J𝒱.cache, z, Ψ, (1, 1)) # This line is different than in _homotopy_equations2
    F[1:m.Nz]          = μ_sss - z
    F[(m.Nz + 1):(m.Nz + m.Ny)] = ξ_sss + m.linearization[:Γ₅] * z + m.linearization[:Γ₆] * y + q * 𝒱_sss
    F[((m.Nz + m.Ny) + 1):end]  = Γ₃ + Γ₄ * Ψ + (m.linearization[:Γ₅] + m.linearization[:Γ₆] * Ψ) * (Γ₁ + Γ₂ * Ψ) + q * JV
end

function _homotopy_equations2(F::AbstractArray, x::AbstractArray, m::RiskAdjustedLinearization, q::Number)

    # Unpack
    z = @view x[1:m.Nz]
    y = @view x[(m.Nz + 1):(m.Nz + m.Ny)]
    Ψ = @view x[(m.Nz + m.Ny + 1):end]
    Ψ = reshape(Ψ, m.Ny, m.Nz)

    # Given coefficients, update the model
    update!(m.nonlinear, z, y, Ψ)
    update!(m.linearization, z, y, Ψ)

    # Calculate residuals
    μ_sss              = get_tmp(m.nonlinear.μ.cache, z, y, (1, 1)) # select the first DiffCache b/c that one
    ξ_sss              = get_tmp(m.nonlinear.ξ.cache, z, y, (1, 1)) # corresponds to autodiffing both z and y
    𝒱_sss              = get_tmp(m.nonlinear.𝒱.cache, z, y, Ψ, z, (1, 1))
    Γ₁                 = get_tmp(m.linearization.μz.cache, z, y, (1, 1))
    Γ₂                 = get_tmp(m.linearization.μy.cache, z, y, (1, 1))
    Γ₃                 = get_tmp(m.linearization.ξz.cache, z, y, (1, 1))
    Γ₄                 = get_tmp(m.linearization.ξy.cache, z, y, (1, 1))
    JV                 = get_tmp(m.linearization.J𝒱.cache, z, y, Ψ, (1, 1))
    F[1:m.Nz]          = μ_sss - z
    F[(m.Nz + 1):(m.Nz + m.Ny)] = ξ_sss + m.linearization[:Γ₅] * z + m.linearization[:Γ₆] * y + q * 𝒱_sss
    F[((m.Nz + m.Ny) + 1):end]  = Γ₃ + Γ₄ * Ψ + (m.linearization[:Γ₅] + m.linearization[:Γ₆] * Ψ) * (Γ₁ + Γ₂ * Ψ) + q * JV
end

mutable struct RALHomotopyError <: Exception
    msg::String
end
Base.showerror(io::IO, ex::RALHomotopyError) = print(io, ex.msg)


text:
using Test
import VTUFileHandler

@test include("test_VTU.jl")

max_stars_repo_name: nagelt/VTUFileHandler

max_issues_repo_name: nagelt/VTUFileHandler

max_forks_repo_name: nagelt/VTUFileHandler


text:
mutable struct Node
    index::Int64
    centroid::Vector{Float64}
    ifaces::Vector{Int64}
    icells::Vector{Int64}
    print::Function
    function Node()
        this = new()
        this.print = function() printnode(this) end
        this
    end
end

mutable struct Face
    index::Int64
    centroid::Vector{Float64}
    inodes::Vector{Int64}
    iowner::Int64
    ineighbour::Int64
    area::Float64
    sf::Vector{Float64}
    cn::Vector{Float64}
    geodiff::Float64
    T::Vector{Float64}
    gf::Float64
    walldist::Float64
    iownerneighbourcoeff::Int64
    ineighbourownercoeff::Int64
    print::Function
    function Face()
        this = new()
        this.print = function() printface(this) end
        this
    end
end

mutable struct Boundary
    index::Int64
    name::String
    type::String
    nfaces::Int64
    startface::Int64
    print::Function
    function Boundary()
        this = new()
        this.print = function() printboundary(this) end
        this
    end
end

mutable struct Cell
    index::Int64
    centroid::Vector{Float64}
    inodes::Vector{Int64}
    ifaces::Vector{Int64}
    ineighbours::Vector{Int64}
    volume::Float64
    facesign::Vector{Int64}
    nneighbours::Int64
    print::Function
    function Cell()
        this = new()
        this.print = function() printcell(this) end
        this
    end
end

mutable struct Mesh
    nodes::Vector{Node}
    faces::Vector{Face}
    boundaries::Vector{Boundary}
    cells::Vector{Cell}
    print::Function
    function Mesh(nodes, faces, boundaries, cells)
        this = new()
        this.nodes = nodes
        this.faces = faces
        this.boundaries = boundaries
        this.cells = cells
        this.print = function() printmesh(this) end
        this
    end
end

function printnode(node::Node)
    print("\n")
    println("        index: ", node.index)
    println("     centroid: ", node.centroid)
    println("       iFaces: ", node.ifaces)
    println("       iCells: ", node.icells)
    print("\n")
end

function printface(face::Face)
    print("\n")
    println("                   index: ", face.index)
    println("                centroid: ", face.centroid)
    println("                  iNodes: ", face.inodes)
    println("                  iOwner: ", face.iowner)
    println("              iNeighbour: ", face.ineighbour)
    println("                    area: ", face.area)
    println("                      Sf: ", face.sf)
    println("                      CN: ", face.cn)
    println("                 geodiff: ", face.geodiff)
    println("                       T: ", face.T)
    println("                      gf: ", face.gf)
    println("                walldist: ", face.walldist)
    println("    iOwnerNeighbourCoeff: ", face.iownerneighbourcoeff)
    println("    iNeighbourOwnerCoeff: ", face.ineighbourownercoeff)
    print("\n")
end

function printboundary(boundary::Boundary)
    print("\n")
    println("             index: ", boundary.index)
    println("              name: ", boundary.name)
    println("              type: ", boundary.type)
    println("    numberOfBFaces: ", boundary.nfaces)
    println("         startFace: ", boundary.startface)
    print("\n")
end

function printcell(cell::Cell)
    print("\n")
    println("                 index: ", cell.index)
    println("              centroid: ", cell.centroid)
    println("                iNodes: ", cell.inodes)
    println("                iFaces: ", cell.ifaces)
    println("           iNeighbours: ", cell.ineighbours)
    println("                volume: ", cell.volume)
    println("              faceSign: ", cell.facesign)
    println("    numberOfNeighbours: ", cell.nneighbours)
    print("\n")
end

function printmesh(mesh::Mesh)
    print("\n")
    println("                    nodes: ", typeof(mesh.nodes))
    println("            numberOfNodes: ", length(mesh.nodes))
    println("                    faces: ", typeof(mesh.faces))
    println("            numberOfFaces: ", length(mesh.faces))
    println("                 elements: ", typeof(mesh.cells))
    println("         numberOfElements: ", length(mesh.cells))
    println("               boundaries: ", typeof(mesh.boundaries))
    println("       numberOfBoundaries: ", length(mesh.boundaries))
    println("    numberOfInteriorFaces: ")
    print("\n")
end


text:
module TestConstraintScalarize

using Test

using VecMathOptInterface
const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_scalarize()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.Scalarize{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = [
            "test_basic_$(F)_$(S)" for F in [
                "VectorOfVariables",
                "VectorAffineFunction",
                "VectorQuadraticFunction",
            ] for S in ["Nonnegatives", "Nonpositives", "Zeros"]
        ],
    )
    MOI.empty!(bridged_mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0, 0.0, 2.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [-3, -1],
        )
    MOI.Test.test_conic_linear_VectorOfVariables(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.Zeros,
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        3,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, 0),
            (MOI.VariableIndex, MOI.GreaterThan{Float64}, 0),
        ),
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorOfVariables,
                MOI.Nonnegatives,
            }(),
        ),
    )
    func = MOI.get(bridged_mock, MOI.ConstraintFunction(), ci)
    MOI.delete(bridged_mock, func.variables[2])
    new_func = MOI.VectorOfVariables(func.variables[[1, 3]])
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) == new_func
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.Nonnegatives(2)
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        values = [1.0, 2.0]
        MOI.set(bridged_mock, attr, ci, values)
        @test MOI.get(bridged_mock, attr, ci) == values
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        2,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, 0),
            (MOI.VariableIndex, MOI.GreaterThan{Float64}, 0),
        ),
    )
    MOI.empty!(bridged_mock)
    # VectorAffineFunction-in-Nonnegatives
    # VectorAffineFunction-in-Zeros
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0, 0.0, 2.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [0, 2, 0],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [-3, -1],
        )
    MOI.Test.test_conic_linear_VectorAffineFunction(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.Zeros,
            }(),
        ),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        values = [1.0, -2.0]
        MOI.set(bridged_mock, attr, ci, values)
        @test MOI.get(bridged_mock, attr, ci) == values
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        3,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
        ),
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.Nonnegatives,
            }(),
        ),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        values = [1.0, -2.0, 3.0]
        MOI.set(bridged_mock, attr, ci, values)
        @test MOI.get(bridged_mock, attr, ci) == values
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        3,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
        ),
    )
    MOI.empty!(bridged_mock)
    # Test setting VectorAffineFunction with nonzero constants
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [0.0, 0.0]),
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [1.0, 0.75]),
    )
    MOI.Test.test_modification_func_vectoraffine_nonneg(bridged_mock, config)
    MOI.empty!(bridged_mock)
    # VectorOfVariables-in-Nonnegatives
    # VectorOfVariables-in-Nonpositives
    # VectorOfVariables-in-Zeros
    # VectorAffineFunction-in-Zeros
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [-4, -3, 16, 0],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [7, 2, -4],
        )
    MOI.Test.test_conic_linear_VectorOfVariables_2(bridged_mock, config)
    MOI.empty!(bridged_mock)
    # VectorAffineFunction-in-Nonnegatives
    # VectorAffineFunction-in-Nonpositives
    # VectorAffineFunction-in-Zeros
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [-4, -3, 16, 0],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [0],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [0],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [7, 2, -4, 7],
        )
    MOI.Test.test_conic_linear_VectorAffineFunction_2(bridged_mock, config)
    MOI.Test.test_constraint_ConstraintPrimalStart(
        bridged_mock,
        MOI.Test.Config(),
    )
    MOI.Test.test_constraint_ConstraintDualStart(
        bridged_mock,
        MOI.Test.Config(),
    )
    return
end

end  # module

TestConstraintScalarize.runtests()


text:
###
### This file is generated by running
### ` julia generate_buildjl.jl L/Libtask/build_tarballs.jl`
### in the Yggdrasil root directory, with 2 updates:
### 1. add prefix tp products, see https://github.com/JuliaPackaging/Yggdrasil#binaryproviderjl,
### 2. products filter
###

using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libtask_v1_0"], :libtask_v1_0),
    LibraryProduct(prefix, ["libtask_v1_1"], :libtask_v1_1),
    LibraryProduct(prefix, ["libtask_v1_2"], :libtask_v1_2),
    LibraryProduct(prefix, ["libtask_v1_3"], :libtask_v1_3),
]

products_tmp = filter(products) do prod
    endswith(prod.libnames[1], "$(VERSION.major)_$(VERSION.minor)")
end
length(products_tmp) == 0 && (products_tmp = [products[end]])
products = products_tmp

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaBinaryWrappers/Libtask_jll.jl/releases/download/Libtask-v0.3.1+0"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc) => ("$bin_prefix/Libtask.v0.3.1.aarch64-linux-gnu.tar.gz", "ad342e90969280fcd198993947fb27a052357c278e40645574b93f13b9f49f25"),
    Linux(:i686, libc=:glibc) => ("$bin_prefix/Libtask.v0.3.1.i686-linux-gnu.tar.gz", "a4061a2098da96cd6c49300b2895a61ce02d822c603ba85995d300149c9cec5e"),
    Windows(:i686) => ("$bin_prefix/Libtask.v0.3.1.i686-w64-mingw32.tar.gz", "e38da986025dac5767da841cc3666973a303061f4b853c0821fc25de4802f674"),
    MacOS(:x86_64) => ("$bin_prefix/Libtask.v0.3.1.x86_64-apple-darwin14.tar.gz", "93bf6be03236252eea33998c5fa4f0a8173d44b20266e8e05473a2ae957eebde"),
    Linux(:x86_64, libc=:glibc) => ("$bin_prefix/Libtask.v0.3.1.x86_64-linux-gnu.tar.gz", "349fc5fdd1d1234f1c1467fd655b9a5e18ac926e5b16db4ba7743b75a9540739"),
    Windows(:x86_64) => ("$bin_prefix/Libtask.v0.3.1.x86_64-w64-mingw32.tar.gz", "837410807715caf15b931f04c968a3ca8967bd0d70df888d1ed535a1679030d0"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
function allexamples2d()
    include("arrangement2d.jl")
    include("biconnectedgraphs.jl")
    include("graphalgorithms.jl")
    include("orient2d.jl")
    include("point_in_polygon.jl")
    include("randomarrangement2d.jl")
    include("randombubbles.jl")
    include("randomlines.jl")
    include("randomshapes.jl")
    include("stresstest2d.jl")
    include("svg2lar.jl")
end

allexamples2d()


text:
using ClimateMachine
using ClimateMachine.ConfigTypes
using ClimateMachine.Mesh.Topologies: StackedCubedSphereTopology, grid1d
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Filters
using ClimateMachine.Atmos: AtmosFilterPerturbations
using ClimateMachine.DGMethods: ESDGModel, init_ode_state, courant
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.ODESolvers
using ClimateMachine.SystemSolvers
using ClimateMachine.VTK: writevtk, writepvtu
using ClimateMachine.GenericCallbacks:
    EveryXWallTimeSeconds, EveryXSimulationSteps
using ClimateMachine.Thermodynamics: soundspeed_air
using ClimateMachine.TemperatureProfiles
using ClimateMachine.VariableTemplates: flattenednames

using CLIMAParameters
using CLIMAParameters.Planet: R_d, cv_d, Omega, planet_radius, MSLP
import CLIMAParameters

using MPI, Logging, StaticArrays, LinearAlgebra, Printf, Dates, Test
using CUDA

const output_vtk = false

struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet();
#const X = 20
const X = 1
CLIMAParameters.Planet.planet_radius(::EarthParameterSet) = 6.371e6 / X
CLIMAParameters.Planet.Omega(::EarthParameterSet) = 7.2921159e-5 * X

# No this isn't great but w/e

include("DryAtmos.jl")

function sphr_to_cart_vec(vec, lat, lon)
    FT = eltype(vec)
    slat, clat = sin(lat), cos(lat)
    slon, clon = sin(lon), cos(lon)
    u = MVector{3, FT}(
        -slon * vec[1] - slat * clon * vec[2] + clat * clon * vec[3],
        clon * vec[1] - slat * slon * vec[2] + clat * slon * vec[3],
        clat * vec[2] + slat * vec[3],
    )
    return u
end

struct BaroclinicWave <: AbstractDryAtmosProblem end

function init_state_prognostic!(
    bl::DryAtmosModel,
    ::BaroclinicWave,
    state,
    aux,
    localgeo,
    t,
)
    coords = localgeo.coord
    FT = eltype(state)

    # parameters
    _grav::FT = grav(param_set)
    _R_d::FT = R_d(param_set)
    _cv_d::FT = cv_d(param_set)
    _Ω::FT = Omega(param_set)
    _a::FT = planet_radius(param_set)
    _p_0::FT = MSLP(param_set)

    k::FT = 3
    T_E::FT = 310
    T_P::FT = 240
    T_0::FT = 0.5 * (T_E + T_P)
    Γ::FT = 0.005
    A::FT = 1 / Γ
    B::FT = (T_0 - T_P) / T_0 / T_P
    C::FT = 0.5 * (k + 2) * (T_E - T_P) / T_E / T_P
    b::FT = 2
    H::FT = _R_d * T_0 / _grav
    z_t::FT = 15e3
    λ_c::FT = π / 9
    φ_c::FT = 2 * π / 9
    d_0::FT = _a / 6
    V_p::FT = 1
    M_v::FT = 0.608
    p_w::FT = 34e3             ## Pressure width parameter for specific humidity
    η_crit::FT = 10 * _p_0 / p_w ## Critical pressure coordinate
    q_0::FT = 0                ## Maximum specific humidity (default: 0.018)
    q_t::FT = 1e-12            ## Specific humidity above artificial tropopause
    φ_w::FT = 2π / 9           ## Specific humidity latitude wind parameter

    # grid
    λ = @inbounds atan(coords[2], coords[1])
    φ = @inbounds asin(coords[3] / norm(coords, 2))
    z = norm(coords) - _a

    r::FT = z + _a
    γ::FT = 1 # set to 0 for shallow-atmosphere case and to 1 for deep atmosphere case

    # convenience functions for temperature and pressure
    τ_z_1::FT = exp(Γ * z / T_0)
    τ_z_2::FT = 1 - 2 * (z / b / H)^2
    τ_z_3::FT = exp(-(z / b / H)^2)
    τ_1::FT = 1 / T_0 * τ_z_1 + B * τ_z_2 * τ_z_3
    τ_2::FT = C * τ_z_2 * τ_z_3
    τ_int_1::FT = A * (τ_z_1 - 1) + B * z * τ_z_3
    τ_int_2::FT = C * z * τ_z_3
    I_T::FT =
        (cos(φ) * (1 + γ * z / _a))^k -
        k / (k + 2) * (cos(φ) * (1 + γ * z / _a))^(k + 2)

    # base state virtual temperature, pressure, specific humidity, density
    T_v::FT = (τ_1 - τ_2 * I_T)^(-1)
    p::FT = _p_0 * exp(-_grav / _R_d * (τ_int_1 - τ_int_2 * I_T))

    # base state velocity
    U::FT =
        _grav * k / _a *
        τ_int_2 *
        T_v *
        (
            (cos(φ) * (1 + γ * z / _a))^(k - 1) -
            (cos(φ) * (1 + γ * z / _a))^(k + 1)
        )
    u_ref::FT =
        -_Ω * (_a + γ * z) * cos(φ) +
        sqrt((_Ω * (_a + γ * z) * cos(φ))^2 + (_a + γ * z) * cos(φ) * U)
    v_ref::FT = 0
    w_ref::FT = 0

    # velocity perturbations
    F_z::FT = 1 - 3 * (z / z_t)^2 + 2 * (z / z_t)^3
    if z > z_t
        F_z = FT(0)
    end
    d::FT = _a * acos(sin(φ) * sin(φ_c) + cos(φ) * cos(φ_c) * cos(λ - λ_c))
    c3::FT = cos(π * d / 2 / d_0)^3
    s1::FT = sin(π * d / 2 / d_0)
    if 0 < d < d_0 && d != FT(_a * π)
        u′::FT =
            -16 * V_p / 3 / sqrt(3) *
            F_z *
            c3 *
            s1 *
            (-sin(φ_c) * cos(φ) + cos(φ_c) * sin(φ) * cos(λ - λ_c)) /
            sin(d / _a)
        v′::FT =
            16 * V_p / 3 / sqrt(3) * F_z * c3 * s1 * cos(φ_c) * sin(λ - λ_c) /
            sin(d / _a)
    else
        u′ = FT(0)
        v′ = FT(0)
    end
    w′::FT = 0
    u_sphere = SVector{3, FT}(u_ref + u′, v_ref + v′, w_ref + w′)
    #u_sphere = SVector{3, FT}(u_ref, v_ref, w_ref)
    u_cart = sphr_to_cart_vec(u_sphere, φ, λ)

    ## temperature & density
    T::FT = T_v
    ρ::FT = p / (_R_d * T)
    ## potential & kinetic energy
    e_pot = aux.Φ
    e_kin::FT = 0.5 * u_cart' * u_cart
    e_int = _cv_d * T

    ## Assign state variables
    state.ρ = ρ
    state.ρu = ρ * u_cart
    if total_energy
        state.ρe = ρ * (e_int + e_kin + e_pot)
    else
        state.ρe = ρ * (e_int + e_kin)
    end

    nothing
end

function main()
    ClimateMachine.init(parse_clargs = true)
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    polynomialorder = 3
    numelem_horz = 8
    numelem_vert = 5

    timeend = 10 * 24 * 3600
    outputtime = 24 * 3600

    FT = Float64
    result = run(
        mpicomm,
        polynomialorder,
        numelem_horz,
        numelem_vert,
        timeend,
        outputtime,
        ArrayType,
        FT,
    )
end

function run(
    mpicomm,
    polynomialorder,
    numelem_horz,
    numelem_vert,
    timeend,
    outputtime,
    ArrayType,
    FT,
)
    _planet_radius::FT = planet_radius(param_set)
    domain_height = FT(30e3)
    vert_range = grid1d(
        _planet_radius,
        FT(_planet_radius + domain_height),
        nelem = numelem_vert,
    )
    topology = StackedCubedSphereTopology(mpicomm, numelem_horz, vert_range)

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = polynomialorder,
        meshwarp = cubedshellwarp,
    )

    T_profile =
        DecayingTemperatureProfile{FT}(param_set, FT(290), FT(220), FT(8e3))

    if total_energy
        sources = (Coriolis(),)
    else
        sources = (Coriolis(), Gravity())
    end
    problem = BaroclinicWave()
    model = DryAtmosModel{FT}(
        SphericalOrientation(),
        problem,
        ref_state = DryReferenceState(T_profile),
        sources = sources,
    )

    esdg = ESDGModel(
        model,
        grid,
        #volume_numerical_flux_first_order = CentralVolumeFlux(),
        #volume_numerical_flux_first_order = EntropyConservative(),
        volume_numerical_flux_first_order = KGVolumeFlux(),
        #surface_numerical_flux_first_order = MatrixFlux(),
        surface_numerical_flux_first_order = RusanovNumericalFlux(),
    )

    linearmodel = DryAtmosAcousticGravityLinearModel(model)
    lineardg = DGModel(
        linearmodel,
        grid,
        RusanovNumericalFlux(),
        #CentralNumericalFluxFirstOrder(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient();
        direction = VerticalDirection(),
        state_auxiliary = esdg.state_auxiliary,
    )

    # determine the time step
    element_size = (domain_height / numelem_vert)
    acoustic_speed = soundspeed_air(param_set, FT(330))

    dx = min_node_distance(grid)
    cfl = 3
    dt = cfl * dx / acoustic_speed

    Q = init_ode_state(esdg, FT(0))

    #odesolver = LSRK144NiegemannDiehlBusch(esdg, Q; dt = dt, t0 = 0)

    linearsolver = ManyColumnLU()
    odesolver = ARK2GiraldoKellyConstantinescu(
        esdg,
        lineardg,
        LinearBackwardEulerSolver(linearsolver; isadjustable = false),
        Q;
        dt = dt,
        t0 = 0,
        split_explicit_implicit = false,
    )

    eng0 = norm(Q)
    @info @sprintf """Starting
                      ArrayType       = %s
                      FT              = %s
                      polynomialorder = %d
                      numelem_horz    = %d
                      numelem_vert    = %d
                      dt              = %.16e
                      norm(Q₀)        = %.16e
                      """ "$ArrayType" "$FT" polynomialorder numelem_horz numelem_vert dt eng0

    # Set up the information callback
    starttime = Ref(now())
    cbinfo = EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            runtime = Dates.format(
                convert(DateTime, now() - starttime[]),
                dateformat"HH:MM:SS",
            )
            @info @sprintf """Update
                              simtime = %.16e
                              runtime = %s
                              norm(Q) = %.16e
                              """ gettime(odesolver) runtime energy
        end
    end
    cbcfl = EveryXSimulationSteps(100) do
        simtime = gettime(odesolver)

        @views begin
            ρ = Array(Q.data[:, 1, :])
            ρu = Array(Q.data[:, 2, :])
            ρv = Array(Q.data[:, 3, :])
            ρw = Array(Q.data[:, 4, :])
        end

        u = ρu ./ ρ
        v = ρv ./ ρ
        w = ρw ./ ρ

        ue = extrema(u)
        ve = extrema(v)
        we = extrema(w)

        @info @sprintf """CFL
                simtime = %.16e
                u = (%.4e, %.4e)
                v = (%.4e, %.4e)
                w = (%.4e, %.4e)
                """ simtime ue... ve... we...
    end
    callbacks = (cbinfo, cbcfl)

    #filterorder = 32
    #filter = ExponentialFilter(grid, 0, filterorder)
    #cbfilter = EveryXSimulationSteps(1) do
    #    Filters.apply!(
    #        Q,
    #        #AtmosFilterPerturbations(model),
    #        :,
    #        grid,
    #        filter,
    #       # state_auxiliary = esdg.state_auxiliary,
    #    )
    #    nothing
    #end
    #callbacks = (callbacks..., cbfilter)

    if output_vtk
        # create vtk dir
        vtkdir =
            "vtk_esdg_total_KG_ncg_hires_baroclinic" *
            "_poly$(polynomialorder)_horz$(numelem_horz)_vert$(numelem_vert)" *
            "_$(ArrayType)_$(FT)"
        mkpath(vtkdir)

        vtkstep = 0
        # output initial step
        do_output(mpicomm, vtkdir, vtkstep, esdg, Q, model)

        # setup the output callback
        cbvtk = EveryXSimulationSteps(floor(outputtime / dt)) do
            vtkstep += 1
            do_output(mpicomm, vtkdir, vtkstep, esdg, Q, model)
        end
        callbacks = (callbacks..., cbvtk)
    end

    ## Create a callback to report state statistics for main MPIStateArrays
    ## every ntFreq timesteps.
    nt_freq = floor(Int, 1 // 10 * timeend / dt)
    cbsc =
        ClimateMachine.StateCheck.sccreate([(Q, "state")], nt_freq; prec = 12)
    callbacks = (callbacks..., cbsc)

    solve!(
        Q,
        odesolver;
        timeend = timeend,
        adjustfinalstep = false,
        callbacks = callbacks,
    )

    # final statistics
    engf = norm(Q)
    @info @sprintf """Finished
    norm(Q)                 = %.16e
    norm(Q) / norm(Q₀)      = %.16e
    norm(Q) - norm(Q₀)      = %.16e
    """ engf engf / eng0 engf - eng0
    engf

    ## Check results against reference if present
    ClimateMachine.StateCheck.scprintref(cbsc)
    #! format: off
    refDat = (
        [
            [ "state",     "ρ",   1.26085994139264381125e-02,  1.51502814805562069367e+00,  3.46330071392505767225e-01,  3.42221017037761976454e-01 ],
            [ "state", "ρu[1]",  -1.26747868050267172180e+02,  1.22735648852872344605e+02, -6.56484249582622303443e-02,  1.07588365672914481053e+01 ],
            [ "state", "ρu[2]",  -1.44635478251794808102e+02,  1.11383888659731638882e+02, -1.53109264073002888581e-03,  1.06376480052955262323e+01 ],
            [ "state", "ρu[3]",  -1.50479775987282266669e+02,  1.62284843398145170568e+02,  5.84596035289077428643e-02,  9.73728925076320095400e+00 ],
            [ "state",    "ρe",   1.73490214503025617887e+03,  2.70352534694924892392e+05,  6.42401799036320589948e+04,  7.10054852130306971958e+04 ],
        ],
        [
            [ "state",     "ρ",    12,    12,    12,    12 ],
            [ "state", "ρu[1]",    12,    12,    12,    12 ],
            [ "state", "ρu[2]",    12,    12,    12,    12 ],
            [ "state", "ρu[3]",    12,    12,    12,    12 ],
            [ "state",    "ρe",    12,    12,    12,    12 ],
        ],
    )
    #! format: on
    if length(refDat) > 0
        @test ClimateMachine.StateCheck.scdocheck(cbsc, refDat)
    end
end

function do_output(
    mpicomm,
    vtkdir,
    vtkstep,
    dg,
    Q,
    model,
    testname = "baroclinicwave",
)
    ## name of the file that this MPI rank will write
    filename = @sprintf(
        "%s/%s_mpirank%04d_step%04d",
        vtkdir,
        testname,
        MPI.Comm_rank(mpicomm),
        vtkstep
    )

    statenames = flattenednames(vars_state(model, Prognostic(), eltype(Q)))
    auxnames = flattenednames(vars_state(model, Auxiliary(), eltype(Q)))
    writevtk(
        filename,
        Q,
        dg,
        statenames,
        dg.state_auxiliary,
        auxnames;
        number_sample_points = 10,
    )

    ## Generate the pvtu file for these vtk files
    if MPI.Comm_rank(mpicomm) == 0
        ## name of the pvtu file
        pvtuprefix = @sprintf("%s/%s_step%04d", vtkdir, testname, vtkstep)

        ## name of each of the ranks vtk files
        prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
            @sprintf("%s_mpirank%04d_step%04d", testname, i - 1, vtkstep)
        end

        writepvtu(pvtuprefix, prefixes, (statenames..., auxnames...), eltype(Q))

        @info "Done writing VTK: $pvtuprefix"
    end
end

@testset "$(@__FILE__)" begin
    tic = Base.time()

    main()

    toc = Base.time()
    time = toc - tic
    println(time)
end