text:
using Test
using MiniTB
using TBComponents

# @test size( ps ) == 4
function getnet( data :: Dict, gross :: Float64 ) :: Float64
     person = data[:person]
     person.wage = gross
     person.hours = gross/DEFAULT_WAGE
     rc = calculate( person, data[:params] )
     return rc[:netincome]
end

function locmakebc( person :: Person, params :: Parameters ) :: BudgetConstraint
    data = Dict(
        :person=>person,
        :params=>params )
    bc = makebc( data, getnet )
    return bc
end

p1 = DEFAULT_PERSON
t1 = deepcopy(ZERO_PARAMS)
p1.wage = 100.0
t1.it_rate = [0.2,0.4]
t1.it_band = [50.0,99999]

@testset begin

    tax1 = calculatetax(p1, t1 )
    @test tax1 ≈ 30
    p1.wage += 1
    tax2 = calculatetax(p1, t1 )
    @test tax2-tax1 ≈ 0.4 # mr 0.4

    bc = locmakebc( p1, t1 )

    println( "bc=$bc" )
    @test size( bc )[1] in [3,4] # allow myself a spare...
end


text:
# Types
export DeterministicStrategy

# Type validation methods
export is_deterministic

# constructor methods
export deterministic_strategies

"""
    DeterministicStrategy(conditionals :: Matrix{Int64}) <: AbstractStrategy{Int64}

A strategy matrix describing the behavior of a deterministic black-box.
A strategy *deterministic* if its elements satisfy ``P(y|x)\\in\\{0,1\\}`` in addition
to the non-negativity and normalization constraints.
By default, the constructor creates a strategy for a 'BlackBox' scenario, however,
a `Scenario` can be passed to the `DeterministicStrategy` constructor.

    DeterministicStrategy(conditionals :: Matrix{Int}, scenario :: Scenario)

The product of two deterministic strategies is a `DeterministicStrategy`.

    *(S1::DeterministicStrategy, S2::DeterministicStrategy) :: DeterministicStrategy

When multiplied a new `Scenario` can be specified.

    *(
        S1::DeterministicStrategy,
        S2::DeterministicStrategy,
        scenario::Scenario
    ) :: Deterministic Strategy

### Errors:

A `DomainError` is thrown if:
* The provided `Scenario` does not match the dimension of the `conditionals` matrix.
* The elements of `conditionals` are not `0` or `1`.
* The strategy elements are not non-negative and normalized.
"""
struct DeterministicStrategy <: AbstractStrategy{Int64}
    conditionals :: Matrix{Int64}
    scenario :: Scenario
    DeterministicStrategy(
        conditionals :: Matrix{<:Real},
        scenario :: Scenario
    ) = (is_deterministic(conditionals) && (size(conditionals) == strategy_dims(scenario))) ? new(conditionals, scenario) : throw(
        DomainError(conditionals, "conditionals are not a deterministic (`{0,1}`) distribution.")
    )
    DeterministicStrategy(
        conditionals :: Matrix{<:Real}
    ) = is_deterministic(conditionals) ? new(conditionals, BlackBox(size(conditionals)...)) : throw(
        DomainError(conditionals, "conditionals are not a deterministic (`{0,1}`) distribution.")
    )
end

*(
    S1::DeterministicStrategy,
    S2::DeterministicStrategy
) = DeterministicStrategy(S1.conditionals*S2.conditionals)
*(
    S1::DeterministicStrategy,
    S2::DeterministicStrategy,
    scenario::Scenario
) = DeterministicStrategy(S1.conditionals*S2.conditionals, scenario)

"""
    is_deterministic( strategy :: AbstractMatrix  ) :: Bool

Returns `true` if all elements of `strategy` are either `0` or `1` and the matrix
is a valid conditional probability distribution.
"""
function is_deterministic(S :: AbstractMatrix) :: Bool
    is_conditionals = is_conditional_distribution(S)
    is_01 = all(i -> (S[i] == 0) || (S[i] == 1), 1:length(S))

    (is_conditionals && is_01)
end

"""
    deterministic_strategies(scenario :: BlackBox) :: Vector{Matrix{Int64}}

    deterministic_strategies(num_out :: Int64, num_in :: Int64) :: Vector{Matrix{Int64}}

Enumerates the set of deterministic strategies for the specified `BlackBox`. For
performance, enumerated deterministic strategies are left as matrices and not
constructed into [`DeterministicStrategy`](@ref) types.
"""
deterministic_strategies(
    num_out :: Int64,
    num_in :: Int64
) = deterministic_strategies(BlackBox(num_out,num_in))

function deterministic_strategies(scenario :: BlackBox) :: Vector{Matrix{Int64}}
    num_in = scenario.num_in
    num_out = scenario.num_out

    strategies = (num_out == 1) ? Matrix{Int64}[fill(1,(1,num_in))] : map( i -> begin
        m = zeros(Int64, num_out, num_in)
        base_n_array = digits(i, base = num_out, pad = num_in) .+ 1

        for j in 1:num_in
            m[base_n_array[j],j] = 1
        end

        m
    end, 0:(num_out^num_in - 1))

    strategies
end

ext: jl


text:
using StructuredSolvers
using Test
using Random

@testset "StructuredSolvers.jl" begin
    @test 1 == 1
    # Write your tests here.

    @testset "Conjugate gradient" begin
        n = 5

        M = l1Manifold(ones(n))
        x = zeros(n)

        Random.seed!(1435)
        A = rand(n, n)
        A += A'
        A += Diagonal((1-minimum(eigvals(A))) .* ones(n))
        b = A * ones(n)

        Aoperator(η) = A * η

        ## Solves A * η + b = 0
        d, j, d_type = StructuredSolvers.solve_tCG(M, x, b, Aoperator; ν=1e-13, ϵ_residual = 1e-13, maxiter=1e5, safe_projection=false, printlev=0)

        @test isapprox(d, -ones(n))
        @test d_type == :Solved
    end

end


text:
module Day1
export solvePart1, solvePart2

import Base.Iterators
import DelimitedFiles

function getInput()
    open("day1input.txt") do f
        return DelimitedFiles.readdlm(f, Int64)
    end
end

function solvePart1()
    pos = sum(getInput())
    println("Final resulting frequency: $pos")
    return pos
end

function solvePart2()
    visited = Set{Int64}()
    pos = 0
    push!(visited, pos)

    for input in Base.Iterators.cycle(getInput())
        pos += input
        if in(pos, visited)
            break
        end
        push!(visited, pos)
    end

    println("First frequency reached twice is: $pos")
    return pos
end

end


text:
module Jaynes

using Reexport

import Base.display
import Base: getindex, haskey, iterate, isempty, convert, collect, getindex, setindex!, push!, merge, merge!, get, filter, length, ndims, keys, +, rand, size
import Base: isless
import Base: Pair

# Jaynes implements the abstract GFI from Gen.
import Gen
import Gen: Selection, ChoiceMap, Trace, GenerativeFunction
import Gen: DynamicChoiceMap, EmptySelection
import Gen: get_value, has_value
import Gen: get_values_shallow, get_submaps_shallow
import Gen: get_args, get_retval, get_choices, get_score, get_gen_fn, has_argument_grads, accepts_output_grad, get_params
import Gen: select, choicemap
import Gen: simulate, generate, project, propose, assess, update, regenerate
import Gen: init_param!, accumulate_param_gradients!, choice_gradients, init_update_state, apply_update!

# Gen diff types.
import Gen: Diff, UnknownChange, NoChange
import Gen: SetDiff, DictDiff, VectorDiff, IntDiff, Diffed

# Yarrrr I'm a com-pirate!
using MacroTools
using IRTools
using IRTools: @dynamo, IR, xcall, arguments, insertafter!, recurse!, isexpr, self, argument!, Variable, meta, renumber, Pipe, finish, blocks, predecessors, dominators, block, successors, Block, block!, branches, Branch, branch!, CFG, stmt
using Random
using Mjolnir: inline_consts!, partials!, ssa!, prune!
using InteractiveUtils: subtypes

# Static selektor.
using StaticArrays

@reexport using Distributions
import Distributions: Distribution
import Distributions: logpdf

# Differentiable.
@reexport using Zygote
using ForwardDiff
using ForwardDiff: Dual
using DistributionsAD

# Plotting.
using UnicodePlots: lineplot

# ------------ Toplevel importants ------------ #

const Address = Union{Int, Symbol, Pair}

# This is primarily used when mapping choice maps to arrays.
isless(::Symbol, ::Pair) = true
isless(::Pair, ::Symbol) = false
isless(::Int, ::Symbol) = true
isless(::Symbol, ::Int) = false
isless(::Int, ::Pair) = true
isless(::Pair, ::Int) = false

include("unwrap.jl")
include("kludges.jl")

# Jaynes introduces a new type of generative function.
abstract type TypedGenerativeFunction{N, R, Tr, T} <: GenerativeFunction{R, Tr} end

# ------------ includes ------------ #

include("core.jl")
export trace

include("compiler.jl")
export NoChange, Change, ScalarDiff, IntDiff, DictDiff, SetDiff, VectorDiff, BoolDiff, Diffed
export Δ, _propagate
export prepare_ir!, infer!, @abstract, InterpretationContext
export detect_switches, detect_kernel

include("pipelines.jl")
export NoStatic, DefaultPipeline, StaticWithLints, SpecializerPipeline, AutomaticAddressingPipeline
export record_cached!

include("macros.jl")
export @primitive, @jaynes

# Tracer language features.
export learnable, fillable, factor

# Compiler.
export NoChange, UndefinedChange
export construct_graph, compile_function

# Vectors to dynamic value address map.
export dynamic

# Selections and parameters.
export select, target, static, array, learnables
export anywhere, intersection, union
export compare, update_learnables, merge!, merge

# Gen compat.
include("gen_fn_interface.jl")

export JFunction, JTrace
export get_analysis, get_ir, get_fn
export init_param!, accumulate_param_gradients!, choice_gradients
export choicemap, select
export get_value, has_value
export get_params_grad

# Contexts.
export generate
export simulate
export update
export propose
export regenerate
export assess
export get_learnable_gradients, get_choice_gradients
export get_learnable_gradient, get_choice_gradient

constrain(v::Vector{Pair{T, K}}) where {T <: Tuple, K} = JChoiceMap(target(v))
export constrain

# Typing rules.
include("typing_rules.jl")
export absint

# Utilities.
export display, getindex, haskey, get_score, get_ret, flatten, lineplot

end # module


text:
"""
    run_two_step_clust(data::ClustData; norm_op::String="zscore", norm_scope::String="full", method::String="kmeans", representation::String="centroid", n_clust_1::Int=5, n_clust_2::Int=3, n_seg::Int=data.T, n_init::Int=100, iterations::Int=300, attribute_weights::Dict{String,Float64}=Dict{String,Float64}(), features_2::Array{String,1}=Array{String,1}(), get_all_clust_results::Bool=false, kwargs...)
Run a first regular cluster over the data and a second about the defined `features_2` within each period.
"""
function run_two_step_clust(data::ClustData;
                            norm_op::String="zscore",
                            norm_scope::String="full",
                            method::String="kmeans",
                            representation::String="centroid",
                            n_clust_1::Int=5,
                            n_clust_2::Int=3,
                            n_seg::Int=data.T,
                            n_init::Int=100,
                            iterations::Int=300,
                            attribute_weights::Dict{String,Float64}=Dict{String,Float64}(),
                            features_2::Array{String,1}=Array{String,1}(),
                            get_all_clust_results::Bool=false,
                            kwargs...)
  clust_data=run_clust(data;norm_op=norm_op,norm_scope=norm_scope,method=method,representation=representation,n_clust=n_clust_1,n_init=n_init,iterations=iterations,attribute_weights=attribute_weights).best_results
  new_k_ids=deepcopy(clust_data.k_ids)
  for i in 1:clust_data.K
      index=findall(clust_data.k_ids.==i)
      dt=Dict{String,Array}()
      for name in keys(data.data)
        tech=split(name,"-")[1]
        if name in features_2 || tech in features_2
          dt[name]=(data.data[name][:,index])
        end
      end
      single_period=ClustData(FullInputData(clust_data.region, clust_data.years, 0, dt),Int64(clust_data.weights[i]),clust_data.T)
      new_k_ids[index]=run_clust(single_period; n_clust=Int64(min(single_period.K,n_clust_2)),n_init=Int64(round(n_init/clust_data.K))).best_ids.+(i-1)*n_clust_2
  end
  return run_clust(data;method="predefined",representation=representation,n_init=1,n_clust=Int64(n_clust_1*n_clust_2),k_ids=new_k_ids,n_seg=n_seg,get_all_clust_results=get_all_clust_results)
end

"""
    run_minmax_clust(data::ClustData; norm_op::String="zscore", norm_scope::String="full", method::String="kmeans", representation::String="centroid", n_clust_1::Int=5, n_clust_2::Int=3, n_seg::Int=data.T, n_init::Int=100, iterations::Int=300, attribute_weights::Dict{String,Float64}=Dict{String,Float64}(), min::Array{String,1}=Array{String,1}(), max::Array{String,1}=Array{String,1}(),, get_all_clust_results::Bool=false, kwargs...)
Choose the min or max of a certain attribute to represent this attribute in this period instead of the centroid or so
"""
function run_minmax_clust(data::ClustData;
                            norm_op::String="zscore",
                            norm_scope::String="full",
                            method::String="kmeans",
                            representation::String="centroid",
                            n_clust::Int=5,
                            n_seg::Int=data.T,
                            n_init::Int=100,
                            iterations::Int=300,
                            attribute_weights::Dict{String,Float64}=Dict{String,Float64}(),
                            min::Array{String,1}=Array{String,1}(),
                            max::Array{String,1}=Array{String,1}(),
                            get_all_clust_results::Bool=false,
                            kwargs...)
  clust_result=run_clust(data;norm_op=norm_op,norm_scope=norm_scope,method=method,representation=representation,n_clust=n_clust,n_init=n_init,iterations=iterations,attribute_weights=attribute_weights)
  clust_data=clust_result.best_results
  for i in 1:clust_data.K
    index=findall(clust_data.k_ids.==i)
    for name in keys(data.data)
      att=split(name,"-")[1]
      if name in min || att in min
        indexminmax=index[findmin(mean(data.data[name][:,index],dims=1))[2][2]]
        clust_data.data[name][:,i]=data.data[name][:,indexminmax]
      elseif name in max || att in max
        indexminmax=index[findmax(mean(data.data[name][:,index],dims=1))[2][2]]
        clust_data.data[name][:,i]=data.data[name][:,indexminmax]
      end
    end
  end
  return clust_result
end

"""
    run_pure_clust(data::ClustData; norm_op::String="zscore", norm_scope::String="full", method::String="kmeans", representation::String="centroid", n_clust_1::Int=5, n_clust_2::Int=3, n_seg::Int=data.T, n_init::Int=100, iterations::Int=300, attribute_weights::Dict{String,Float64}=Dict{String,Float64}(), clust::Array{String,1}=Array{String,1}(), get_all_clust_results::Bool=false, kwargs...)
Replace the original timeseries of the attributes in clust with their clustered value
"""
function run_pure_clust(data::ClustData;
                            norm_op::String="zscore",
                            norm_scope::String="full",
                            method::String="kmeans",
                            representation::String="centroid",
                            n_clust::Int=5,
                            n_seg::Int=data.T,
                            n_init::Int=100,
                            iterations::Int=300,
                            attribute_weights::Dict{String,Float64}=Dict{String,Float64}(),
                            clust::Array{String,1}=Array{String,1}(),
                            get_all_clust_results::Bool=false,
                            kwargs...)
  clust_result=run_clust(data;norm_op=norm_op,norm_scope=norm_scope,method=method,representation=representation,n_clust=n_clust,n_init=n_init,iterations=iterations,attribute_weights=attribute_weights)
  clust_data=clust_result.best_results
  mod_data=deepcopy(data.data)
  for i in 1:clust_data.K
    index=findall(clust_data.k_ids.==i)
    for name in keys(mod_data)
      att=split(name,"-")[1]
      if name in clust || att in clust
        mod_data[name][:,index]=repeat(clust_data.data[name][:,i], outer=(1,length(index)))
      end
    end
  end
  return ClustResultSimple(ClustData(data.region, data.years, data.K, data.T, mod_data, data.weights, data.deltas, data.k_ids), clust_result.clust_config)
end

"""
    run_darkperiod_clust(data::ClustData, ev_arr::Array{SimpleExtremeValueDescr,1}; norm_op::String="zscore", norm_scope::String="full", method::String="kmeans", representation::String="centroid", n_clust_1::Int=5, n_clust_2::Int=3, n_seg::Int=data.T, n_init::Int=100, iterations::Int=300, attribute_weights::Dict{String,Float64}=Dict{String,Float64}(), get_all_clust_results::Bool=false, kwargs...)
Choose Seasonal extremes based on the provided Simple Extreme Value Descriptions
"""
function run_darkperiod_clust(data::ClustData,
                                ev_arr::Array{SimpleExtremeValueDescr,1};
                                norm_op::String="zscore",
                                norm_scope::String="full",
                                method::String="kmeans",
                                representation::String="centroid",
                                n_clust::Int=5,
                                n_clust_dark::Int=5,
                                n_seg::Int=data.T,
                                n_init::Int=100,
                                n_dark_init::Int=100,
                                iterations::Int=300,
                                attribute_weights::Dict{String,Float64}=Dict{String,Float64}(),
                                get_all_clust_results::Bool=false,
                                config::String="",
                                kwargs...)
    extr_results=Array{ClustData,1}()
    for ev in ev_arr
        data_mod,extr_val,x = simple_extr_val_sel(data,ev;rep_mod_method="append")
        if config=="2extr"
            ev=SimpleExtremeValueDescr(ev.data_type,ev.extremum,ev.peak_def,n_clust_dark)
            x,extr,x = simple_extr_val_sel(extr_val,ev;rep_mod_method="append")
            k_ids=zeros(Int64,size(extr_val.k_ids))
            k_ids[findall(extr_val.k_ids.!=0)].= 1
            extr_result=ClustData(extr.region, extr.years, extr.K, extr.T, extr.data, extr.weights, extr.deltas, k_ids)
        else
            extr_result=run_clust(extr_val;norm_op=norm_op,norm_scope=norm_scope,method=method,representation=representation,n_clust=n_clust_dark,n_init=n_dark_init,iterations=iterations,attribute_weights=attribute_weights,kwargs...).best_results
        end
        push!(extr_results,extr_result)
        data=data_mod
    end
    # run clustering
    data_res = run_clust(data;norm_op=norm_op,norm_scope=norm_scope,method=method,representation=representation,n_clust=n_clust,n_init=n_init,iterations=iterations,attribute_weights=attribute_weights) # default
    # representation modification
    return ClustResultSimple(representation_modification(extr_results, data_res.best_results), data_res.clust_config)
end


text:
### visualize a solution.

function visualizefiltersolution(𝓧, xq, wfunc, dfunc,f)
    𝑤_new = get𝑤scaled(𝓧)
    h_new = geth(𝑤_new, wfunc, dfunc, 𝓧)

    f_itp = xx->getp(xx, 𝑤_new, h_new, wfunc, dfunc, 𝓧)
    f_itp_xq = collect( f_itp(xq[i]) for i = 1:length(xq) )
    f_itp_𝓧 = collect( f_itp(𝓧[i]) for i = 1:length(𝓧) )

    𝑒 = xx->(f(xx)-f_itp(xx))
    𝑒_xq = collect( 𝑒(xq[i]) for i = 1:length(xq) )
    𝑒_𝓧 = collect( 𝑒(𝓧[i]) for i = 1:length(𝓧) )

    return f_itp_xq, f_itp_𝓧, 𝑒_xq, 𝑒_𝓧
end

function packagefiltersolutionfast( 𝓧_new,
                                    h_history,
                                    wfunc,
                                    zfunc,
                                    dfunc,
                                    filter_type_num::Int)

##
    c,ν, discrepancy, xq, yq_bary, yq_cheby = BarycentrictoChebyshevfast(𝓧_new, wfunc, dfunc, filter_type_num)
    println("conversion discrepany between barycenter and Chebyshev is ", discrepancy)

    ######
    L = length(c)

    L_old = L

    # get impulse response.
    𝑏 = c
    if filter_type_num == 2
        𝑏 = converttotype2(c)
    elseif filter_type_num == 3
        𝑏 = converttotype3(c)
    elseif filter_type_num ==4
        𝑏 = converttotype4(c)
    end

    q = 1
    if filter_type_num == 1 || filter_type_num == 2
        q = 0
    end

    gethfunc = gethtype4
    if filter_type_num == 1
        L = length(c)-1

        gethfunc = gethtype1
    elseif filter_type_num == 2
        gethfunc = gethtype2
    elseif filter_type_num ==3
        gethfunc = gethtype3
    end

    h = gethfunc(𝑏)

    return h
end

function packagefiltersolution(𝓧_new,
                                h_history,
                                wfunc,
                                zfunc,
                                dfunc,
                                filter_type_num::Int,
                                fig_num::Int,
                                savefig_flag::Bool,
                                plot_flag::Bool)

    @time c,ν, discrepancy, xq, yq_bary, yq_cheby = BarycentrictoChebyshev(𝓧_new, wfunc, dfunc, filter_type_num)
    println("conversion discrepany between barycenter and Chebyshev is ", discrepancy)

    ######
    L = length(c)

    L_old = L

    # get impulse response.
    𝑏 = c
    if filter_type_num == 2
        𝑏 = converttotype2(c)
    elseif filter_type_num == 3
        𝑏 = converttotype3(c)
    elseif filter_type_num ==4
        𝑏 = converttotype4(c)
    end

    q = 1
    if filter_type_num == 1 || filter_type_num == 2
        q = 0
    end

    gethfunc = gethtype4
    if filter_type_num == 1
        L = length(c)-1

        gethfunc = gethtype1
    elseif filter_type_num == 2
        gethfunc = gethtype2
    elseif filter_type_num ==3
        gethfunc = gethtype3
    end

    h = gethfunc(𝑏)

    if !plot_flag
        return h, fig_num
    end

    ######### plot to interpret results.

    PyPlot.figure(fig_num)
    fig_num += 1
    PyPlot.plot(xq, yq_bary, "--", label = "yq_bary")
    PyPlot.plot(xq, yq_cheby, label = "yq_cheby")
    PyPlot.title("Barycenter vs. Chebyshev")
    PyPlot.legend()

    if savefig_flag
        save_name = Printf.@sprintf("./outputs/type_%s_Barycenter_to_Chebyshev.png",
                                        filter_type_num)
        PyPlot.savefig(save_name)
        sleep(save_delay)
        PyPlot.close(fig_num)
    end

    ## plot level error history.

    PyPlot.figure(fig_num)
    fig_num += 1

    PyPlot.plot(collect(1:length(h_history)), abs.(h_history), ".-")
    PyPlot.title("level error history")

    if savefig_flag
        save_name = Printf.@sprintf("./outputs/type_%s_level_error_history.png",
                                        filter_type_num)
        PyPlot.savefig(save_name)
        sleep(save_delay)
        PyPlot.close(fig_num)
    end

    ## plot resultant filter's response.
    f(xx) = wfunc(xx)*dfunc(xx)
    L = length(𝓧_new)-2
    xq = convert(Vector{BigFloat},collect(LinRange(-1,1,5*L)))

    f_itp_xq, f_itp_𝓧_new, 𝑒_xq, 𝑒_𝓧_new = visualizefiltersolution(𝓧_new,xq,wfunc,dfunc,f)

    PyPlot.figure(fig_num)
    fig_num += 1
    PyPlot.plot(xq, f.(xq), "--", label = "f = w*d")
    PyPlot.plot(xq, f_itp_xq, label = "solution")
    PyPlot.plot(xq, dfunc.(xq), label = "d")
    PyPlot.plot(𝓧_new, f_itp_𝓧_new,".", label = "reference positions")
    PyPlot.title("solution filter's response")
    PyPlot.legend()

    if savefig_flag
        save_name = Printf.@sprintf("./outputs/type_%s_solution.png",
                                        filter_type_num)
        PyPlot.savefig(save_name)
        sleep(save_delay)
        PyPlot.close(fig_num)
    end

    PyPlot.figure(fig_num)
    fig_num += 1
    PyPlot.plot(xq, 𝑒_xq, label = "error")
    PyPlot.plot(𝓧_new, 𝑒_𝓧_new,".", label = "reference positions")
    PyPlot.title("error")
    PyPlot.legend()

    if savefig_flag
        save_name = Printf.@sprintf("./outputs/type_%s_solution_error.png",
                                        filter_type_num)
        PyPlot.savefig(save_name)
        sleep(save_delay)
        PyPlot.close(fig_num)
    end

    ###### analyze

    # prepare.
    N = length(h)

    a = 0.0 #-π
    b = 2*π # π
    ω_set = convert(Vector{BigFloat}, collect( LinRange(a,b,N*20) ) )
    ω_set_fft = collect( LinRange(a,b-2*π/N,N))

    # this is the verified one.
    y = Vector{Complex{BigFloat}}(undef,0)
    y2 = Vector{BigFloat}(undef,0)

    if filter_type_num == 1
        y = collect( evalunivariateChebyshevpolynomial(c,cos(ω_set[i]))*exp( im*(q*π - (N-1)*ω_set[i])/2 ) for i = 1:length(ω_set) )
        y2 = collect( sum( c[k+1]*cos(k*ω_set[i]) for k = 0:L ) for i = 1:length(ω_set) )

    elseif filter_type_num == 2
        y = collect( cos(ω_set[i]/2)*evalunivariateChebyshevpolynomial(c,cos(ω_set[i]))*exp( im*(q*π - (N-1)*ω_set[i])/2 ) for i = 1:length(ω_set) )
        y2 = collect( sum( 𝑏[k]*cos((k-0.5)*ω_set[i]) for k = 1:L ) for i = 1:length(ω_set) )

    elseif filter_type_num == 3
       y = collect( sin(ω_set[i])*evalunivariateChebyshevpolynomial(c,cos(ω_set[i]))*exp( im*(q*π - (N-1)*ω_set[i])/2 ) for i = 1:length(ω_set) )
       y2 = collect( sum( 𝑏[k]*sin(k*ω_set[i]) for k = 1:L ) for i = 1:length(ω_set) )

    elseif filter_type_num == 4
       y = collect( sin(ω_set[i]/2)*evalunivariateChebyshevpolynomial(c,cos(ω_set[i]))*exp( im*(q*π - (N-1)*ω_set[i])/2 ) for i = 1:length(ω_set) )
       y2 = collect( sum( 𝑏[k]*sin((k-0.5)*ω_set[i]) for k = 1:L ) for i = 1:length(ω_set) )
    end

    y = abs.(y)
    y2 = abs.(y2)
    discrepancy = LinearAlgebra.norm(y-y2)
    println("conversion discrepany between y and y2 is ", discrepancy)

    ##### Visualize impulse response.
    PyPlot.figure(fig_num)
    fig_num += 1
    PyPlot.plot(h, label = "h")
    PyPlot.plot(h, ".")
    title_string = Printf.@sprintf("check if the recovered h has the characteristics of type %d", filter_type_num)
    PyPlot.title(title_string)
    PyPlot.legend()

    if savefig_flag
        save_name = Printf.@sprintf("./outputs/type_%s_solution_impulse_rsp.png",
                                        filter_type_num)
        PyPlot.savefig(save_name)
        sleep(save_delay)
        PyPlot.close(fig_num)
    end

    ##### prepare fft and DFT from impulse response.

    impulse_rsp = fft(convert(Vector{Float64},h))
    mag_rsp_fft = abs.(impulse_rsp)
    phase_rsp_fft = angle.(impulse_rsp)

    impulse_rsp_G = collect( naiveimpulsersp(h,ω_set[i]) for i = 1:length(ω_set) )
    mag_rsp_G = abs.(impulse_rsp_G)
    phase_rsp_G = angle.(impulse_rsp_G)

    println("discrepancy between y and G is ", LinearAlgebra.norm(y-mag_rsp_G))

    ##### visualize.

    ## plot the polynomial specified by b, see if it is in agreement with
    #   the fft/DFT response.
    PyPlot.figure(fig_num)
    fig_num += 1
    PyPlot.plot(ω_set_fft, mag_rsp_fft, label = "fft: solution discrete")
    PyPlot.plot(ω_set_fft, mag_rsp_fft, ".")

    PyPlot.plot(ω_set, mag_rsp_G, label = "G: solution continuous")

    y_d = abs.(wfunc.(cos.(ω_set)).*dfunc.(cos.(ω_set)))
    PyPlot.plot(ω_set, y_d, label = "y_d: target")

    # # same as y_d but without the abs().
    # zq = convert(Vector{BigFloat}, collect( LinRange(0,π,5*L) ) )
    # z = wfunc.(cos.(zq)).*zfunc.(zq)
    # PyPlot.plot(zq, z, "--", label = "(w∘cos).*targetfunc")

    PyPlot.title("Magnitude response")
    PyPlot.legend()

    if savefig_flag
        save_name = Printf.@sprintf("./outputs/type_%s_solution_magnitude_rsp.png",
                                        filter_type_num)
        PyPlot.savefig(save_name)
        sleep(save_delay)
        PyPlot.close(fig_num)
    end

    return h, fig_num
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
## ----------------------------------------------------------------------------------
function collect_ME_DAT!(DAT, iJR; nths = nthreads())

    project = string(parentmodule(iJR))

    costider = iJR.COST_IDER
    objider = iJR.BIOMASS_IDER

    rxns_map = iJR.load_rxns_map()

    # util stuff
    WLOCK = ReentrantLock()
    isexpdep(method) = !(method in DAT[:EXP_INDEPENDENT])
    depks(method, typ, Ed_met, exp) =
        isexpdep(method) ? (method, typ, Ed_met, exp) : (method, typ, Ed_met)

    function dat_file(;method, exp)
        kwargs = isexpdep(method) ? (;method, exp) : (;method)
        procdir(iJR, "maxent_ep_dat", kwargs, ".jls")
    end

    # Feed jobs
    Ch = Channel(nths) do ch
        for method in DAT[:ME_METHODS]
            for exp in DAT[:EXPS]
                put!(ch, (exp, method))
                !isexpdep(method) && break # Do only once
            end
        end
    end

    @threads for thid in 1:nths
        thid = threadid()
        for (exp, method) in Ch

            # ME data
            datfile = dat_file(;method, exp)
            # !isfile(datfile) && continue
            dat = ldat(datfile)

            model = dat[:model]
            epouts = dat[:epouts]
            exp_beta = dat[:exp_beta] # maximum(keys(epouts))
            epout = epouts[exp_beta]

            lock(WLOCK) do
                @info("Doing",
                    project,
                    exp, method,
                    length(epouts),
                    epout.iter, thid
                ); println()
            end

            # Biomass and cost
            ep_biom = ChU.av(model, epout, objider)
            ep_biomstd = sqrt(ChU.va(model, epout, objider))
            ep_cost = ChU.av(model, epout, costider)
            ep_coststd = sqrt(ChU.va(model, epout, costider))

            # store
            lock(WLOCK) do
                DAT[depks(method, :flx, "D", exp)...] = ep_biom
                DAT[depks(method, :err, "D", exp)...] = ep_biomstd
                DAT[depks(method, :flx, "cost", exp)...] = ep_cost
                DAT[depks(method, :err, "cost", exp)...] = ep_coststd
            end

            # exchanges
            for Ed_met in DAT[:EXCH_FLX_IDERS]
                model_exch = rxns_map[Ed_met]

                # flxs
                ep_av = ChU.av(model, epout, model_exch)
                ep_std = sqrt(ChU.va(model, epout, model_exch))

                lock(WLOCK) do
                    DAT[depks(method, :flx, Ed_met, exp)...] = ep_av
                    DAT[depks(method, :err, Ed_met, exp)...] = ep_std
                end
            end

            # inner flxs
            idermap = iJR.load_inners_idermap()
            for exider in DAT[:INNER_FLX_IDERS]

                model_iders = idermap[exider]

                # flxs
                ep_av = ChU.av(model, epout, model_iders[1])
                ep_std = sqrt(ChU.va(model, epout, model_iders[1]))
                if length(model_iders) == 2 # reversible
                    # r = r+ - r-
                    ep_av -= ChU.av(model, epout, model_iders[2])
                    # ep_std += sqrt(ChU.va(model, epout, model_iders[2]))
                    ep_std = NaN
                end

                lock(WLOCK) do
                    DAT[depks(method, :flx, exider, exp)...] = ep_av
                    DAT[depks(method, :err, exider, exp)...] = ep_std
                end
            end

        end # for (exp, method)
    end # for thid

    return DAT
end


text:
module PDKernels

using Reexport

include("./kernels/Kernels.jl")
include("./interpolation/Interpolation.jl")
include("./parameters/Parameters.jl")

@reexport using .Kernels
@reexport using .Interpolation
@reexport using .Parameters

end # module


text:
using DenseSynchrony
using Test

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


text:
# NOTE: Should always be accessed through index(::CVI)
struct VariableIndex
    _value::Int
end
const CVI = VariableIndex
index(cvi::CVI) = cvi._value
Base.convert(::Type{CVI}, vi::VI) = CVI(vi.value)

struct ConstraintIndex{S<:_C_SETS}
    _value::Int
end
const CCI = ConstraintIndex
index(cci::CCI) = cci._value
function Base.convert(::Type{CCI{S}}, ci::CI{SAF,S}) where {S<:_C_SETS}
    return CCI{S}(ci.value)
end

struct ScalarAffineFunction
    coeffs::Vector{Float64}
    indices::Vector{CVI}
    constant::Float64

    function ScalarAffineFunction(
        coeffs::Vector{Float64},
        indices::Vector{CVI},
        constant::Float64,
    )
        @assert length(coeffs) == length(indices)
        return new(coeffs, indices, constant)
    end
end
const CSAF = ScalarAffineFunction
CSAF() = CSAF(Float64[], CVI[], 0.0)

# TODO: Unit test
# NOTE: WE can convert this way (SAF --> CSAF). This is useful for the MOI
# wrapper. However, we are NOT allowed to go the other way: you must use the
# `instantiate` function to make sure that the MOI indices are proper.
function Base.convert(::Type{CSAF}, saf::SAF)
    return CSAF(
        [term.coefficient for term in saf.terms],
        [convert(CVI, term.variable_index) for term in saf.terms],
        saf.constant,
    )
end
function Base.convert(::Type{CSAF}, sv::SV)
    return CSAF([1.0], [convert(CVI, sv.variable)], 0.0)
end

"""
Internal representation for an affine constraint.
"""
struct AffineConstraint{S<:_C_SETS}
    f::CSAF
    s::S
end

# TODO: Unit test
function AffineConstraint(_f::SAF, _s::S) where {S<:_C_SETS}
    f, s = MOIU.normalize_constant(_f, _s)
    @assert f.constant == 0
    return AffineConstraint{S}(convert(CSAF, f), s)
end

# TODO: Unit test
"""
Compute the largest variable index in an affine expression (or constraint).
"""
_max_var_index(ac::AffineConstraint) = _max_var_index(ac.f)
function _max_var_index(csaf::CSAF)
    Base.isempty(csaf.coeffs) && return 0
    return maximum(index(cvi) for cvi in csaf.indices)
end

"""
Internal representation for a polyhedron. We disambiguate different signs in
the linear constraints for type stability reasons (in particular, for working
with bases). Variable bounds are also separately stored as intervals.

It is presumed that `length(bounds)` is equal to the number of variables (i.e.
the ambient dimension), and so the constraints may not contain an index greater
than this value.
"""
mutable struct Polyhedron
    lt_constrs::Vector{AffineConstraint{LT}}
    gt_constrs::Vector{AffineConstraint{GT}}
    et_constrs::Vector{AffineConstraint{ET}}
    bounds::Vector{IN}
    function Polyhedron(aff_constrs::Vector, bounds::Vector{IN})
        n = length(bounds)
        p = new([], [], [], bounds)
        for aff_constr in aff_constrs
            @assert _max_var_index(aff_constr) <= n
            add_constraint(p, aff_constr)
        end
        return p
    end
end

function Polyhedron()
    return Polyhedron(AffineConstraint[], IN[])
end

"Query the ambient dimension (i.e. number of variables)."
ambient_dim(p::Polyhedron) = length(p.bounds)

"""
Increase the ambient dimension of a polyhedron by one. The new variable does
not have bounds, nor does it appear in any of the constraints.
"""
function add_variable(p::Polyhedron)
    push!(p.bounds, IN(-Inf, Inf))
    return nothing
end

# TODO: Unit test
"""
Add a single linear constraint to a polyhedron. The maximum variable index in
the constraint must be less than the ambient dimension of the polyhedron.
"""
function add_constraint(
    p::Polyhedron,
    aff_constr::AffineConstraint{S},
) where {S<:_C_SETS}
    n = ambient_dim(p)
    @assert _max_var_index(aff_constr) <= n
    if S == LT
        push!(p.lt_constrs, aff_constr)
    elseif S == GT
        push!(p.gt_constrs, aff_constr)
    else
        @assert S == ET
        push!(p.et_constrs, aff_constr)
    end
    return nothing
end

"""
Query a single linear constraint from a polyhedron. Returns the `i`-th affine
constraint with "sense" set `T` added to `p`.
"""
get_constraint(p::Polyhedron, T::Type{LT}, i::Int) = p.lt_constrs[i]
get_constraint(p::Polyhedron, T::Type{GT}, i::Int) = p.gt_constrs[i]
get_constraint(p::Polyhedron, T::Type{ET}, i::Int) = p.et_constrs[i]
get_constraints(p::Polyhedron, T::Type{LT}) = p.lt_constrs
get_constraints(p::Polyhedron, T::Type{GT}) = p.gt_constrs
get_constraints(p::Polyhedron, T::Type{ET}) = p.et_constrs

"""
Reports the total number of linear constraints describing a polyhedron. Does
not count any bounds on the variables.
"""
function num_constraints(p::Polyhedron)
    return num_constraints(p, LT) +
           num_constraints(p, GT) +
           num_constraints(p, ET)
end
"""
Reports the number of linear constraints with "sense" set `T` describing a
polyhedron.
"""
num_constraints(p::Polyhedron, T::Type{LT}) = length(p.lt_constrs)
num_constraints(p::Polyhedron, T::Type{GT}) = length(p.gt_constrs)
num_constraints(p::Polyhedron, T::Type{ET}) = length(p.et_constrs)

# TODO: Unit test
"""
Returns `true` if the ambient dimension is 0 and there are no constraints,
and false otherwise.
"""
function Base.isempty(p::Polyhedron)
    return ambient_dim(p) == 0 && num_constraints(p) == 0 && isempty(p.bounds)
end

struct Disjunction
    f::Vector{CSAF}
    s::DisjunctiveConstraints.DisjunctiveSet
end

struct DisjunctiveFormulater{
    T<:DisjunctiveConstraints.AbstractDisjunctiveFormulation,
}
    disjunction::Disjunction
    method::T
end

mutable struct DMIPFormulation
    _feasible_region::Polyhedron
    disjunction_formulaters::Dict{DisjunctiveFormulater,Vector{CVI}}
    _variable_kind::Vector{_V_INT_SETS}
    obj::CSAF

    function DMIPFormulation(
        _feasible_region::Polyhedron,
        _variable_kind::Vector,
        obj::CSAF,
    )
        n = ambient_dim(_feasible_region)
        @assert length(_variable_kind) == n
        @assert _max_var_index(obj) <= n
        return new(_feasible_region, Dict(), _variable_kind, obj)
    end
end

function DMIPFormulation()
    return DMIPFormulation(Polyhedron(), _V_INT_SETS[], CSAF())
end

function all_variables(fm::DMIPFormulation)
    return [CVI(i) for i in 1:length(fm._variable_kind)]
end
num_variables(fm::DMIPFormulation) = ambient_dim(fm._feasible_region)

function add_variable(fm::DMIPFormulation, kind::_V_INT_SETS = nothing)
    add_variable(fm._feasible_region)
    push!(fm._variable_kind, kind)
    return nothing
end

function get_variable_kind(form::DMIPFormulation, cvi::CVI)
    return form._variable_kind[index(cvi)]
end
function set_variable_kind!(form::DMIPFormulation, cvi::CVI, kind::_V_INT_SETS)
    form._variable_kind[index(cvi)] = kind
    return nothing
end

num_constraints(form::DMIPFormulation) = num_constraints(form._feasible_region)
function num_constraints(
    form::DMIPFormulation,
    ::Type{CCI{S}},
) where {S<:_C_SETS}
    return num_constraints(form._feasible_region, S)
end

function num_disjunctive_constraints(form::DMIPFormulation)
    return length(form.disjunction_formulaters)
end

function get_constraint(form::DMIPFormulation, cci::CCI{S}) where {S<:_C_SETS}
    return get_constraint(form._feasible_region, S, index(cci))
end
function get_constraints(
    form::DMIPFormulation,
    ::Type{CCI{S}},
) where {S<:_C_SETS}
    return get_constraints(form._feasible_region, S)
end

function add_constraint(fm::DMIPFormulation, aff_constr::AffineConstraint)
    add_constraint(fm._feasible_region, aff_constr)
    return nothing
end

function Base.isempty(form::DMIPFormulation)
    return isempty(form._feasible_region) &&
           isempty(form.disjunction_formulaters) &&
           isempty(form._variable_kind) &&
           isempty(form.obj.indices) &&
           form.obj.constant == 0
end

function get_bounds(form::DMIPFormulation, cvi::CVI)
    bound = form._feasible_region.bounds[index(cvi)]
    l, u = bound.lower, bound.upper
    if get_variable_kind(form, cvi) isa ZO
        l = max(0, l)
        u = min(1, u)
    end
    return l, u
end

function set_bounds!(form::DMIPFormulation, cvi::CVI, bounds::IN)
    form._feasible_region.bounds[index(cvi)] = bounds
    return nothing
end

function attach_formulater!(
    form::DMIPFormulation,
    formulater::DisjunctiveFormulater,
)
    if haskey(form.disjunction_formulaters, formulater)
        throw(ArgumentError("Formulater cannot be attached twice to a model."))
    end
    start_index = num_variables(form) + 1
    for var_kind in new_variables_to_attach(formulater)
        add_variable(form, var_kind)
    end
    cvis = [CVI(i) for i in start_index:num_variables(form)]
    form.disjunction_formulaters[formulater] = cvis
    return nothing
end


text:
"""
A package for fitting data from auditory evidence accumulation task (Poisson clicks task)
to evidence accumulation model.
"""
__precompile__()

module GeneralUtils

# 3rd party
# using PyCall
# using PyPlot
# using ForwardDiff

export

# non-GUI utilities
    append_to_file,
    replacer,
    print_vector,
    print_vector_g,
    two_level_copy,
    next_file,
    num2fixed_string,
    ascii_key_ize,
    symbol_key_ize,
    vectorize_dict,

# GUI utilities
    safe_axes,
    axisWidthChange,
    axisHeightChange,
    axisMove,
    remove_xtick_labels,
    remove_ytick_labels,
    install_nearest_point_callback,
    install_callback_reporter,
    remove_BP,
    remove_all_BPs,
    get_current_fig_position,
    set_current_fig_position,
    capture_current_figure_configuration,

# Gradient utilities
    FDversion,
    make_dict,
    ForwardDiffZeros,
    get_eltype,
    get_value,
    vgh,
    keyword_vgh

include("non_gui_utils.jl")
# include("gui_utils.jl")
include("gradient_utils.jl")

end # module


text:
@testset "MagneticFieldMeasurement" begin
  # Generate data
  description_ = "Just some testset."

  shp = [3, 3, 3]
  fov = [3.0, 3.0, 3.0]u"mm"
  ctr = [0.0, 0.0, 0.0]u"mm"
  positions_ = RegularGridPositions(shp, fov, ctr)

  gradient = [1.0u"T/m", 1.5u"T/m", 2.0u"T/m"]
  fields_ = fill(0.0u"T", (length(positions_), 3))
  for (idx, pos) in enumerate(positions_)
    fields_[idx, :] = pos.*gradient
  end

  fieldsError_ = fill(0.1u"mT" .|> u"T", (length(positions_), 3))
  fieldsFrequency_ = fill(0.0u"Hz", length(positions_))

  currents_ = fill(1.0u"A", (length(positions_), 3))
  timestamp_ = [DateTime("2021-08-25T07:08:21.881"), DateTime("2021-08-25T07:08:23.881")]
  sensorOffset_ = [0.1u"mm", 0.2u"mm", 0.3u"mm"] .|> u"m"
  temperature_ = fill(20.0u"°C", length(positions_))

  # Set data
  measurement = MagneticFieldMeasurement()

  # Test setters
  description(measurement, description_)
  positions(measurement, positions_)
  fields(measurement, fields_)
  fieldsError(measurement, fieldsError_)
  fieldsFrequency(measurement, fieldsFrequency_)
  currents(measurement, currents_)
  timestamp(measurement, timestamp_)
  sensorOffset(measurement, sensorOffset_)
  temperature(measurement, temperature_)

  @test measurement.description == description_
  @test measurement.positions == positions_
  @test measurement.fields == fields_
  @test measurement.fieldsError == fieldsError_
  @test measurement.fieldsFrequency == fieldsFrequency_
  @test measurement.currents == currents_
  @test measurement.timestamp == timestamp_
  @test measurement.sensorOffset == sensorOffset_
  @test measurement.temperature == temperature_

  # Test getters
  @test description(measurement) == description_
  @test positions(measurement) == positions_
  @test fields(measurement) == fields_
  @test fieldsError(measurement) == fieldsError_
  @test fieldsFrequency(measurement) == fieldsFrequency_
  @test currents(measurement) == currents_
  @test timestamp(measurement) == timestamp_
  @test sensorOffset(measurement) == sensorOffset_
  @test temperature(measurement) == temperature_

  # Save measurement
  magneticMeasurementBaseDir = joinpath(tmpdir, "magnetic")
  mkpath(magneticMeasurementBaseDir)
  magneticMeasurementPath = joinpath(magneticMeasurementBaseDir, "measurement.h5")
  saveMagneticFieldAsHDF5(measurement, magneticMeasurementPath)
  @test isfile(magneticMeasurementPath)

  # Read measurement back
  measurementRead = MagneticFieldMeasurement(magneticMeasurementPath)

  # Test read data
  @test description(measurementRead) == description_
  @test positions(measurementRead) == positions_
  @test fields(measurementRead) == fields_
  @test fieldsError(measurementRead) == fieldsError_
  @test fieldsFrequency(measurementRead) == fieldsFrequency_
  @test currents(measurementRead) == currents_
  @test timestamp(measurementRead) == timestamp_
  @test sensorOffset(measurementRead) == sensorOffset_
  @test temperature(measurementRead) == temperature_
end


text:
import RoME
#include(joinpath(pkgdir(RoME), "test", "runtests.jl")) # February 2021: some tests fail here
#include(joinpath(pkgdir(RoMEPlotting), "test", "runtests.jl"))


text:
module Clustering
  export UPGMA, WPGMA, NeighbourJoining

  import ProfileAligner.Profile

  #typealias DistanceHash{T} Dict{Profile{T}, Dict{Profile{T}, T}}

  type DistanceMatrix{T}
    keys :: Dict{Profile{T}, Int64}
    hsh :: Array{T, 2}
    K :: Vector{T}
    size :: Int
    DistanceMatrix(n :: Int) = new(Dict{Profile{T}, Int64}(), zero(Array(T, 2*n, 2*n)), zero(Array(T, 2*n)), 0)
  end

  function haskey{T}(matrix :: DistanceMatrix{T}, p :: Profile{T})
    Base.haskey(matrix.keys, p)
  end

  function setValue{T}(
      distancesInfo :: DistanceMatrix{T},
      vi :: Profile{T},
      vj :: Profile{T},
      value :: T
      )
    if ! Base.haskey(distancesInfo.keys, vi)
      distancesInfo.keys[vi] = distancesInfo.size + 1
      distancesInfo.size += 1
    end
    if ! Base.haskey(distancesInfo.keys, vj)
      distancesInfo.keys[vj] = distancesInfo.size + 1
      distancesInfo.size += 1
    end
    distancesInfo.hsh[distancesInfo.keys[vi], distancesInfo.keys[vj]] = value
    distancesInfo.hsh[distancesInfo.keys[vj], distancesInfo.keys[vi]] = value
  end

  function getUPGMAWeight{T}(distancesInfo :: DistanceMatrix{T},
                          k :: Profile{T}, f :: Profile{T}, g :: Profile{T})
    (f.numberofstrings*get(distancesInfo, f, k) + g.numberofstrings*get(distancesInfo, g, k))/(g.numberofstrings + f.numberofstrings)
  end

  getCurrentPGMADistance{T}(
          distancesInfo :: DistanceMatrix{T},
          vertices :: Vector{Profile{T}},
          vertsize :: Int64,
          i :: Int,
          j :: Int
          ) = get(distancesInfo, vertices[i], vertices[j])

  function getCurrentNJDistance{T}(
          distancesInfo :: DistanceMatrix{T},
          vertices :: Vector{Profile{T}},
          vertsize :: Int64,
          i :: Int,
          j :: Int
          )
    (vertsize - 2)*get(distancesInfo, vertices[i], vertices[j]) - distancesInfo.K[distancesInfo.keys[vertices[i]]] - distancesInfo.K[distancesInfo.keys[vertices[j]]]
  end

  getWPGMAWeight{T}(
          distancesInfo :: DistanceMatrix{T},
          k :: Profile{T},
          f :: Profile{T},
          g :: Profile{T}
          ) = (get(distancesInfo, f, k) + get(distancesInfo, g, k)) * 0.5

  UPGMA{T}(startvertices :: Vector{Profile{T}}, scoreFunc, mergeFunc) =
          buildTree(startvertices,
              scoreFunc, mergeFunc, getUPGMAWeight)

  WPGMA{T}(startvertices :: Vector{Profile{T}}, scoreFunc, mergeFunc) = buildTree(
          startvertices, scoreFunc, mergeFunc, getWPGMAWeight)

  NeighbourJoining{T}(startvertices :: Vector{Profile{T}}, scoreFunc, mergeFunc) =
    buildTree(startvertices, scoreFunc, mergeFunc,
        getNJWeight, getNJParentDist,
        getCurrentNJDistance)

  function getNJWeight{T}(distancesInfo :: DistanceMatrix{T},
          k :: Profile{T},
          f :: Profile{T},
          g :: Profile{T})
    (
      get(distancesInfo, k, f) +
      get(distancesInfo, k, g) -
      get(distancesInfo, f, g)
    ) * 0.5
  end

  function getPGMAParentDist{T}(
          distancesInfo :: DistanceMatrix{T},
          vertices :: Vector{Profile{T}},
          vertsize :: Int64,
          i :: Int,
          j :: Int)
    get(distancesInfo, vertices[i], vertices[j]) * 0.5

  end

  function getNJParentDist{T}(
          distancesInfo :: DistanceMatrix{T},
          vertices :: Vector{Profile{T}},
          vertsize :: Int64,
          i :: Int,
          j :: Int)
    get(distancesInfo, vertices[i], vertices[j]) / 2.0
      + (distancesInfo.K[distancesInfo.keys[vertices[i]]]
      - distancesInfo.K[distancesInfo.keys[vertices[j]]])/(2.0*(vertsize - 2))
  end

  function buildTree{T}(startvertices :: Vector{Profile{T}},
          scoreFunc, mergeFunc, weightFunc,
          parentWeightFunc = getPGMAParentDist,
          currentDistFunc = getCurrentPGMADistance
          )
    distancesInfo = initDistanceMatrix(startvertices, scoreFunc)
    while(length(startvertices) > 1)
      simplifyTree!(
          distancesInfo,
          startvertices,
          mergeFunc,
          parentWeightFunc,
          weightFunc,
          currentDistFunc
          )
    end
    startvertices[1]
  end

  function simplifyTree!{T}(
          distancesInfo :: DistanceMatrix{T},
          vertices :: Vector{Profile{T}},
          mergeFunc,
          parentWeightFunc,
          weightFunc,
          currentDistFunc
          )
    n = length(vertices)
    #println("in simplify")
    #1. calculate Q

    (mini, minj) = (1, 2)
    Qmin = currentDistFunc(distancesInfo, vertices, n, 1, 2)

    for i = 1 : n - 1
      for j = i + 1 : n
        #println(vertices[i], vertices[j])
        Qcurrent = currentDistFunc(distancesInfo, vertices, n, i, j)
        if Qmin > Qcurrent
          (mini, minj) = (i, j)
          Qmin = Qcurrent
        end
      end
    end
    #
    newT = mergeFunc(vertices[mini], vertices[minj])
    computeDistances!(distancesInfo, vertices, n, mini, minj, newT, parentWeightFunc, weightFunc)

    #
    deleteat!(vertices, [mini, minj])
    push!(vertices, newT)
  end

  function computeDistances!{T}(
          distancesInfo :: DistanceMatrix{T},
          vertices :: Vector{Profile{T}},
          vertsize :: Int64,
          mini :: Int,
          minj :: Int,
          newT :: Profile{T},
          parentWeightFunc,
          weightFunc
          )
    #compute new distances from new to previous 2
    #distancesInfo[newT] = Dict{Profile{T}, T}()

    f = vertices[mini]
    g = vertices[minj]
    t2 = one(T) + one(T)
    setValue(distancesInfo, newT, f, parentWeightFunc(distancesInfo, vertices, vertsize, mini, minj))
    setValue(distancesInfo, newT, g, get(distancesInfo, f, g) - get(distancesInfo, newT, f))
    setValue(distancesInfo, newT, newT, zero(T)) #assume that profile aligns to profile with zero distance
    #compute to all others
    for i in 1 : length(vertices)
      i == mini && continue
      i == minj && continue
      k = vertices[i]
      setValue(distancesInfo, newT, k, weightFunc(distancesInfo, k, f, g))
      distancesInfo.K[distancesInfo.keys[k]] += (-get(distancesInfo, k, g) - get(distancesInfo, k, f) + get(distancesInfo, k, newT))
      distancesInfo.K[distancesInfo.keys[newT]] += get(distancesInfo, k, newT)
    end
  end

  function initDistanceMatrix{T}(startvertices :: Vector{Profile{T}}, scoreFunc)
    n = length(startvertices)
    distancesInfo = DistanceMatrix{T}(n)
    for i in 1 : n
      vi = startvertices[i]
      for j in i : n
        vj = startvertices[j]
        setValue(distancesInfo, vi,vj,
                scoreFunc(vi, vj))
        distancesInfo.K[distancesInfo.keys[vi]] += get(distancesInfo, vi, vj)
        if i != j
          distancesInfo.K[distancesInfo.keys[vj]] += get(distancesInfo, vi, vj)
        end
      end
    end
    distancesInfo
  end

  function get2{T}(distancesInfo :: DistanceMatrix{T}, k1 :: Profile{T}, k2 :: Profile{T})
    distancesInfo.hsh[distancesInfo.keys[k1], distancesInfo.keys[k2]]
  end

  function getKValue{T}(distancesInfo :: DistanceMatrix{T}, k1 :: Profile{T})
    distancesInfo.K[distancesInfo.keys[k1]]
  end

  function get{T}(distancesInfo :: DistanceMatrix{T}, k1 :: Profile{T}, k2 :: Profile{T})
    Base.haskey(distancesInfo.keys, k1) && Base.haskey(distancesInfo.keys, k2) && return get2(distancesInfo, k1, k2)
    Base.haskey(distancesInfo.keys, k1) && error("invalid k2 key in distance hash")
    Base.haskey(distancesInfo.keys, k2) && error("invalid k1 key in distance hash")
    error("invalid keys in distance hash")
  end
end


text:
__precompile__()

module OpenSpecFun

deps = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(deps)
    include(deps)
else
    error("OpenSpecFun.jl not properly installed. Please run Pkg.build(\"OpenSpecFun\").")
end

import Base.Math: @horner, nan_dom_err

# export
#     # bessel
#     airy,
#     airyai,
#     airyaiprime,
#     airybi,
#     airybiprime,
#     airyprime,
#     airyx,
#     besselh,
#     besselhx,
#     besseli,
#     besselix,
#     besselj,
#     besselj0,
#     besselj1,
#     besseljx,
#     besselk,
#     besselkx,
#     bessely,
#     bessely0,
#     bessely1,
#     besselyx,
#     hankelh1,
#     hankelh1x,
#     hankelh2,
#     hankelh2x,
#     # erf
#     erf,
#     erfc,
#     erfcinv,
#     erfcx,
#     erfi,
#     erfinv

const libm = Base.libm_name

for f in ["bessel", "erf", "gamma"]
    include("$f.jl")
end

end # module


text:
using Base: String
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Interpreter exceptions

"""
    InterpreterException <: Exception

Abstract type for all SceneLang interpreter errors.

See also:
[`BadCharacter`](@ref),
[`UnfinishedExpression`](@ref),
[`UndefinedIdentifier`](@ref),
[`WrongTokenType`](@ref),
[`WrongValueType`](@ref),
[`InvalidKeyword`](@ref),
[`InvalidType`](@ref),
[`InvalidCommand`](@ref),
[`InvalidExpression`](@ref),
[`InvalidSymbol`](@ref),
[`InvalidNumber`](@ref),
[`InvalidSize`](@ref),
[`InvalidFilePath`](@ref),
[`IdentifierRedefinition`](@ref)
[`SettingRedefinition`](@ref)
[`UndefinedSetting`](@ref)
"""
abstract type InterpreterException <: Exception end

macro make_exception(name::Symbol, descriptive_str::AbstractString)
    struct_doc =
        """
            $name <: InterpreterException

        $descriptive_str

        See also: [`InterpreterException`](@ref)

        # Fields

        - `location::SourceLocation`: location of the error
        - `msg::AbstractString`: descriptive error message
        - `len::Int`: how many characters are involved in the error
        """
    constructor_doc =
        """
            $name(location::SourceLocation, msg::AbstractString)

        Construct an instance of [`$name`](@ref) with `len = 1`.
        """

    quote
        @doc $struct_doc $name

        struct $(esc(name)) <: InterpreterException
            location::SourceLocation
            msg::AbstractString
            len::Int
        end

        @doc $constructor_doc $name(::SourceLocation, ::AbstractString)

        function $(esc(name))(location::SourceLocation, msg::AbstractString)
            $(esc(name))(location, msg, 1)
        end
    end
end

function Base.showerror(io::IO, e::InterpreterException)
    print(io, typeof(e))
    printstyled(io, " @ ", e.location, color=:light_black)
    println(io)
    printstyled(io, e.msg, color=:red)
    iszero(e.location.line_num) &&
        return
    println(io)
    printstyled(io, "source: ", color=:light_black)
    println(io, read_at_line(e.location.file_name, e.location.line_num))
    printstyled(io, " " ^ (e.location.col_num + 7), color=:light_black)
    printstyled(io, "^" ^ e.len, color=:red)
end

function Base.showerror(io::IO, e::InterpreterException, bt; backtrace = false)
    try
        showerror(io, e)
    finally
        nothing
    end
end

@make_exception BadCharacter           "There is an invalid character in the SceneLang script."
@make_exception UnfinishedExpression   "A special environment (e.g. a string, mathematical expression, list...) has been opened and not closed."
@make_exception UndefinedIdentifier    "The given identifier has not been defined in the script."
@make_exception WrongTokenType         "The given value is of a different type than expected by the syntax."
@make_exception WrongValueType         "The given value has a different type than expected by the syntax."
@make_exception InvalidKeyword         "The given keyword is not valid in the given context."
@make_exception InvalidType            "The given type does not exist."
@make_exception InvalidCommand         "The given command does not exist."
@make_exception InvalidExpression      "The given expression contains invalid elements. Capabilities are restrained to contain malicious code injection."
@make_exception InvalidSymbol          "The given symbol is not valid in the given context."
@make_exception InvalidNumber          "The token has an invalid numerical format."
@make_exception InvalidSize            "The given collection has an invalid size in the given context."
@make_exception InvalidFilePath        "The given file path is invalid."
@make_exception IdentifierRedefinition "An identifier is being redefined without being unset first."
@make_exception SettingRedefinition    "A rendering setting is being defined multiple times."
@make_exception UndefinedSetting       "One or more necessary rendering settings have not been set up in the SceneLang script."

function Base.showerror(io::IO, e::UndefinedSetting)
    print(io, typeof(e))
    printstyled(io, " @ ", e.location.file_name, color=:light_black)
    println(io)
    printstyled(io, e.msg, color=:red)
end


text:
import GroupTools.FiniteGroup

import GroupTools.AbstractGroup
import GroupTools.AbstractSymmetry

import GroupTools.element
import GroupTools.elements
import GroupTools.element_name
import GroupTools.element_names

import GroupTools.isidentity
import GroupTools.apply_operation

import GroupTools.group
import GroupTools.group_order
import GroupTools.group_multiplication_table
import GroupTools.generate_multiplication_table

import GroupTools.⋊
import GroupTools.⋉

export ⋊
export ⋉
export group_order
export group_multiplication_table
export generate_multiplication_table
export apply_operation
export isidentity
export element
export elements
export element_names
export element_names

export AbstractGroup
export AbstractSymmetry
export FiniteGroup


text:
module Turing2MonteCarloMeasurements

using MonteCarloMeasurements, Turing, AxisArrays

export Particles, truthplot

function MonteCarloMeasurements.Particles(t::NTuple{N, <:AxisArrays.AxisArray}; crop=0) where N
    [adapted_particles(t[i].data[crop+1:end,:]) for i in 1:N]
end

function MonteCarloMeasurements.Particles(a::AxisArrays.AxisArray; crop=0)
    adapted_particles(a.data[crop+1:end,:])
end

"""
    Particles(chain::Chains; crop=0)

Return a named tuple of particles or vector of particles where the keys are the symbols in the Turing model that produced the chain. `crop>0` can be given to discard a number of samples from the beginning of the chain.

# Example:
In the example below, the model contained variables named `reviewer_pop_bias` and `reviewer_pop_gain`
```
julia> cp = Particles(chain, crop=500);

julia> cp.reviewer_pop_bias
Part1000(0.2605 ± 0.72)

julia> cp.reviewer_pop_gain
Part1000(0.1831 ± 0.62)
```
"""
function MonteCarloMeasurements.Particles(chain::Chains; crop=0)
    p = get_params(chain)
    (;collect((k => Particles(getfield(p,k); crop=crop) for k in keys(p)))...)
end

function adapted_particles(v)
    T = float(typeof(v[1]))
    Particles(vec(T.(v)))
end

function truthplot(truth::AbstractVector, p::AbstractVector, fn; kwargs...)
    figs = truthplot.(truth, p, fn, axis=false, title=string(fn), titlefont=8)
    plot(figs...)
end

function truthplot(truth::Number, p::AbstractParticles, fn; kwargs...)
    fig = density(p; title=string(fn), lab="ess: $(round(Int,MonteCarloMeasurements.ess(p)))", kwargs...)
    vline!(fig,[truth], lab="")
    fig
end

function truthplot(truth::NamedTuple, p::NamedTuple; kwargs...)
    plots = []
    for fn in fieldnames(typeof(truth))
        f   = getfield(truth,fn)
        local pf
        try
            pf  = getfield(p,fn)
        catch
            continue
        end
        fig = truthplot(f, pf, fn; kwargs...)
        push!(plots,fig)
    end
    display.(plots)
    plots
end

end # module


text:
@tm1def mutable struct Dimension
  Name::Union{String,Nothing}
  UniqueName::Union{String,Nothing}
  Attributes::Union{Dict,Nothing}
  AllLeavesHierarchyName::Union{String,Nothing}
  Hierarchies::Union{Vector{Hierarchy},Nothing}
end

Dimension(name::AbstractString) = Dimension(Dict("Name" => name))

namefield(dimension::Dimension) = dimension.Name

@api_default function get_all_dimensions(api::TM1API; options...)
  params = Dict("\$expand" => "Hierarchies(\$expand=*)")
  result, page_data = tm1_get_paged_json(api, "Dimensions"; params = params, options...)
  map(Dimension, get(result, "value", [])), page_data
end

@api_default function get_all_model_dimensions(api::TM1API; options...)
  params = Dict("\$expand" => "Hierarchies(\$expand=*)")
  result, page_data =
    tm1_get_paged_json(api, "ModelDimensions()"; params = params, options...)
  map(Dimension, get(result, "value", [])), page_data
end

@api_default function get_all_control_dimensions(api::TM1API; options...)
  params = Dict("\$expand" => "Hierarchies(\$expand=*)")
  result, page_data =
    tm1_get_paged_json(api, "ControlDimensions()"; params = params, options...)
  map(Dimension, get(result, "value", [])), page_data
end

@api_default function get_dimension(api::TM1API, dimension_name::AbstractString; options...)
  # I'm not entirely sure if this is behaving as desired
  params = Dict("\$expand" => "Hierarchies(\$expand=*)")
  result =
    tm1_get_json(api, "Dimensions('" * dimension_name * "')"; params = params, options...)
  Dimension(result)
end

@api_default function delete_dimension(
  api::TM1API,
  dimension_name::AbstractString;
  options...,
)
  # not sure what to return here
  tm1_delete(api, "Dimensions('" * dimension_name * "')"; options...)
end


text:
# TODO: implement writing a dict of graphs

function _graphml_read_one_graph(e::XMLElement, isdirected::Bool)
    nodes = Dict{AbstractString,Int}()
    edges = Vector{Edge}()

    nodeid = 1
    for f in child_elements(e)
        if name(f) == "node"
            nodes[attribute(f, "id")] = nodeid
            nodeid += 1
        elseif name(f) == "edge"
            n1 = attribute(f, "source")
            n2 = attribute(f, "target")
            push!(edges, Edge(nodes[n1], nodes[n2]))
        else
            warn("Skipping unknown node '$(name(f))'")
        end
    end
    #Put data in graph
    g = (isdirected ? DiGraph : Graph)(length(nodes))
    for edge in edges
        add_edge!(g, edge)
    end
    return g
end

function loadgraphml(io::IO, gname::AbstractString)
    xdoc = parse_string(readall(io))
    xroot = root(xdoc)  # an instance of XMLElement
    name(xroot) == "graphml" || error("Not a GraphML file")

    # traverse all its child nodes and print element names
    for c in child_nodes(xroot)  # c is an instance of XMLNode
        if is_elementnode(c)
            e = XMLElement(c)  # this makes an XMLElement instance
            if name(e) == "graph"
                edgedefault = attribute(e, "edgedefault")
                isdir = edgedefault=="directed" ? true :
                             edgedefault=="undirected" ? false : error("Unknown value of edgedefault: $edgedefault")
                if has_attribute(e, "id")
                    graphname = attribute(e, "id")
                else
                    graphname =  isdir ? "digraph" : "graph"
                end
                gname == graphname && return _graphml_read_one_graph(e, isdir)
            else
                warn("Skipping unknown XML element '$(name(e))'")
            end
        end
    end
    error("Graph $gname not found")
end

function loadgraphml_mult(io::IO)
    xdoc = parse_string(readall(io))
    xroot = root(xdoc)  # an instance of XMLElement
    name(xroot) == "graphml" || error("Not a GraphML file")

    # traverse all its child nodes and print element names
    graphs = Dict{AbstractString, SimpleGraph}()
    for c in child_nodes(xroot)  # c is an instance of XMLNode
        if is_elementnode(c)
            e = XMLElement(c)  # this makes an XMLElement instance
            if name(e) == "graph"
                edgedefault = attribute(e, "edgedefault")
                isdir = edgedefault=="directed" ? true :
                             edgedefault=="undirected" ? false : error("Unknown value of edgedefault: $edgedefault")
                if has_attribute(e, "id")
                    graphname = attribute(e, "id")
                else
                    graphname =  isdir ? "digraph" : "graph"
                end
                graphs[graphname] =  _graphml_read_one_graph(e, isdir)
            else
                warn("Skipping unknown XML element '$(name(e))'")
            end
        end
    end
    return graphs
end

function savegraphml_mult(io::IO, graphs::Dict)
    xdoc = XMLDocument()
    xroot = create_root(xdoc, "graphml")
    set_attribute(xroot,"xmlns","http://graphml.graphdrawing.org/xmlns")
    set_attribute(xroot,"xmlns:xsi","http://www.w3.org/2001/XMLSchema-instance")
    set_attribute(xroot,"xsi:schemaLocation","http://graphml.graphdrawing.org/xmlns/1.0/graphml.xsd")

    for (gname, g) in graphs
        xg = new_child(xroot, "graph")
        set_attribute(xg,"id",gname)
        strdir = is_directed(g) ? "directed" : "undirected"
        set_attribute(xg,"edgedefault",strdir)

        for i=1:nv(g)
            xv = new_child(xg, "node")
            set_attribute(xv,"id","n$(i-1)")
        end

        m = 0
        for e in edges(g)
            xe = new_child(xg, "edge")
            set_attribute(xe,"id","e$m")
            set_attribute(xe,"source","n$(src(e)-1)")
            set_attribute(xe,"target","n$(dst(e)-1)")
            m += 1
        end
    end
    show(io, xdoc)
    return length(graphs)
end

savegraphml(io::IO, g::SimpleGraph, gname::AbstractString) =
    savegraphml_mult(io, Dict(gname=>g))

filemap[:graphml] = (loadgraphml, loadgraphml_mult, savegraphml, savegraphml_mult)


text:
#EasyConf: Interface to read/write simple configuration files

module EasyConf

export EasyConfFile
export getvalue

# Types
################################################################################

mutable struct ConfigDict
	d::Dict{String, Any}
end
ConfigDict() = ConfigDict(Dict{String, Any}())

mutable struct EasyConfFile
	s::IOStream
end

# Read/Write functions
################################################################################

#-------------------------------------------------------------------------------
function Base.open(::Type{EasyConfFile}, path::String, args...)
	s = open(path, args...)
	return EasyConfFile(s)
end

#-------------------------------------------------------------------------------
Base.close(f::EasyConfFile) = close(f.s)

#-------------------------------------------------------------------------------
function Base.read(f::EasyConfFile)
	cfg = ConfigDict()

	for line in eachline(f.s)
		elements = split(line, "=")
		if length(elements) != 2
			continue
		end
		key = String(strip(elements[1])) #Need string (not SubString)
		cfg[key] = strip(elements[2])
	end

	return cfg
end

#-------------------------------------------------------------------------------
function Base.write(f::EasyConfFile, cfg::ConfigDict)
	for (k,v) in cfg.d
		write(f.s, "$k = $v")
	end
end

#-------------------------------------------------------------------------------
function Base.read(::Type{EasyConfFile}, path::String)
	local f, cfg
	try
		f = open(EasyConfFile, path, "r")
	catch
		#Cannot open; Return empty dictionary
		return ConfigDict()
	end

	try
		cfg = read(f)
	finally
		close(f)
	end
	return cfg
end

#-------------------------------------------------------------------------------
function Base.write(::Type{EasyConfFile}, path::String, cfg::ConfigDict)
	f = open(EasyConfFile, path, "w")
	try
		write(f, cfg)
	finally
		close(f)
	end
end

# Other tools
################################################################################

Base.setindex!(c::ConfigDict, v, k::String) = (c.d[k] = v)
Base.getindex(c::ConfigDict, k::String) = c.d[k]

# Safe means of reading from ConfigDict
#-------------------------------------------------------------------------------
function getvalue(c::ConfigDict, k::String, vdefault)
	local result
	try
		result = c[k]
	catch
		return vdefault
	end
	return result
end

end #EasyConf

max_stars_repo_name: ma-laforge/AudioStick

max_issues_repo_name: ma-laforge/AudioStick

max_forks_repo_name: ma-laforge/AudioStick


text:
include("../Resistances/heat.jl")
include("../../Utilities/summary.jl")

using Measurements, Modia
#ENV["MODIA_PLOT"] = "PyPlot"
Modia.usePlotPackage("PyPlot")
@usingModiaPlot

# Boundary model(s)
T_a = 293.15u"K"
T_b = 313.15u"K"

r_inner = measurement(1.0,0.1)u"m"
r_outer = 2.0u"m"
lambda = 5.0u"W/(m*K)"

# Base Classes
PartialTest = Model(
    Q_total=Var(init=0.0u"J"),
    partialEquations =:[
        der(Q_total) = Qb_flow
    ]
)

# Sphere Model
sphere_model = PartialTest | Model(
    boundaryT_a = Temperature | Map(T=T_a),
    boundaryT_b = Temperature | Map(T=T_b),
    resistor = Sphere | Map(r_inner=r_inner,r_outer=r_outer,lambda=lambda),
    equations =:[
        Qb_flow = resistor.port_a.Q_flow
    ],
    connect =:[
        (boundaryT_a.port,resistor.port_a)
        (resistor.port_b,boundaryT_b.port)
    ]
)

result = @instantiateModel(sphere_model, log=false,FloatType=Measurement{Float64})
@time simulate!(result, stopTime=1.0, tolerance=1e-4, log=false)
y_min, t_min, y_max, t_max = reportMinMax("Q_total","Sphere Model",result)
errorCheck(-2513.2742u"J"±0.0u"J", y_min)
plot(result, ["Q_total"], figure=1)

# Convection model
surfaceArea = 2u"m^2"
alpha = 15u"W/(m^2*K)"

convection_model = PartialTest | Model(
    boundaryT_a = Temperature | Map(T=T_a),
    boundaryT_b = Temperature | Map(T=T_b),
    resistor = Convection | Map(surfaceArea = surfaceArea, alpha=alpha),
    equations =:[
        Qb_flow = resistor.port_a.Q_flow
    ],
    connect =:[
        (boundaryT_a.port,resistor.port_a)
        (resistor.port_b,boundaryT_b.port)
    ]
)

result = @instantiateModel(convection_model, log=false,FloatType=Measurement{Float64})
@time simulate!(result, stopTime=1.0, tolerance=1e-4, log=false)
y_min, t_min, y_max, t_max = reportMinMax("Q_total","Convection Model",result)
errorCheck(-600u"J"±0.0u"J", y_min)
plot(result, ["Q_total"], figure=1)

# Contact model
surfaceArea = 2u"m^2" # units hid until constant unit issue solved
Rc_pp = 5u"m^2*K/W"

contact_model = PartialTest | Model(
    boundaryT_a = Temperature | Map(T=T_a),
    boundaryT_b = Temperature | Map(T=T_b),
    resistor = Contact | Map(surfaceArea = surfaceArea, Rc_pp=Rc_pp),
    equations =:[
        Qb_flow = resistor.port_a.Q_flow
    ],
    connect =:[
        (boundaryT_a.port,resistor.port_a)
        (resistor.port_b,boundaryT_b.port)
    ]
)

result = @instantiateModel(contact_model, log=false,FloatType=Measurement{Float64})
@time simulate!(result, stopTime=1.0, tolerance=1e-4, log=false)
y_min, t_min, y_max, t_max = reportMinMax("Q_total","Contact Model",result)
errorCheck(-8u"J"±0.0u"J", y_min)
plot(result, ["Q_total"], figure=1)

# Radiation exact model
surfaceArea = 2u"m^2"
epsilon = 0.5u"m/m"

radiationExact_model = PartialTest | Model(
    boundaryT_a = Temperature | Map(T=T_a),
    boundaryT_b = Temperature | Map(T=T_b),
    resistor = Radiation | Map(surfaceArea = surfaceArea, epsilon=epsilon),
    equations =:[
        Qb_flow = resistor.port_a.Q_flow
    ],
    connect =:[
        (boundaryT_a.port,resistor.port_a)
        (resistor.port_b,boundaryT_b.port)
    ]
)

result = @instantiateModel(radiationExact_model, log=false,FloatType=Measurement{Float64})
@time simulate!(result, stopTime=1.0, tolerance=1e-4, log=false)
y_min, t_min, y_max, t_max = reportMinMax("Q_total","Radiation Exact Model",result)
errorCheck(-126.51638u"J"±0.0u"J", y_min)
plot(result, ["Q_total"], figure=1)

# Radiation approximate model
radiationApproximate_model = PartialTest | Model(
    boundaryT_a = Temperature | Map(T=T_a),
    boundaryT_b = Temperature | Map(T=T_b),
    resistor = Radiation | Map(surfaceArea = surfaceArea, epsilon=epsilon, useExact=false),
    equations =:[
        Qb_flow = resistor.port_a.Q_flow
    ],
    connect =:[
        (boundaryT_a.port,resistor.port_a)
        (resistor.port_b,boundaryT_b.port)
    ]
)

result = @instantiateModel(radiationApproximate_model, log=false,FloatType=Measurement{Float64})
@time simulate!(result, stopTime=1.0, tolerance=1e-4, log=false)
y_min, t_min, y_max, t_max = reportMinMax("Q_total","Radiation Approximate Model",result)
errorCheck(-126.51638u"J"±0.0u"J", y_min)
plot(result, ["Q_total"], figure=1)


text:
# adapted from papiStdEventDefs.h
@enum(Preset::Cuint,
    L1_DCM = PAPI_PRESET_MASK,  # Level 1 data cache misses
    L1_ICM,  # Level 1 instruction cache misses
    L2_DCM,  # Level 2 data cache misses
    L2_ICM,  # Level 2 instruction cache misses
    L3_DCM,  # Level 3 data cache misses
    L3_ICM,  # Level 3 instruction cache misses
    L1_TCM,  # Level 1 total cache misses
    L2_TCM,  # Level 2 total cache misses
    L3_TCM,  # Level 3 total cache misses
    CA_SNP,  # Snoops
    CA_SHR,  # Request for shared cache line (SMP)
    CA_CLN,  # Request for clean cache line (SMP)
    CA_INV,  # Request for cache line Invalidation (SMP)
    CA_ITV,  # Request for cache line Intervention (SMP)
    L3_LDM,  # Level 3 load misses
    L3_STM,  # Level 3 store misses
    BRU_IDL, # Cycles branch units are idle
    FXU_IDL, # Cycles integer units are idle
    FPU_IDL, # Cycles floating point units are idle
    LSU_IDL, # Cycles load/store units are idle
    TLB_DM,  # Data translation lookaside buffer misses
    TLB_IM,  # Instr translation lookaside buffer misses
    TLB_TL,  # Total translation lookaside buffer misses
    L1_LDM,  # Level 1 load misses
    L1_STM,  # Level 1 store misses
    L2_LDM,  # Level 2 load misses
    L2_STM,  # Level 2 store misses
    BTAC_M,  # BTAC miss
    PRF_DM,  # Prefetch data instruction caused a miss
    L3_DCH,  # Level 3 Data Cache Hit
    TLB_SD,  # Xlation lookaside buffer shootdowns (SMP)
    CSR_FAL, # Failed store conditional instructions
    CSR_SUC, # Successful store conditional instructions
    CSR_TOT, # Total store conditional instructions
    MEM_SCY, # Cycles Stalled Waiting for Memory Access
    MEM_RCY, # Cycles Stalled Waiting for Memory Read
    MEM_WCY, # Cycles Stalled Waiting for Memory Write
    STL_ICY, # Cycles with No Instruction Issue
    FUL_ICY, # Cycles with Maximum Instruction Issue
    STL_CCY, # Cycles with No Instruction Completion
    FUL_CCY, # Cycles with Maximum Instruction Completion
    HW_INT,  # Hardware interrupts
    BR_UCN,  # Unconditional branch instructions executed
    BR_CN,   # Conditional branch instructions executed
    BR_TKN,  # Conditional branch instructions taken
    BR_NTK,  # Conditional branch instructions not taken
    BR_MSP,  # Conditional branch instructions mispred
    BR_PRC,  # Conditional branch instructions corr. pred
    FMA_INS, # FMA instructions completed
    TOT_IIS, # Total instructions issued
    TOT_INS, # Total instructions executed
    INT_INS, # Integer instructions executed executed
    FP_INS,  # Floating point instructions executed
    LD_INS,  # Load instructions executed
    SR_INS,  # Store instructions executed
    BR_INS,  # Total branch instructions executed
    VEC_INS, # Vector/SIMD instructions executed (could include integer)
    RES_STL, # Cycles processor is stalled on resource
    FP_STAL, # Cycles any FP units are stalled
    TOT_CYC, # Total cycles executed
    LST_INS, # Total load/store inst. executed
    SYC_INS, # Sync. inst. executed
    L1_DCH,  # L1 D Cache Hit
    L2_DCH,  # L2 D Cache Hit
    L1_DCA,  # L1 D Cache Access
    L2_DCA,  # L2 D Cache Access
    L3_DCA,  # L3 D Cache Access
    L1_DCR,  # L1 D Cache Read
    L2_DCR,  # L2 D Cache Read
    L3_DCR,  # L3 D Cachee Read
    L1_DCW,  # L1 D Cache Write
    L2_DCW,  # L2 D Cache Write
    L3_DCW,  # L3 D Cache Write
    L1_ICH,  # L1 instruction cache hits
    L2_ICH,  # L2 instruction cache hits
    L3_ICH,  # L3 instruction cache hits
    L1_ICA,  # L1 instruction cache accesses
    L2_ICA,  # L2 instruction cache accesses
    L3_ICA,  # L3 instruction cache accesses
    L1_ICR,  # L1 instruction cache reads
    L2_ICR,  # L2 instruction cache reads
    L3_ICR,  # L3 instruction cache reads
    L1_ICW,  # L1 instruction cache writes
    L2_ICW,  # L2 instruction cache writes
    L3_ICW,  # L3 instruction cache writes
    L1_TCH,  # L1 total cache hits
    L2_TCH,  # L2 total cache hits
    L3_TCH,  # L3 total cache hits
    L1_TCA,  # L1 total cache accesses
    L2_TCA,  # L2 total cache accesses
    L3_TCA,  # L3 total cache accesses
    L1_TCR,  # L1 total cache reads
    L2_TCR,  # L2 total cache reads
    L3_TCR,  # L3 total cache reads
    L1_TCW,  # L1 total cachee writes
    L2_TCW,  # L2 total cache writes
    L3_TCW,  # L3 total cache writes
    FML_INS, # FM ins
    FAD_INS, # FA ins
    FDV_INS, # FD ins
    FSQ_INS, # FSq ins
    FNV_INS, # Finv ins
    FP_OPS,  # Floating point operations executed
    SP_OPS,  # Floating point operations executed; optimized to count scaled single precision vector operations
    DP_OPS,  # Floating point operations executed; optimized to count scaled double precision vector operations
    VEC_SP,  # Single precision vector/SIMD instructions
    VEC_DP,  # Double precision vector/SIMD instructions
    REF_CYC, # Reference clock cycles
    END      # This should always be last!
)

is_preset(evt::Cuint) = ((evt & PAPI_PRESET_MASK) != 0) && ((evt & PAPI_NATIVE_MASK) == 0)
is_native(evt::Cuint) = ((evt & PAPI_PRESET_MASK) == 0) && ((evt & PAPI_NATIVE_MASK) != 0)

primitive type Native sizeof(Cuint)*8 end
Native(x::Integer) = Base.bitcast(Native, convert(Cuint, x))
Base.cconvert(::Type{T}, x::Native) where {T<:Integer} = Base.bitcast(T, x)::T

const Event = Union{Preset, Native}

Base.promote_rule(::Type{Native}, ::Type{Preset}) = Event

function exists(evt::Event)
	errcode = ccall((:PAPI_query_event, libpapi), Cint, (Cuint,), evt)
	if errcode != PAPI_OK && errcode != PAPI_ENOEVNT
		throw(PAPIError(errcode))
	end

    return errcode == PAPI_OK
end

"""
    event_to_name(evt::Event)

Converts an event into a name.
"""
function event_to_name(evt::Event)
	str_buf = Vector{UInt8}(undef,PAPI_MAX_STR_LEN)
	@papichk ccall((:PAPI_event_code_to_name, libpapi), Cint, (Cuint, Ptr{UInt8}), evt, str_buf)
	unsafe_string(pointer(str_buf))
end

"""
    name_to_event(evt::Event)

Converts event name into an event
"""
function name_to_event(name::AbstractString)
    evt = Ref{Cuint}()
    @papichk ccall((:PAPI_event_name_to_code, libpapi), Cint, (Cstring, Ptr{Cuint}), name, evt)
    if is_preset(evt[])
        Preset(evt[])
    else
        Native(evt[])
    end
end

function Base.show(io::IO, evt::Native)
  name = event_to_name(evt)
  print(io, name)
end

"""
    available_presets()

Returns a list of the presets available on the current platform.
"""
available_presets() = filter(exists, instances(Preset))

macro event_str(n)
    quote name_to_event($n) end
end


text:
module Default

export default, Optional

struct default{T}
    value::T
end

function (def::default)(::Nothing)
    return def.value
end

function (def::default)(val::Any)
    return val
end

Optional{T} = Union{Nothing,T}
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    Coordinate transformations related with the geodetic and geocentric
#    coordinates.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [2] ESA Navipedia: http://www.navipedia.net/
#
#   [3] mu-blox ag (1999). Datum Transformations of GPS Positions. Application
#       Note.
#
#   [4] ISO TC 20/SC 14 N (2011). Geomagnetic Reference Models.
#
#   [5] Borkowski, K. M (1987). Transformation of geocentric to geodetic
#       coordinates without approximations. Astrophysics and Space Science, vol.
#       139, pp. 1-4.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export ecef_to_geodetic, geodetic_to_ecef
export geocentric_to_geodetic, geodetic_to_geocentric

"""
    ecef_to_geodetic(r_e::AbstractVector)

Convert the vector `r_e` [m] represented in the Earth-Centered, Earth-Fixed
(ECEF) reference frame into Geodetic coordinates (WGS-84).

!!! info
    The algorithm is based in **[3]**.

# Returns

* Latitude [rad].
* Longitude [rad].
* Altitude [m].

# Reference

- **[3]**: mu-blox ag (1999). Datum Transformations of GPS Positions.
    Application Note.
"""
function ecef_to_geodetic(r_e::AbstractVector)
    # Auxiliary variables.
    X = r_e[1]
    Y = r_e[2]
    Z = r_e[3]

    # Auxiliary variables.
    p = √(X^2 + Y^2)
    θ = atan(Z * a_wgs84, p * b_wgs84)
    e_wgs84² = e_wgs84^2
    sin_θ, cos_θ = sincos(θ)

    # Compute Geodetic.
    lon = atan(Y, X)
    lat = atan(
        Z + el_wgs84^2 * b_wgs84 * sin_θ^3,
        p -   e_wgs84² * a_wgs84 * cos_θ^3
    )

    sin_lat, cos_lat = sincos(lat)

    N = a_wgs84/√(1 - e_wgs84² * sin_lat^2)

    # Avoid singularity if we are near the poles (~ 1 deg according to [1,
    # p.172]). Note that `cosd(1) = -0.01745240643728351`.
    if !(-0.01745240643728351 < cos_lat < 0.01745240643728351)
        h = p/cos_lat - N
    else
        h = Z/sin_lat - N*(1 - e_wgs84²)
    end

    return lat, lon, h
end

"""
    geodetic_to_ecef(lat::Number, lon::Number, h::Number)

Convert the latitude `lat` [rad], longitude `lon` [rad], and altitude `h` \\[m]
(WGS-84) into a vector represented on the Earth-Centered, Earth-Fixed (ECEF)
reference frame.

!!! info
    The algorithm is based in **[3]**.

# Reference

- **[3]**: mu-blox ag (1999). Datum Transformations of GPS Positions.
    Application Note.
"""
function geodetic_to_ecef(lat::Number, lon::Number, h::Number)
    # Auxiliary variables.
    sin_lat, cos_lat = sincos(lat)
    sin_lon, cos_lon = sincos(lon)

    # Radius of curvature [m].
    N = a_wgs84/√(1 - e_wgs84^2 * sin_lat^2)

    # Compute the position in ECEF frame.
    return SVector(
        (                        N + h) * cos_lat * cos_lon,
        (                        N + h) * cos_lat * sin_lon,
        ((b_wgs84 / a_wgs84)^2 * N + h) * sin_lat
    )
end

"""
    geocentric_to_geodetic(ϕ_gc::Number, r::Number)

Compute the geodetic latitude and altitude (WGS-84) from the geocentric latitude
`ϕ_gc` (-π/2, π/2) [rad] and radius `r` [m]. Notice that the longitude is the
same in both geocentric and geodetic coordinates.

!!! info
    The algorithm is based in **[5]**.

# Returns

* Geodetic latitude [rad].
* Altitude above the reference ellipsoid (WGS-84) [m].

# References

- **[5]** Borkowski, K. M (1987). Transformation of geocentric to geodetic
    coordinates without approximations. Astrophysics and Space Science, vol.
    139, pp. 1-4.
"""
function geocentric_to_geodetic(ϕ_gc::Number, r::Number)
    # Obtain the `z` component and the equatorial component `re`.
    sin_ϕ_gc, cos_ϕ_gc = sincos(ϕ_gc)
    re = r * cos_ϕ_gc
    z  = r * sin_ϕ_gc
    sign_z = z == 0 ? +1 : sign(z)

    # Auxiliary variables.
    a  = a_wgs84
    a² = a^2
    b  = sign_z * b_wgs84
    b² = b^2

    # Compute the parameters.
    E  = (b * z - (a² - b²)) / (a * re)
    E² = E^2
    F  = (b * z + (a² - b²)) / (a * re)
    P  = (4 / 3) * (E * F + 1)
    Q  = 2 * (E² - F^2)
    D  = P^3 + Q^2

    if D ≥ 0
        aux = √D
        v = (aux - Q)^(1/3) - (Q + aux)^(1/3)
    else
        aux = √(-P)
        v = 2 * aux * cos(acos(Q / (P * aux)) / 3)
    end

    G = (√(E² + v) + E)/2

    # NOTE: Reference [5] appears to have an error in Eq. (13), where we must
    # have G^2 instead of G inside the square root. The correct version can be
    # seen here:
    #
    #   https://www.astro.uni.torun.pl/~kb/Papers/geod/Geod-BG.htm
    #
    t = √(G^2 + (F - v * G) / (2 * G - E)) - G

    # Compute the geodetic latitude and altitude.
    ϕ_gd = atan(a * (1 - t^2)/(2b * t))
    sin_ϕ_gd, cos_ϕ_gd = sincos(ϕ_gd)
    h = (re - a * t) * cos_ϕ_gd + (z - b) * sin_ϕ_gd

    return ϕ_gd, h
end

"""
    geodetic_to_geocentric(ϕ_gd::Number, h::Number)

Compute the geocentric latitude and radius from the geodetic latitude `ϕ_gd`
(-π/2, π/2) [rad] and height above the reference ellipsoid `h` \\[m] (WGS-84).
Notice that the longitude is the same in both geocentric and geodetic
coordinates.

!!! info
    The algorithm is based in **[4]**(p. 3).

# Returns

* Geocentric latitude [rad].
* Radius from the center of the Earth [m].

# References

- **[4]** ISO TC 20/SC 14 N (2011). Geomagnetic Reference Models.
"""
function geodetic_to_geocentric(ϕ_gd::Number, h::Number)
    # Auxiliary variables to decrease computational burden.
    sin_ϕ_gd, cos_ϕ_gd = sincos(ϕ_gd)
    sin²_ϕ_gd = sin_ϕ_gd^2
    e²_wgs84  = e_wgs84^2

    # Radius of curvature in the prime vertical [m].
    N    = a_wgs84/√(1 - e²_wgs84 * sin²_ϕ_gd )

    # Compute the geocentric latitude and radius from the Earth center.
    ρ    = (N + h) * cos_ϕ_gd
    z    = (N * (1 - e²_wgs84) + h) * sin_ϕ_gd
    r    = √(ρ^2 + z^2)
    ϕ_gc = asin(z / r)

    return ϕ_gc, r
end


text:
using BioinformaticsBISC195
using Test

@testset "BioinformaticsBISC195" begin

@testset "Using Strings" begin

    @testset "normalizeDNA" begin
        @test normalizeDNA("aatgn") == "AATGN"
        @test_throws Exception normalizeDNA("ZCA")
        @test_throws Exception normalizeDNA(42)
        c = normalizeDNA('C')
        @test c == "C"
        @test typeof(c) == String
    end # normalizeDNA

    @testset "composition" begin
        seq = rand(['A','T','G','C','N'], 20) |> join
        bc = composition(seq)
        @test bc isa Dict

        @test bc['A'] == count(x-> x == 'A', seq)
        @test bc['C'] == count(x-> x == 'C', seq)
        @test bc['G'] == count(x-> x == 'G', seq)
        @test bc['T'] == count(x-> x == 'T', seq)
        @test bc['N'] == count(x-> x == 'N', seq)

        bc = composition(lowercase(seq))

        @test bc['A'] == count(x-> x == 'A', seq)
        @test bc['C'] == count(x-> x == 'C', seq)
        @test bc['G'] == count(x-> x == 'G', seq)
        @test bc['T'] == count(x-> x == 'T', seq)
        @test bc['N'] == count(x-> x == 'N', seq)
    end # composition

    @testset "gc_content" begin
        @test gc_content("ANTG") == 1/3
        @test gc_content("cccggg") * 100 == 100.0
        @test gc_content("ATta") == 0.0
    end # gc_content

    @testset "complement" begin
        @test complement("ATTAN") == "TAATN"
        @test complement("gcta") == "CGAT"
        @test complement("nnnnnnn") == "NNNNNNN"
        @test_throws Exception complement("ABC")
    end # complement

    @testset "reverse_complement" begin
        @test reverse_complement("ATTAN") == "NTAAT"
        @test reverse_complement("gcta") == "TAGC"
        @test reverse_complement("nnnnnnn") == "NNNNNNN"
        @test_throws Exception reverse_complement("ABC")
    end # reverse_complement

    @testset "parse_fasta" begin
        #testpath = normpath(joinpath(@__DIR__, "..", "BioinformaticsBISC195/data"))
        testpath = normpath(joinpath(@__DIR__, "..", "data"))
        genomes = joinpath(testpath, "cov2_genomes.fasta")
        ex1_path = joinpath(testpath, "ex1.fasta")
        ex2_path = joinpath(testpath, "ex2.fasta")

        ex1 = parse_fasta(ex1_path)
        @test ex1 isa Tuple
        @test all(x-> x isa String, ex1[1])
        @test all(x-> x isa String, ex1[2])

        @test ex1[1] == ["ex1.1 | easy", "ex1.2 | multiline"]
        @test ex1[2] == ["AATTATAGC", "CGCCCCCCAGTCGGATT"]

        @test_throws Exception parse_fasta(ex2_path)

        cov2 = parse_fasta(genomes)
        @test length(cov2[1]) == 8
        @test length(cov2[2]) == 8
    end #parse_fasta

    @testset "fasta_header" begin
        @test fasta_header("M0002 |China|Homo sapiens") == ("M0002", "China", "Homo sapiens")
        @test fasta_header("Another sequence") == ("Another sequence",)
        @test fasta_header("headers| can | have| any number | of | info blocks") == ("headers", "can", "have", "any number", "of", "info blocks")
    end

    #test isDNA, kmercount, kmerdistance
    @testset "isDNA" begin
        @test isDNA("AcGt") == true
        @test isDNA("Ackk") == false #returns false bc the normalizeDNA is not used to fix it beforehand
        @test isDNA(normalizeDNA("Ackk")) == true # returns true because normalizeDNA converts the 'k's into 'N's, which is a valid ambiguous base
        @test isDNA("zzcGt") == false
    end

    @testset "kmercount" begin
        @test kmercount("ggg", 3) == Dict("GGG" => 1)
        @test kmercount("ATATATATA", 4) == Dict("TATA" => 3, "ATAT" => 3)
        @test_throws Exception kmercount("ACGGGqq", 4)
        @test_throws Exception kmercount("a", 2)
    end

    @testset "kmerdistance" begin
        @test kmerdistance(["AAA", "AGT"], ["AAA", "GGG"]) == 1 - 1/3
        @test kmerdistance(["AAA", "AGT"], ["TTT", "GGG"]) == 1
        @test kmerdistance(["AAA", "ATT"], ["AGT", "AAA", "GTGG"]) == 0.75
        @test kmerdistance(["GGT", "CAT"], ["CAT", "GGT"]) == 0
    end

    @testset "kmercollecting" begin
        ex1 = parse_fasta(normpath(joinpath(@__DIR__, "..", "data/ex1.fasta")))
        #ex1 = parse_fasta(joinpath(@__DIR__, "data", "ex1.fasta"))
        kc1 = kmercollecting(ex1[2], 3)
        # I changed these all to Sets because the order of a set does not matter for comparison
        v1 = Set(["AAT", "ATT", "TAT", "ATA", "TAG", "TTA", "AGC"])
        v2 = Set(["AGT", "CCC", "CGG", "ATT", "GAT", "GCC", "CAG", "GTC", "CGC", "GGA", "CCA", "TCG"])
        @test kc1[1] == v1
        @test kc1[2] == v2

        set2 = ("ACCGGTT", "AAA")
        kc2 = kmercollecting(set2, 3)
        v3 = Set(["GGT", "ACC", "GTT", "CCG", "CGG"])
        v4 = Set(["AAA"])
        @test kc2[1] == v3
        @test kc2[2] == v4
    end

    @testset "shortendates" begin
        asia = parse_fasta(normpath(joinpath(@__DIR__, "..", "data", "covgen_asiashort.fasta")))
        #asia = parse_fasta(joinpath(@__DIR__, "data", "covgen_asiashort.fasta"))
        @test shortendates(asia[1]) == ["2019-12", "2020-01", "2020-02", "2020-03"]

        dates1 = ["Should Work | 2031-04", "Should work|2000-12-01| extra", "should not work | 2020", "should work|2001-01"]
        @test shortendates(dates1) == ["2031-04", "2000-12", "2001-01"]

        dates2 = ["should not work| XX-XX", "should work|XXXX-XX", "should work   |  XXXXXXXextra"] # the function does not check validity of dates, so these should still work
        @test shortendates(dates2) == ["XXXX-XX", "XXXXXXX"]
    end

    @testset "monthlycomparisons" begin
        headers1 = ["Should Work | 2019-12", "Should work|2020-01-01| extra", "should not work | 2020", "should work|2020-02", "should work|2020-03-2"]
        seq1 = ["ACCTTGGA", "ACCCTGGA", "ACCCTGG", "ACCCTGC", "AACCTGC"]
        data1 = [headers1, seq1]
        @test monthlycomparison(data1[2][1], data1, 3) == [maximum(swscorematrix(data1[2][1], data1[2][2])), maximum(swscorematrix(data1[2][1], data1[2][3])), maximum(swscorematrix(data1[2][1], data1[2][4]))]

        headers2 = ["Should Work | 2019-12", "Should work|2020-01-01| extra", "should skip | 2020-01", "should work|2020-02", "should work|2020-03-2"]
        seq2 = ["CCATAGG", "CCATTGG", "CCATGG", "CAATGG", "CAATGT"]
        data2 = [headers2, seq2]
        @test monthlycomparison(data2[2][1], data2, 3) == [maximum(swscorematrix(data2[2][1], data2[2][2])), maximum(swscorematrix(data2[2][1], data2[2][4])), maximum(swscorematrix(data2[2][1], data2[2][5]))]

        # covgen_asiashort.fasta is a version of the covgen_asia file with only 4 data points listed, from month 2019-21 thru 2020-03
        asia = parse_fasta(normpath(joinpath(@__DIR__, "..", "data", "covgen_asiashort.fasta")))
        #asia = parse_fasta(joinpath(@__DIR__, "data", "covgen_asiashort.fasta"))

        # this line below can be uncommented to run the test, but for some reason it does not work on github (I believe it times out and sends a kill signal)
        # this line does run on my computer and the test passes, but I am going to leave it uncommented for the push so all the other tests can run
        #@test monthlycomparison(asia[2][1], asia, 3) == [maximum(swscorematrix(asia[2][1], asia[2][2])), maximum(swscorematrix(asia[2][1], asia[2][3])), maximum(swscorematrix(asia[2][1], asia[2][4]))]
        @test_throws Exception monthlycomparison(asia[2][1], asia, 4)
        @test_throws Exception monthlycomparison("Isn't real DNA", asia, 3)

        # this file has the same data as covgen_asiashort, but it skips data for Feb 2020, so it should also throw an error
        asiabad = parse_fasta(normpath(joinpath(@__DIR__, "..", "data", "covgen_asiabad.fasta")))
        #asiabad = parse_fasta(joinpath(@__DIR__, "data", "covgen_asiabad.fasta"))
        @test_throws Exception monthlycomparison(asiabad[2][1], asiabad, 3)
    end

end # strings

# @testset "Using BioSequences" begin

#     @testset "normalizeDNA" begin
#         @test normalizeDNA("aatgn") == dna"AATGN"
#         @test_throws Exception normalizeDNA("ZCA")
#         @test_throws Exception normalizeDNA(42)
#         c = normalizeDNA('C')
#         @test c == dna"c"
#         @test c isa LongSequence
#     end #  normalizeDNA

#     @testset "gc_content" begin
#         @test gc_content(dna"ANTG") == 0.25
#         @test gc_content(dna"cccggg") * 100 == 100.0
#         @test gc_content(dna"ATta") == 0.0
#     end #  composition

#     @testset "composition" begin
#         seq = rand(['A','T','G','C','N'], 20) |> join |> LongDNASeq
#         bc = composition(seq)

#         @test bc[DNA_A] == count(==(DNA_A), collect(seq))
#         @test bc[DNA_C] == count(==(DNA_C), collect(seq))
#         @test bc[DNA_G] == count(==(DNA_G), collect(seq))
#         @test bc[DNA_T] == count(==(DNA_T), collect(seq))
#         @test bc[DNA_N] == count(==(DNA_N), collect(seq))
#     end #  gc_content

#     @testset "complement" begin
#         @test complement(dna"ATTAN") == dna"TAATN"
#         @test complement(dna"gcta") == dna"CGAT"
#         @test complement(dna"nnnnnnn") == dna"NNNNNNN"
#     end #  complement

#     @testset "reverse_complement" begin
#         @test reverse_complement(dna"ATTAN") == dna"NTAAT"
#         @test reverse_complement(dna"gcta") == dna"TAGC"
#         @test reverse_complement(dna"nnnnnnn") == dna"NNNNNNN"
#     end #  reverse_complement

#     @testset "parse_fasta" begin
#         testpath = normpath(joinpath(@__DIR__, "..", "data"))
#         genomes = joinpath(testpath, "cov2_genomes.fasta")
#         ex1_path = joinpath(testpath, "ex1.fasta")
#         ex2_path = joinpath(testpath, "ex2.fasta")

#         ex1 = parse_fasta(ex1_path)
#         @test ex1 isa Tuple
#         @test all(x-> x isa String, ex1[1])
#         @test all(x-> x isa LongSequence, ex1[2])

#         @test ex1[1] == ["ex1.1 | easy", "ex1.2 | multiline"]
#         @test ex1[2] == [dna"AATTATAGC", dna"CGCCCCCCAGTCGGATT"]

#         @test_throws Exception parse_fasta(ex2_path)

#         cov2 = parse_fasta(genomes)
#         @test length(cov2[1]) == 8
#         @test length(cov2[2]) == 8
#     end # parse_fasta

# end # BioSequences

end # BioinformaticsBISC195


text:
"""
    solve!(solution::AbstractMovingWindowSolution, problem, solver::AbstractMovingWindowSolver; kwargs...) :: AbstractMovingWindowSolution

returns the [`AbstractMovingWindowSolution`](@ref) of an [`AbstractInitialValueProblem`](@ref).
"""
function NSDEBase.solve!(solution::AbstractMovingWindowSolution, problem::AbstractInitialValueProblem, solver::AbstractMovingWindowSolver; kwargs...)
    return solver(solution, problem; kwargs...)
end

"""
    solve(problem, solver::AbstractMovingWindowSolver; kwargs...) :: AbstractMovingWindowSolution

returns the [`AbstractMovingWindowSolution`](@ref) of an [`AbstractInitialValueProblem`](@ref).
"""
function NSDEBase.solve(problem::AbstractInitialValueProblem, solver::AbstractMovingWindowSolver; kwargs...)
    return solver(problem; kwargs...)
end


text:
module NveData

using DataStructures
using DataFrames
using CSV
using Dates
using DelimitedFiles
using Statistics

export watershed_info

export read_esri_raster
export raster2vec
export read_bil, read_bil!

export read_dbk_ind, subset_dbk_ind
export read_metadata, subset_metadata

export read_senorge_data
export aggregate_elevations, aggregate_elevations!
export senorge_info

export write_met_data
export write_update
export metadata_elevbands
export process_runoff_data

export init_dataset, update_dataset

export find_last_senorge

export read_old_data

include("utils.jl")
include("utils_raster.jl")
include("init_dataset.jl")
include("update_dataset.jl")

end

max_forks_repo_name: jmgnve/NveData


text:
using Test
using Solvent
using StaticArrays
using LinearAlgebra
using SparseArrays
using Random
using Printf

# Creates a Laplacian matrix based on the code from: http://math.mit.edu/~stevenj/18.303/lecture-10.html
# construct the (M+1)xM matrix D, not including the 1/dx factor
sdiff1(M) = [ [1.0 spzeros(1, M-1)]; spdiagm(1=>ones(M-1)) - I ]

# make the discrete -Laplacian in 2d, with Dirichlet boundaries
function Laplacian(Nx, Ny, Lx, Ly)
    dx = Lx / (Nx+1)
    dy = Ly / (Ny+1)
    Dx = sdiff1(Nx) / dx
    Dy = sdiff1(Ny) / dy
    Ax = Dx' * Dx
    Ay = Dy' * Dy
    return kron(spdiagm(0=>ones(Ny)), Ax) + kron(Ay, spdiagm(0=>ones(Nx)))
end

@testset "Solvent Laplacian: CG" begin
    T = Float64
    Lx = 1
    Ly = 1
    Nx = 10
    Ny = 10
    A = Laplacian(Nx, Ny, Lx, Ly)
    n, _ = size(A)
    @info "Size of matrix A: ($n, $n)"
    b = rand(T, n)
    x = rand(T, n)

    mulbyA!(y, x) = (y .= A * x)

    rtol = sqrt(eps(T))
    atol = eps(T)
    solver_type = ConjugateGradientMethod(M = 100)
    linearsolver = LinearSolver(
        mulbyA!,
        solver_type,
        x;
        rtol = rtol,
        atol = atol,
    )

    x0 = copy(x)
    linearsolve!(linearsolver, x, b)
    @test norm(A * x - b) < rtol * norm(A * x0 - b)

end

# Tests effectiveness of Jacobi preconditioner using specially designed matrix from:
# http://www.math.iit.edu/~fass/477577_Chapter_16.pdf
@testset "Solvent Preconditioned Jacobi: CG" begin
    T = Float64
    n = 1000
    dv0 = .5 .+ .√(1:n)
    dv1 = ones(n - 1)
    dv2 = ones(n - 100)
    A = spdiagm(0=>dv0) + spdiagm(1=>dv1) + spdiagm(-1=>dv1) +
        spdiagm(100=>dv2) + spdiagm(-100=>dv2)
    @info "Size of matrix A: ($n, $n)"
    b = ones(T, 1000)
    x = rand(T, n)

    mulbyA!(y, x) = (y .= A * x)

    rtol = sqrt(eps(T))
    atol = eps(T)
    solver_type = ConjugateGradientMethod(M = 100)
    identitysolver = LinearSolver(
        mulbyA!,
        solver_type,
        x;
        rtol = rtol,
        atol = atol,
    )
    jacobisolver = LinearSolver(
        mulbyA!,
        solver_type,
        x;
        pc_alg = Jacobi(),
        rtol = rtol,
        atol = atol,
    )

    x0 = copy(x)
    jacobiiters = linearsolve!(jacobisolver, x, b)
    @test norm(A * x - b) < rtol * norm(A * x0 - b)
    identityiters = linearsolve!(identitysolver, x0, b)
    @test jacobiiters < identityiters
end

max_stars_repo_name: CliMA/SystemSolvers

max_issues_repo_name: CliMA/SystemSolvers