text:
module BenchPerfConfigSweepsTests

include("test_aqua.jl")

end  # module BenchPerfConfigSweepsTests


text:
using Plots
outs = (out1,out2);

nrs = ([o isa Dict ? BM.parse_rgi(o[:rgi])[3] : -1 for o in out1],
       [o isa Dict ? BM.parse_rgi(o[:rgi])[3] : -1 for o in out2])

v = ([o isa Dict ? o[:vol_km3][:median] : NaN for o in out1],
     [o isa Dict ? o[:vol_km3][:median] : NaN for o in out2])

a = ([o isa Dict ? o[:area_km2] : NaN for o in out1],
     [o isa Dict ? o[:area_km2] : NaN for o in out2])

h = (1000*v[1]./a[1], 1000*v[2]./a[2])

va = ([o isa Dict ? o[:volarea_km3] : NaN for o in out1],
      [o isa Dict ? o[:volarea_km3] : NaN for o in out2])

hva = (1000*va[1]./a[1], 1000*va[2]./a[2])

ratioh = (h[1]./hva[1], h[2]./hva[2])

# radar/fitting
fls = filter(x->startswith(x, region_str), readdir("data/GlaThiDa_3.0/mhuss/"))
h_nrs = map(x->parse(Int, split(split(x,'.')[1],'_')[2]), fls)

p1 = Plots.scatter(hva[1], h[1], label="with radar", xlabel="mean h of vol-area (m)",
                   ylabel= "mean h BITE (m)", reuse=false, title="Region $region_str") #, ticks=:native
top = max(VAWTools.maximumnan(hva[1]), VAWTools.maximumnan(h[1]))
Plots.plot!(linspace(0, top), linspace(0, top), label="1:1")
Plots.plot!(linspace(0, top), linspace(0, 2*top), label="1:2")
Plots.plot!(linspace(0, top), linspace(0, 0.5*top), label="1:0.5")

p2 = Plots.scatter(hva[2], h[2], label="without radar", xlabel="mean h of vol-area (m)",
                   ylabel= "mean h BITE (m)", reuse=false) # , ticks=:native
top = VAWTools.maximumnan(hva[2])
Plots.plot!(linspace(0, top), linspace(0, top), label="1:1")
Plots.plot!(linspace(0, top), linspace(0, 2*top), label="1:2")
Plots.plot!(linspace(0, top), linspace(0, 0.5*top), label="1:0.5")
Plots.plot!(linspace(0, top), linspace(0, 5*top), label="1:5")
Plots.plot!(linspace(0, top), linspace(0, top/5), label="1:0.2")

#display(plot(p1, p2, layout=(2,1)))
display(p1); display(p2)

#savefig("figs/scatter-$(region_str)_$(split(dir_,'/')[end-1]).png")

### Find outliers above

    fac = 5.0
    outl = []
    outln = []
    for i=1:2
        inds = ratioh[i].>fac
        push!(outl, [(n,ii,o) for (n,ii,o) in zip(nrs[i][inds], (1:length(outs[i]))[inds], outs[i][inds])])
        push!(outln, nrs[i][inds])
    end
    @show length.(outl), length.(nrs)

    # find outliers below
    fac = 0.5
    outlb = []
    for i=1:2
        inds = ratioh[i].<fac
        push!(outlb, [(n,o) for (n,o) in zip(nrs[i][inds], outs[i][inds])])
    end
    @show length.(outlb), length.(nrs)

if false
    mode = 1
    ii = length(outl)
    update_cache=false
    nr, ind, radar, o = outl[mode][ii];
    @show nr, ii
    @show radar, h[mode][ind], hva[mode][ind]

    (gid, gl,gb,pp,pm,pn,pl,th0d,
     theta0, logposterior, logposterior1d, logprior, logprior_debug,
     loglikelihood, loglikelihood1d, fwdm_fn, fwdm1d_fn,
     pmcmc_defaults, fit_target) = BM.init(o[:rgi], fit_vars, fit_target, sigma_of_model, runtyp, update_cache=update_cache);
    fwdsol = fwdm_fn(o[:thetas_mode]);
    display(BM.plot_run(fwdsol))
    display(Plots.plot(gl, field=gl.dhdt.vals, reuse=false))
    display(Plots.plot(gb.fields1d_original[:dhdt], reuse=false))

end
nothing

# Region 7, mode 1 TOO BIG, mode 1
#  261 odd dhdt
#  437 ?
#  532 synthetic dhdt is bad: positive in ablation area?! -> fixed dhdt.  But still bad...; not fitting to IV does not improve it.
#  657 synthetic dhdt is bad: positive in ablation area?!
#  664 synth, smb>0
#  858 synth but ok dhdt ???
#  887 synth, smb>0
#  890 has no ablation area, not calving?!
# > 893 ??? IV problem?
#  895 synth, smb>0
#  896 no ablation
# > 902 no ablation, odd DEM, --> use other DEM
#  903 high SMB
#  932 SMB>0, iscalving==false
#  951 synthetic dhdt is bad: positive in ablation area?!
#  984 synthetic dhdt is bad: positive in ablation area?!
# > 1396 The fitting is bad here with the inferred dhdt opposite to the measured one?! IV?
# 1464 (Kronebreen) has no radar?! (only one point)
# 1471 fixed(?)
# 1476
# 1492 Julibreen
# 1541
# 1542
# 1558

ext: jl


text:
#Ensemble Kalman Inversion: specific structures and function definitions

"""
    Inversion <: Process

An ensemble Kalman Inversion process
"""
struct Inversion <: Process end

"""
    find_ekp_stepsize(ekp::EnsembleKalmanProcess{FT, IT, Inversion}, g::Array{FT, 2}; cov_threshold::FT=0.01) where {FT}

Find largest stepsize for the EK solver that leads to a reduction of the determinant of the sample
covariance matrix no greater than cov_threshold.
"""
function find_ekp_stepsize(
    ekp::EnsembleKalmanProcess{FT, IT, Inversion},
    g::Array{FT, 2};
    cov_threshold::FT = 0.01,
) where {FT, IT}
    accept_stepsize = false
    if !isempty(ekp.Δt)
        Δt = deepcopy(ekp.Δt[end])
    else
        Δt = FT(1)
    end
    # final_params [N_par × N_ens]
    cov_init = cov(get_u_final(ekp), dims = 2)
    while accept_stepsize == false
        ekp_copy = deepcopy(ekp)
        update_ensemble!(ekp_copy, g, Δt_new = Δt)
        cov_new = cov(get_u_final(ekp_copy), dims = 2)
        if det(cov_new) > cov_threshold * det(cov_init)
            accept_stepsize = true
        else
            Δt = Δt / 2
        end
    end

    return Δt

end

"""
     split_indices_by_success(output)

Obtain the successful/failed particle split. Failure, if any isnan present in the output along dimension 1.
"""
function split_indices_by_success(output)
    failed_ens = [i for i = 1:size(output,2) if any(isnan.(output[:,i]))]
    successful_ens = filter(x-> !(x in failed_ens), collect(1:size(output,2)))
    if length(failed_ens) > length(successful_ens)
            @warn string("More than 50% of runs produced NaN.",
                         "\nIterating... \nbut consider increasing model stability.",
                         "\nThis will effect optimization result.")
        end

    return successful_ens, failed_ens
end

function update_successful_ens(u,g,y,obs_noise_cov)
    #update successful ones
    cov_ug = cov(u, g, dims=2, corrected=false)
    cov_gg = cov(g, g, dims=2, corrected=false)

    tmp = (cov_gg + obs_noise_cov) \ (y - g)
    return u + (cov_ug * tmp) # [N_par × N_ens]
end

function update_failed_ens(u_old_fail, u_succ,fail_method)

        if fail_method == "shift"
            cov_u_new = cov(u_succ,u_succ, dims=2)
            #given that the param space is not too large these operations are fine, otherwise solve linear systems
            #get sqrt(cov_u_new)
            eval_new,evec_new = eigen(cov_u_new)
            sqrtcov_new = evec_new*Diagonal(sqrt.(eval_new))*evec_new'
            eval_init,evec_init = eigen(cov_init)
            invsqrtcov_init = evec_init*Diagonal(sqrt.(eval_init))*evec_init'

            return mean(u_succ) .+ sqrtcov_new*invsqrtcov_init*(u_old_fail .- mean(u_old,dims=2))
        elseif fail_method == "sample_succ_gauss"

            cov_u_new = cov(u_succ,u_succ, dims=2)
            mean_u_new = mean(u_succ,dims=2)

            return rand(MvNormal(mean_u_new[:],cov_u_new),size(u_old_fail,2))
        end

end

"""
    update_ensemble!(ekp::EnsembleKalmanProcess{FT, IT, <:Inversion}, g::Array{FT,2} cov_threshold::FT=0.01, Δt_new=nothing) where {FT, IT}

Updates the ensemble according to which type of Process we have. Model outputs `g` need to be a `N_obs × N_ens` array (i.e data are columms).
"""
function update_ensemble!(
    ekp::EnsembleKalmanProcess{FT, IT, Inversion},
    g::Array{FT, 2};
    cov_threshold::FT = 0.01,
    Δt_new = nothing,
    deterministic_forward_map = true,
    fail_method="sample_succ_gauss"
) where {FT, IT}

    # Update follows eqns. (4) and (5) of Schillings and Stuart (2017)

    #catch works when g non-square
    if !(size(g)[2] == ekp.N_ens)
        throw(DimensionMismatch("ensemble size in EnsembleKalmanProcess and g do not match, try transposing g or check ensemble size"))
    end

    #get successes and failures
    successful_ens, failed_ens = split_indices_by_success(g)

    # u: N_par × N_ens
    # g: N_obs × N_ens
    u = deepcopy(get_u_final(ekp))
    u_old = deepcopy(u)
    N_obs = size(g, 1)

    cov_init = cov(u, dims=2)
    cov_ug = cov(u,g, dims = 2, corrected=false) # [N_par × N_obs]
    cov_gg = cov(g,g, dims = 2, corrected=false) # [N_obs × N_obs]

    if !isnothing(Δt_new)
        push!(ekp.Δt, Δt_new)
    elseif isnothing(Δt_new) && isempty(ekp.Δt)
        push!(ekp.Δt, FT(1))
    else
        push!(ekp.Δt, ekp.Δt[end])
    end

    # Scale noise using Δt
    scaled_obs_noise_cov = ekp.obs_noise_cov / ekp.Δt[end]
    noise = rand(MvNormal(zeros(N_obs), scaled_obs_noise_cov), ekp.N_ens)

    # Add obs_mean (N_obs) to each column of noise (N_obs × N_ens) if
    # G is deterministic
    y = deterministic_forward_map ? (ekp.obs_mean .+ noise) : (ekp.obs_mean .+ zero(noise))

    #batch particle failures

    if length(failed_ens) == 0
        # N_obs × N_obs \ [N_obs × N_ens]
        # --> tmp is [N_obs × N_ens]
        tmp = (cov_gg + scaled_obs_noise_cov) \ (y - g)
        u += (cov_ug * tmp) # [N_par × N_ens]
    else

        u[:,successful_ens] = update_successful_ens(
            u[:,successful_ens],
            g[:,successful_ens],
            y[:,successful_ens],
            scaled_obs_noise_cov)

        u[:,failed_ens] = update_failed_ens(
            u[:,failed_ens],       # NB: at t^n
            u[:,successful_ens],   # NB: at t^{n+1}
            fail_method)

        println("successful mean: ", mean(u[:,successful_ens],dims=2))
        println("failed mean: ", mean(u[:,failed_ens],dims=2))
        println("overall mean: ", mean(u,dims=2))

    end

    diff = ekp.obs_mean .- g[:,successful_ens]
    X = ekp.obs_noise_cov \ diff # diff: column vector
    newerr = 1/size(diff)[2]*dot(diff, X)
    println("loss function", newerr)

    # store new parameters (and model outputs)
    push!(ekp.u, DataContainer(u, data_are_columns = true))
    push!(ekp.g, DataContainer(g, data_are_columns = true))

    # Store error
    compute_error!(ekp)

    # Check convergence
    cov_new = cov(get_u_final(ekp), dims = 2)
    cov_ratio = det(cov_new) / det(cov_init)
    if cov_ratio < cov_threshold
        @warn string(
            "New ensemble covariance determinant is less than ",
            cov_threshold,
            " times its previous value.",
            "\nConsider reducing the EK time step.",
        )
    end
end


text:
export units2ucum, validate_units, vucum
# Units issues
# (1) UCUM has not standardized the abbreviations of units used in information technology. For now, use "By" for Bytes and "bit" for bits, cf. https://hl7.org/fhir/2017Jan/valueset-ucum-units.html ; but neither appears in the most recent HL7 UCUM value set.

# Fixing a few common units mistakes that I keep seeing
const units_table = Dict{String,String}(
  "bits/sec" => "bit/s",
  "pascals" => "Pa",
  "degrees" => "deg",
  "reboots" => "{reboots}",
  "percent" => "%",
  "gaps" => "{gaps}",
  "bytes" => "By",
  "cycles" => "{cycles}",
  "counts" => "{counts}",
  "COUNTS" => "{counts}",
  "M**3/M**3" => "m3/m3",
  "M/S" => "m/s",
  "PA" => "Pa",
  "M/M" => "m/m",
  "M/S**2" => "m/s2",
  "M/S/S" => "m/s2",
  "M" => "m"
  )
fix_units(s::AbstractString) = get(units_table, s, s)

function strip_exp!(k::Array{Int64,1}, u::Array{UInt32,1}, j::Int64)
  (j == length(u)) && return j
  i = j

  # ^+, ^-, **+, **-
  if u[i] == 0x0000002b || u[i] == 0x0000002d
    push!(k, i)
    i = i+1
  end

  if i ≤ length(u)
    while i < length(u)
      if u[i] in 0x00000030:0x00000039 #'0':'9'
        i = i+1
      else # end of digits reached; break
        break
      end
    end
    # end
  end
  return i
end

# Home of a string function to convert units to SI
"""
    units2ucum(str)

Convert unit string `str` to UCUM syntax.

Rules applied to strings:
* Deletes "power" notation for integer powers: (^, ^+, ^-, **, **+, **-)
* Deletes "power" notation for special cases: (^+1, ^-1, **+1, **-1)
* Replaces multiplication symbols (*, ⋅, ×) with .
"""
function units2ucum(str::String)
  u = [UInt32(str[i]) for i in eachindex(str)]

  # set denominator position if '-' && no '/' to before the last alphabetical sequence
  # preceding
  m = 0
  d = false
  for i = 1:length(u)
    if u[i] == 0x0000002f
      d = false
      break
    elseif u[i] == 0x0000002d
      if m == 0
        m = i
      end
      d = true
    end
  end
  if d
    j = m
    f = false
    while j > 1
      j -= 1
      if (u[j] in 0x00000061:0x0000007a) || (u[j] in 0x00000041:0x0000005a)
        # start of alphabetical unit sequence
        if f == false
          f = true
        end
      elseif f == true
        insert!(u, j, 0x0000002f)
        break
      end
    end
  end

  i = 1
  k = Int64[]
  while i < length(u)

    # expressions that start with ^
    if u[i] == 0x0000005e
      #= Delete:
          ^
          ^+
          ^-
          ^+1
          ^-1
      =#
      push!(k, i)
      i = strip_exp!(k, u, i+1)

    # expressions that start with *
    elseif u[i] == 0x0000002a
      if i+1 ≤ length(u)
        #= Delete:
            **
            **+
            **-
            **+1
            **-1
        =#
        if u[i+1] == 0x0000002a
          push!(k, i)
          push!(k, i+1)
          i = strip_exp!(k, u, i+2)
        else
          # Change to "."
          u[i] = 0x0000002e
          if i > 1
            if u[i-1] == 0x0000002f
              push!(k, i)
            end
          end
          i = i+1
        end
      end

    # expressions that start with /, (space)
    elseif u[i] in (0x0000002f, 0x00000020)
      # delete whitespace after a / or space
      j = i
      while j < length(u)
        j = j+1
        if u[j] in (0x000000d7, 0x00000020, 0x0000002e, 0x000022c5)
          push!(k, j)
        else
          i = j-1
          break
        end
      end
      i = i+1
    else
      i = i+1
    end
  end
  deleteat!(u, k)

  # Replace multiplication operators with periods
  for i = 1:length(u)
    if u[i] in (0x000000d7, 0x00000020, 0x000022c5)
      u[i] = 0x0000002e
    end
  end

  # Get the position of the "/", if one exists
  m = 0
  for i = 1:length(u)
    if u[i] == 0x0000002f
      m = i
      break
    end
  end

  if m > 0
    k = Int64[]

    # delete "." in "./"
    if m > 1
      if u[m-1] == 0x0000002e
        push!(k, m-1)
      end
    end

    # delete all powers of "1" after the /
    for i = m+1:length(u)
      if u[i] == 0x00000031
        if i == length(u)
          push!(k, i)
          break
        elseif (u[i+1] in 0x00000030:0x00000039) == false
          push!(k, i)
        end
      end
    end

    deleteat!(u, k)
  end
  return String(Char.(u))
end

"""
    vucum(str::String)

Test whether `str` is a valid UCUM unit string.
"""
function vucum(str::String)
  rstr = identity(str)
  if any([occursin(c, rstr) for c in ('%', '^', '[', ']', '{', '}')])
    rstr = replace(rstr, "%" => "%25")
    rstr = replace(rstr, "^" => "%5E")
    rstr = replace(rstr, "[" => "%5B")
    rstr = replace(rstr, "]" => "%5D")
    rstr = replace(rstr, "{" => "%7B")
    rstr = replace(rstr, "}" => "%7D")
  end
  url = "https://ucum.nlm.nih.gov/ucum-service/v1/isValidUCUM/" * rstr
  tf = try
    req = request("GET", url)
    b = String(req.body)
    if b == "true"
      true
    else
      false
    end
  catch err
    @warn(string("Error thrown for unit string: ", str))
    false
  end
  return tf
end

"""
    validate_units(S::GphysData)

Test whether unit strings in S.units are valid under the UCUM standard.

    validate_units(C::GphysChannel)

Test whether C.units is valid under the UCUM standard.
"""
function validate_units(S::GphysData)
  isv = falses(S.n)
  for i = 1:S.n
    isv[i] = vucum(S.units[i])
  end
  return isv
end

validate_units(C::GphysChannel) = vucum(C.units)

ucum_to_seed = Dict{String,String}(
  "{counts}" => "counts",
  "m/s2" => "M/S**2"
)


text:
"""
    Conv(kernel, in=>out, σ; stride = 1, pading = 0, set_w = "Xavier")
This is the traditional convolution layer. `kernel` is a tuple of integers that specifies the kernel size, it must have one or two elements. And, `in` and `out` specifies number of input and out channels.

The input data must have a dimensions WHCB(width, height, channel, batch). If you want to use a data which has a dimentions WHC, you must be add a dimentioins of B.

# Parameters
- `stride`: specify stride of Convolution layer. This is a integer or tuple of 2 elements
- `padding`: specify padding of COnvolution layer. This is a integer or tuple of 2 or 4 elements. If you pecidy `KeepSize` for this parameter, we adjust the size of input and return a matrix which has the same sizes.
- `set_w`: `Xavier` or `He`, it decide a method to create a first parameter. This parameter is the same as `Dense()`.

# Example
```
julia> C = Conv((2, 2), 2=>2, relu)
Convolution(k:(2, 2), IO:2 => 2, σ:relu)

julia> C(rand(10, 10, 2, 5)) |> size
(9, 9, 2, 5)
```

!!! warning
    When you specidies `same` to `padding`, in some cases, it will be returned one size smaller.
    Because of its expression.
```
julia> C = Conv((2, 2), 2=>2, relu, padding = KeepSize)
Convolution(k:(2, 2), IO:2 => 2, σ:relu

julia> C(rand(10, 10, 2, 5)) |> size
(9, 9, 2, 5)
```
"""
struct Conv{N, M, K, I, O}
    σ
    weight
    bias
    stride::NTuple{N, Int}
    padding::NTuple{M, Int}
    dilation::NTuple{N, Int}
    high_accuracy::Bool
end

Conv(σ, weight, bias, stride, padding, dilation, high_accuracy) = Conv{length(stride), length(padding), size(weight)[1:2], size(weight)[3:4]...}(σ, weight, bias, stride, padding, dilation, high_accuracy)

function Conv(weight::AbstractArray{T, N}, σ, k::Tuple, io::Pair{<:Integer, <:Integer};
        stride = 1, padding = 0, dilation=1, high_accuracy::Bool=false) where {N, T}
    stride = expand(Val(N-2), stride)
    dilation = expand(Val(N-2), dilation)
    padding = expand_padding(padding, size(weight)[1:N-2], dilation, stride)
    bias = create_bias(weight, size(weight, N))
    return Conv{length(stride), length(padding), k, io...}(σ, weight, bias,
        stride, padding, dilation, high_accuracy)
end

function Conv(k::NTuple{N, Int}, io::Pair{<:Integer, <:Integer}, σ;
        stride = 1, padding = 0, dilation=1, set_w = "Xavier", high_accuracy::Bool=false) where {N}
    weight = conv_w(k, io, set_w, high_accuracy=high_accuracy)
    return Conv(weight, σ, k, io, stride = stride, padding = padding, dilation = dilation, high_accuracy=high_accuracy)
end

function Base.show(io::IO, C::Conv{N, M, K, I, O}) where {N, M, K, I, O}
    σ = C.σ
    print("Convolution(k:$K, IO:$I => $O, σ:$σ)")
end

function (C::Conv)(x)
    σ, b, w = C.σ, C.bias, C.weight
    dim = DenseConvDims(x, w, stride = C.stride, padding = C.padding, dilation = C.dilation)
    if !C.high_accuracy && isa(x, AbstractArray{Float64})
        x = Float32.(x)
    end
    return σ.(conv(x, w, dim).+b)
end

trainable(C::Conv) = C.weight, C.bias


text:
# ---
# title: 557. Reverse Words in a String III
# id: problem557
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/reverse-words-in-a-string-iii/description/>
# hidden: true
# ---
#
# Given a string, you need to reverse the order of characters in each word
# within a sentence while still preserving whitespace and initial word order.
#
# **Example 1:**
#
#
#
#     Input: "Let's take LeetCode contest"
#     Output: "s'teL ekat edoCteeL tsetnoc"
#
#
# **Note:** In the string, each word is separated by single space and there will
# not be any extra space in the string.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Global Right Hand Side function cache for dynamic rules.
# Now we're talking.
# TODO use a LRUCache?
const RHS_FUNCTION_CACHE = IdDict{Tuple{DynamicRule, Module}, Function}()
const RHS_FUNCTION_CACHE_LOCK = ReentrantLock()

function getrhsfun(r::DynamicRule, m::Module)
    lock(RHS_FUNCTION_CACHE_LOCK) do
        p = (r,m)
        if !haskey(RHS_FUNCTION_CACHE, p)
            z = genrhsfun(r, m)
            RHS_FUNCTION_CACHE[p] = z
        end
        return RHS_FUNCTION_CACHE[p]
    end
end

"""
Generates a [`RuntimeGeneratedFunction`](https://github.com/SciML/RuntimeGeneratedFunctions.jl/) corresponding to the right hand
side of a `:dynamic` [`Rule`](@ref).
"""
function genrhsfun(r::DynamicRule, mod::Module)
    params = Expr(:tuple, :_lhs_expr, :_subst, :_egraph, r.patvars...)
    ex = :($params -> $(r.right))
    closure_generator(mod, ex)
end


text:
# ------------------------------------------------------------------------------
# Local Fourier Analysis Toolkit
# ------------------------------------------------------------------------------

module LFAToolkit

# ------------------------------------------------------------------------------
# standard libraries
# ------------------------------------------------------------------------------

using LinearAlgebra
using Printf
using SparseArrays

# ------------------------------------------------------------------------------
# user available types and methods
# ------------------------------------------------------------------------------

export EvaluationMode, MultigridType, BDDCInjectionType
# mesh
export Mesh1D, Mesh2D, Mesh3D
# bases
export TensorBasis,
    NonTensorBasis,
    TensorH1LagrangeBasis,
    TensorH1UniformBasis,
    TensorMacroElementBasisFrom1D,
    TensorH1LagrangeMacroBasis,
    TensorH1UniformMacroBasis,
    TensorHProlongationBasis,
    TensorH1LagrangeHProlongationBasis,
    TensorH1UniformHProlongationBasis,
    TensorH1LagrangePProlongationBasis,
    TensorHProlongationMacroBasisFrom1D,
    TensorH1LagrangeHProlongationMacroBasis,
    TensorH1UniformHProlongationMacroBasis
# operato oields
export OperatorField
# operators
export Operator,
    GalleryOperator, GalleryVectorOperator, GalleryMacroElementOperator, computesymbols
# preconditioners
export IdentityPC
export NoPreco
export Jacobi
export Schwarz
export ASM
export ChebSchwarz, seteigenvalueestimatescalingchebschwarz
export Chebyshev, seteigenvalueestimatescaling
# -- multigrid
export Multigrid, PMultigrid, HMultigrid
# -- BDDC
export BDDC, LumpedBDDC, DirichletBDDC

# ------------------------------------------------------------------------------
# imports
# ------------------------------------------------------------------------------

include("Enums.jl")
# mesh
include("Mesh.jl")
# bases
include("Basis/Base.jl")
include("Basis/Constructors.jl")
# operator fields
include("OperatorField.jl")
# operators
include("Operator/Base.jl")
include("Operator/Constructors.jl")
# preconditioners
include("PC/Base.jl")
include("PC/Identity.jl")
include("PC/NoPreco.jl")
include("PC/Jacobi.jl")
include("PC/Schwarz.jl")
include("PC/ASM.jl")
include("PC/Chebyshev.jl")
include("PC/ChebSchwarz.jl")
# -- multigrid
include("PC/Multigrid/Base.jl")
include("PC/Multigrid/Constructors.jl")
# -- BDDC
include("PC/BDDC/Base.jl")
include("PC/BDDC/Constructors.jl")

end # module

# ------------------------------------------------------------------------------


text:
mutable struct Buffer
    uri::Union{Nothing,String}
    byteLength::Int
    name::Union{Nothing,String}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Buffer()
        obj = new()
        obj.uri = nothing
        obj.name = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.setproperty!(obj::Buffer, sym::Symbol, x)
    if sym === :byteLength
        x ≥ 1 || throw(ArgumentError("byteLength should be ≥ 1"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{Buffer}) = JSON3.Mutable()
JSON3.omitempties(::Type{Buffer}) = (:uri, :name, :extensions, :extras)


text:
#=NOTE:The purpose of these specifications is to create a container for holding
the building blocks of regression specifications=#

#holds a collection of specs for FM models in vector form
#NOTE: if this construct works well I may add it to Finometrics
struct FMSpecs{T<:Any}
  N::Ref{Int} #number of specs, reference construct allows for immutability
  aggfunc::Function
  specnames::Vector{String}
  yspecs::Vector{Symbol}
  xspecs::Vector{FMExpr}
  xnames::Vector{Vector{Symbol}}
  withinspecs::Vector{FMExpr}
  clusterspecs::Vector{FMExpr}
  results::Vector{T}
end

#creates the bare constructor
function FMSpecs(sizehint::NInt = nothing,
    ::Type{T} = FMLM;
    aggfunc::Function = (f(m::FMLM)::T=m) #default result function is the full regression model
  )  where T
  local specnames::Vector{String} = Vector{String}()
  local yspecs::Vector{Symbol} = Vector{Symbol}()
  local xspecs::Vector{FMExpr} = Vector{FMExpr}()
  local xnames::Vector{Vector{Symbol}} = Vector{Vector{Symbol}}()
  local withinspecs::Vector{FMExpr} = Vector{FMExpr}()
  local clusterspecs::Vector{FMExpr} = Vector{FMExpr}()
  local results::Vector{T} = Vector{T}()

  #makes minimal difference, but this could be handy
  specvector::Vector = [specnames, yspecs, xspecs, xnames, withinspecs, clusterspecs, results]
  if !isnothing(sizehint)
    (v::Vector->sizehint!(v, sizehint)).(specvector)
  end

  return FMSpecs{T}(Ref{Int}(0), aggfunc,  specvector...)
end

Base.show(io::IO, fs::FMSpecs) = print(io, "N: $(fs.N)
  specnames: $(fs.specnames)
  yspecs: $(fs.yspecs)
  xspecs: $(fs.xspecs)
  xnames: $(fs.xnames)
  withinspecs: $(fs.withinspecs)
  clusterspecs: $(fs.clusterspecs)
  results: $(fs.results)")

#this applies the results function to each specification
function computeFMLMresults!(dfs::Union{GroupedDataFrame, Vector{SubDataFrame}, Vector{DataFrame}},
    specs::FMSpecs; parallel::Bool=false, containsmissings::Bool=true)::Nothing where T<:AbstractDataFrame

  #make sure the dimensions are corred
  (specs.N[] == length(specs.specnames) &&
    specs.N[] == length(specs.yspecs) &&
    specs.N[] == length(specs.xspecs) &&
    specs.N[] == length(specs.xnames) &&
    specs.N[] == length(specs.withinspecs) &&
    specs.N[] == length(specs.clusterspecs)
    ) || error("Dimension mismatch error 43434 FMSpecs")

  (specs.N[] == length(dfs)) && "Input df views dimenion mismatch:
    dfs length: $(length(dfs)), specs.N[]: $(specs.N[])"

  #runs the regressions
  #could conceivably parallelize this at some point
  @mpar parallel for i::Int ∈ 1:specs.N[]
    m::FMLM = FMLM(dfs[i], specs.xspecs[i],  specs.yspecs[i],
      withinSym = specs.withinspecs[i], clusterSym = specs.clusterspecs[i],
      XNames=specs.xnames[i], YName = specs.yspecs[i],
      containsmissings=containsmissings)

    push!(specs.results, specs.aggfunc(m))
  end

  return nothing
end

#convenience method for providing a single dataframe
computeFMLMresults!(df::T where T<:AbstractDataFrame, specs::FMSpecs;
  parallel::Bool=false, containsmissings::Bool = true) =
  computeFMLMresults!(
    (i::Int->view(df, :, :)).(1:(specs.N[])), specs,
    parallel=parallel, containsmissings=containsmissings)

#provides a keyword access method for creating specs
function Base.push!(specs; specname::String = "($(specs.N[]+1))",
  yspec::Symbol = len > 1 ? specs.yspecs[end] : :Y,
  xspec::FMExpr = len > 1 ? specs.xspecs[end] : :X,
  xnames::Vector{Symbol} = len > 1 ? specs.xspecs[end] : [:intercept, :X],
  withinspec::FMExpr = nothing,
  clusterspec::FMExpr = nothing)

  #generate the specs
  push!(specs.specnames, specname)
  push!(specs.yspecs, yspec)
  push!(specs.xspecs, xspec)
  push!(specs.xnames, xnames)
  push!(specs.withinspecs, withinspec)
  push!(specs.clusterspecs, clusterspec)

  specs.N[] += 1

  return specs
end

Base.length(specs::FMSpecs) = specs.N[]

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/MortarContact2DAD.jl/blob/master/LICENSE

"""
2D mortar contact mechanics for JuliaFEM using Automatic Differentiation
"""
module MortarContact2DAD

using FEMBase, ForwardDiff, LinearAlgebra, SparseArrays, Statistics

const MortarElements2D = Union{Seg2,Seg3}

include("mortar2dad.jl")
include("contact2dad.jl")
export Mortar2DAD, Contact2DAD

end


text:
function U2NetP(in_ch=3, out_ch=1)
    stage1 = rsu7(in_ch, 16, 64)
    pool12 = MaxPool((2, 2),stride=2, pad=SamePad())

    stage2 = rsu6(64, 16, 64)
    pool23 = MaxPool((2, 2),stride=2, pad=SamePad())

    stage3 = rsu5(64, 16, 64)
    pool34 = MaxPool((2, 2),stride=2, pad=SamePad())

    stage4 = rsu4(64, 16, 64)
    pool45 = MaxPool((2, 2),stride=2, pad=SamePad())

    stage5 = rsu4f(64, 16, 64)
    MaxPool((2, 2),stride=2, pad=SamePad())

    stage6 = rsu4f(64, 16, 64)

    stage5d = rsu4f(128, 16, 64)
    stage4d = rsu4(128, 16, 64)
    stage3d = rsu5(128, 16, 64)
    stage2d = rsu6(128, 16, 64)
    stage1d = rsu7(128, 16, 64)

    side1 = Conv((3, 3), 64 => out_ch; pad=1)
    side2 = Conv((3, 3), 64 => out_ch; pad=1)
    side3 = Conv((3, 3), 64 => out_ch; pad=1)
    side4 = Conv((3, 3), 64 => out_ch; pad=1)
    side5 = Conv((3, 3), 64 => out_ch; pad=1)
    side6 = Conv((3, 3), 64 => out_ch; pad=1)

    outconv = side1 = Conv((1, 1), 6 => out_ch)

    layer = x -> begin
        hx1 = stage1(x)
        hx = pool12(hx1)

        hx2 = stage2(hx)
        hx = pool23(hx2)

        hx3 = stage3(hx)
        hx = pool34(hx3)

        hx4 = stage4(hx)
        hx = pool45(hx4)

        hx5 = stage5(hx)
        hx = pool56(hx5)

        hx6 = stage6(hx)
        hx6up = upsample(hx6,hx5)

        hx5d = stage5d(cat(hx6up,hx5; dims=3))
        hx5dup = upsample(hx5d,hx4)

        hx4d = stage4d(cat(hx5up,hx4; dims=3))
        hx4dup = upsample(hx4d,hx3)

        hx3d = stage3d(cat(hx4up,hx3; dims=3))
        hx3dup = upsample(hx3d,hx2)

        hx2d = stage2d(cat(hx3up,hx2; dims=3))
        hx2dup = upsample(hx2d,hx1)

        hx1d = stage1d(cat(hx2up,hx1; dims=3))

        d1 = side1(hx1d)

        d2 = side2(hx2d)
        d2 = upsample(d2,d1)

        d3 = side3(hx3d)
        d3 = upsample(d3,d1)

        d4 = side4(hx4d)
        d4 = upsample(d4,d1)

        d5 = side5(hx5d)
        d5 = upsample(d5,d1)

        d6 = side6(hx6)
        d6 = upsample(d6,d1)

        d0 = outconv(cat(d1,d2,d3,d4,d5,d6; dims=3))

        return σ.(d0), σ.(d1), σ.(d2), σ.(d3), σ.(d4), σ.(d5), σ.(d6)
    end
end


text:
using OrdinaryDiffEq, PyPlot, LinearAlgebra

const tol = 1e-8

function compute_max_lyapunov_exponent(rhs_fn, x0, params, δ0, Δt, N)
    # solve the reference trajectory
    problem = ODEProblem(rhs_fn, x0, (0, N*Δt), params)
    solver = AutoTsit5(Rosenbrock23())
    ref_sol = solve(problem, solver, reltol=tol, abstol=tol, dense=true)
    # neighboring trajectory
    du = ref_sol.u[2] - ref_sol.u[1] # vector along trajectory
    v = rand(size(du, 1))
    pert = v - dot(du, v) / dot(du, du) * du # vector orthogonal to trajectory
    pert *= δ0 / norm(pert) # make pertubation small
    x_test = x0 + pert
    # result
    s = 0.0
    for i in 1:N
        t = i*Δt
        # advance neighboring trajectory
        p = ODEProblem(rhs_fn, x_test, ((i - 1)*Δt, t), params)
        neighboring_sol = solve(p, solver,  reltol=tol, abstol=tol, dense=true)
        # compute seperation
        δ = norm(neighboring_sol(t) - ref_sol(t))
        # update result
        s += log(δ / δ0)
        # rescale the neighboring trajectory
        x_ref = ref_sol(t)
        x_test = x_ref + δ0 / δ * (x_test - x_ref)
    end
    # result
    return s / (Δt * N)
end

function lorentz_rhs(du, u, p, t)
    β, σ, ρ = p
    x, y, z = u
    du[1] = σ * (y - x)
    du[2] = - x * z - ρ * x - y
    du[3] = x * y - β * z
end

β = 8 / 3.
σ = 10.0
ρ_range = 20:0.1:100
λs = [compute_max_lyapunov_exponent(lorentz_rhs, [1.0, 1.0, 1.0],
           [β, σ, ρ], 1e-5, 0.001, 20)
     for ρ = ρ_range]
plot(ρ_range, λs)
xlabel("\$\\rho\$")
ylabel("\$\\lambda_1\$")

max_stars_repo_name: cosmo-jana/numerics-physics-stuff

max_issues_repo_name: cosmo-jana/numerics-physics-stuff

max_forks_repo_name: cosmo-jana/numerics-physics-stuff


text:
# All FiniteElement code dealing with implicit time integration.
# Directly included by FiniteElement.jl and inside module FiniteElement

"""
    function jac_assemble(solver, t)
Assembles the Jacobian matrix d(dg_rhs)/du at the current time step. Assumes
dg_rhs!() has just been called at this time step/stage and solver.u_interpV and
solver.u_interpF are precomputed.
"""
function jac_assemble(solver, t)
    if is_second_order(solver.app)
        # d∇u/du = (I⊗M) \ (d F_l(u)/du - Gs)
        compute_dgrad_du!(solver.dgrad_du, solver)
        # JD
    end
    jac_convect_volume(JD, solver)
    jac_convect_flux(JD, solver)
end
"""
    function compute_dgrad_du!(dgrad_du, solver)
Assembles the Jacobian matrix d(∇u)/du into dgrad_du. Computed as
d∇u/du = (I⊗M) \\ (d F_l(u)/du - Gs)
"""
function compute_dgrad_du!(dgrad_du, solver)
    for l=1:Geometry.DIM
        # TODO: Remap Gs to matrix array: 1 for each dimension
        dgrad_du .= -solver.Gels
        dlocal_flux_du!(dgrad_du, solver)
        # dgrad_du = M\dgrad_du
        ldiv!(solver.Mel, dgrad_du[l])
    end
end
"""
    function dlocal_flux_du!(dgrad_du, solver)
Updates dgrad_du += d(F_l(u))/du
"""
function dlocal_flux_du!(dgrad_du, solver)
    # TODO
end
"""
    function jac_convect_volume(solver)
Computes the jacobian of the convect_volume! function wrt u and adds it to JD,
the diagonal block matrix of the Jacobian so far.
"""
function jac_convect_volume(JD, solver)

end

max_stars_repo_name: NoseKnowsAll/DGToolkit

max_issues_repo_name: NoseKnowsAll/DGToolkit

max_forks_repo_name: NoseKnowsAll/DGToolkit


text:
# Julia wrapper for header: /usr/include/libavfilter/avfilter.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    avfilter_version,
    avfilter_configuration,
    avfilter_license,
    avfilter_pad_count,
    avfilter_pad_get_name,
    avfilter_pad_get_type,
    avfilter_link,
    avfilter_link_free,
    avfilter_link_get_channels,
    avfilter_link_set_closed,
    avfilter_config_links,
    avfilter_process_command,
    av_filter_iterate,
    avfilter_register_all,
    avfilter_register,
    avfilter_next,
    avfilter_get_by_name,
    avfilter_init_str,
    avfilter_init_dict,
    avfilter_free,
    avfilter_insert_filter,
    avfilter_get_class,
    avfilter_graph_alloc,
    avfilter_graph_alloc_filter,
    avfilter_graph_get_filter,
    avfilter_graph_create_filter,
    avfilter_graph_set_auto_convert,
    avfilter_graph_config,
    avfilter_graph_free,
    avfilter_inout_alloc,
    avfilter_inout_free,
    avfilter_graph_parse,
    avfilter_graph_parse_ptr,
    avfilter_graph_parse2,
    avfilter_graph_send_command,
    avfilter_graph_queue_command,
    avfilter_graph_dump,
    avfilter_graph_request_oldest

function avfilter_version()
    ccall((:avfilter_version, libavfilter), UInt32, ())
end

function avfilter_configuration()
    ccall((:avfilter_configuration, libavfilter), Cstring, ())
end

function avfilter_license()
    ccall((:avfilter_license, libavfilter), Cstring, ())
end

function avfilter_pad_count(pads)
    ccall((:avfilter_pad_count, libavfilter), Cint, (Ptr{AVFilterPad},), pads)
end

function avfilter_pad_get_name(pads, pad_idx::Integer)
    ccall((:avfilter_pad_get_name, libavfilter), Cstring, (Ptr{AVFilterPad}, Cint), pads, pad_idx)
end

function avfilter_pad_get_type(pads, pad_idx::Integer)
    ccall((:avfilter_pad_get_type, libavfilter), AVMediaType, (Ptr{AVFilterPad}, Cint), pads, pad_idx)
end

function avfilter_link(src, srcpad::Integer, dst, dstpad::Integer)
    ccall((:avfilter_link, libavfilter), Cint, (Ptr{AVFilterContext}, UInt32, Ptr{AVFilterContext}, UInt32), src, srcpad, dst, dstpad)
end

function avfilter_link_free(link)
    ccall((:avfilter_link_free, libavfilter), Cvoid, (Ptr{Ptr{AVFilterLink}},), link)
end

function avfilter_link_get_channels(link)
    ccall((:avfilter_link_get_channels, libavfilter), Cint, (Ptr{AVFilterLink},), link)
end

function avfilter_link_set_closed(link, closed::Integer)
    ccall((:avfilter_link_set_closed, libavfilter), Cvoid, (Ptr{AVFilterLink}, Cint), link, closed)
end

function avfilter_config_links(filter)
    ccall((:avfilter_config_links, libavfilter), Cint, (Ptr{AVFilterContext},), filter)
end

function avfilter_process_command(filter, cmd, arg, res, res_len::Integer, flags::Integer)
    ccall((:avfilter_process_command, libavfilter), Cint, (Ptr{AVFilterContext}, Cstring, Cstring, Cstring, Cint, Cint), filter, cmd, arg, res, res_len, flags)
end

function av_filter_iterate(opaque)
    ccall((:av_filter_iterate, libavfilter), Ptr{AVFilter}, (Ptr{Ptr{Cvoid}},), opaque)
end

function avfilter_register_all()
    ccall((:avfilter_register_all, libavfilter), Cvoid, ())
end

function avfilter_register(filter)
    ccall((:avfilter_register, libavfilter), Cint, (Ptr{AVFilter},), filter)
end

function avfilter_next(prev)
    ccall((:avfilter_next, libavfilter), Ptr{AVFilter}, (Ptr{AVFilter},), prev)
end

function avfilter_get_by_name(name)
    ccall((:avfilter_get_by_name, libavfilter), Ptr{AVFilter}, (Cstring,), name)
end

function avfilter_init_str(ctx, args)
    ccall((:avfilter_init_str, libavfilter), Cint, (Ptr{AVFilterContext}, Cstring), ctx, args)
end

function avfilter_init_dict(ctx, options)
    ccall((:avfilter_init_dict, libavfilter), Cint, (Ptr{AVFilterContext}, Ref{Ptr{AVDictionary}}), ctx, options)
end

function avfilter_free(filter)
    ccall((:avfilter_free, libavfilter), Cvoid, (Ptr{AVFilterContext},), filter)
end

function avfilter_insert_filter(link, filt, filt_srcpad_idx::Integer, filt_dstpad_idx::Integer)
    ccall((:avfilter_insert_filter, libavfilter), Cint, (Ptr{AVFilterLink}, Ptr{AVFilterContext}, UInt32, UInt32), link, filt, filt_srcpad_idx, filt_dstpad_idx)
end

function avfilter_get_class()
    ccall((:avfilter_get_class, libavfilter), Ptr{AVClass}, ())
end

function avfilter_graph_alloc()
    ccall((:avfilter_graph_alloc, libavfilter), Ptr{AVFilterGraph}, ())
end

function avfilter_graph_alloc_filter(graph, filter, name)
    ccall((:avfilter_graph_alloc_filter, libavfilter), Ptr{AVFilterContext}, (Ptr{AVFilterGraph}, Ptr{AVFilter}, Cstring), graph, filter, name)
end

function avfilter_graph_get_filter(graph, name)
    ccall((:avfilter_graph_get_filter, libavfilter), Ptr{AVFilterContext}, (Ptr{AVFilterGraph}, Cstring), graph, name)
end

function avfilter_graph_create_filter(filt_ctx, filt, name, args, opaque, graph_ctx)
    ccall((:avfilter_graph_create_filter, libavfilter), Cint, (Ptr{Ptr{AVFilterContext}}, Ptr{AVFilter}, Cstring, Cstring, Ptr{Cvoid}, Ptr{AVFilterGraph}), filt_ctx, filt, name, args, opaque, graph_ctx)
end

function avfilter_graph_set_auto_convert(graph, flags::Integer)
    ccall((:avfilter_graph_set_auto_convert, libavfilter), Cvoid, (Ptr{AVFilterGraph}, UInt32), graph, flags)
end

function avfilter_graph_config(graphctx, log_ctx)
    ccall((:avfilter_graph_config, libavfilter), Cint, (Ptr{AVFilterGraph}, Ptr{Cvoid}), graphctx, log_ctx)
end

function avfilter_graph_free(graph)
    ccall((:avfilter_graph_free, libavfilter), Cvoid, (Ptr{Ptr{AVFilterGraph}},), graph)
end

function avfilter_inout_alloc()
    ccall((:avfilter_inout_alloc, libavfilter), Ptr{AVFilterInOut}, ())
end

function avfilter_inout_free(inout)
    ccall((:avfilter_inout_free, libavfilter), Cvoid, (Ptr{Ptr{AVFilterInOut}},), inout)
end

function avfilter_graph_parse(graph, filters, inputs, outputs, log_ctx)
    ccall((:avfilter_graph_parse, libavfilter), Cint, (Ptr{AVFilterGraph}, Cstring, Ptr{AVFilterInOut}, Ptr{AVFilterInOut}, Ptr{Cvoid}), graph, filters, inputs, outputs, log_ctx)
end

function avfilter_graph_parse_ptr(graph, filters, inputs, outputs, log_ctx)
    ccall((:avfilter_graph_parse_ptr, libavfilter), Cint, (Ptr{AVFilterGraph}, Cstring, Ptr{Ptr{AVFilterInOut}}, Ptr{Ptr{AVFilterInOut}}, Ptr{Cvoid}), graph, filters, inputs, outputs, log_ctx)
end

function avfilter_graph_parse2(graph, filters, inputs, outputs)
    ccall((:avfilter_graph_parse2, libavfilter), Cint, (Ptr{AVFilterGraph}, Cstring, Ptr{Ptr{AVFilterInOut}}, Ptr{Ptr{AVFilterInOut}}), graph, filters, inputs, outputs)
end

function avfilter_graph_send_command(graph, target, cmd, arg, res, res_len::Integer, flags::Integer)
    ccall((:avfilter_graph_send_command, libavfilter), Cint, (Ptr{AVFilterGraph}, Cstring, Cstring, Cstring, Cstring, Cint, Cint), graph, target, cmd, arg, res, res_len, flags)
end

function avfilter_graph_queue_command(graph, target, cmd, arg, flags::Integer, ts::Cdouble)
    ccall((:avfilter_graph_queue_command, libavfilter), Cint, (Ptr{AVFilterGraph}, Cstring, Cstring, Cstring, Cint, Cdouble), graph, target, cmd, arg, flags, ts)
end

function avfilter_graph_dump(graph, options)
    ccall((:avfilter_graph_dump, libavfilter), Cstring, (Ptr{AVFilterGraph}, Cstring), graph, options)
end

function avfilter_graph_request_oldest(graph)
    ccall((:avfilter_graph_request_oldest, libavfilter), Cint, (Ptr{AVFilterGraph},), graph)
end


text:
module MQTT

import Base: connect, read, write, get

include("utils.jl")
include("packet.jl")
include("packets/connect.jl")
include("packets/publish.jl")
include("packets/subscribe.jl")
include("packets/unsubscribe.jl")
include("packets/disconnect.jl")
include("packets/ping.jl")
include("net.jl")
include("client.jl")

export
AT_MOST_ONCE,
AT_LEAST_ONCE,
EXACTLY_ONCE,
Client,
ConnectOpts,
get,
connect,
disconnect,
subscribe,
unsubscribe,
publish

end

max_issues_repo_name: kivaari/jlMQTT

max_forks_repo_name: kivaari/jlMQTT


text:
using Mathepia
using Test

@testset "Mathepia.jl" begin
    @test media_impact_fun(1,2) == exp(-2)
end


text:
using Pope, ZMQ, Base.Test

# prep for listening to the ZMQDataSink
Pope.init_for_zmqdatasink(Pope.ZMQ_PORT)

sub_socket = ZMQ.Socket(Pope.ZMQDataSinkConsts.context, ZMQ.SUB)
ZMQ.subscribe(sub_socket)
ZMQ.connect(sub_socket, "tcp://localhost:$(Pope.ZMQ_PORT)")
function recv_multipart(socket)
  out = ZMQ.Message[]
  push!(out, recv(socket))
  while ZMQ.ismore(socket)
    push!(out, recv(socket))
  end
  out
end
function recv_dataproduct(socket)
  out = recv_multipart(socket)
  channel_string = String(take!(convert(IOStream, out[1])))
  channel_number = parse(Int, channel_string)
  dp = read(seekstart(convert(IOStream,out[2])),Pope.MassCompatibleDataProductFeb2017,1)[1]
  return channel_string, channel_number, dp
end

zds = Pope.make_zmqdatasink(1)
dp = Pope.MassCompatibleDataProductFeb2017(1,2,3,4,5,6,7,8,9,10,11,12,13)
Pope.write_header(zds,nothing)
out_header = recv_multipart(sub_socket)
Pope.write_header_allchannel(zds,nothing)
out_header_allchannel = recv_multipart(sub_socket)
write(zds,dp)
channel_string, channel_number, dp_out = recv_dataproduct(sub_socket)
Pope.write_header_end(zds,nothing,nothing)
out_header_end = recv_multipart(sub_socket)

@testset "zmq_datasink" begin
  @test repr(zds.channel_number) == channel_string
  @test dp_out == dp
end

# test multiple sink
datasinks = Pope.MultipleDataSink((Pope.make_zmqdatasink(1),Pope.make_zmqdatasink(2)))
@testset "MultipleDataSink" begin
  write(datasinks, dp)
  channel_string1, channel_number1, dp_out1 = recv_dataproduct(sub_socket)
  channel_string2, channel_number2, dp_out2 = recv_dataproduct(sub_socket)
  @test channel_string1 == repr(datasinks.t[1].channel_number)
  @test channel_string2 == repr(datasinks.t[2].channel_number)
  @test dp_out1 == dp
  @test dp_out2 == dp
  Pope.write_header(datasinks, nothing)
  Pope.write_header_end(datasinks, nothing,nothing)
  Pope.write_header_allchannel(datasinks, nothing)
end

close(zds.s) # close the socket directly, close(zds) is currently a no-op
close(sub_socket)


text:
function CalibrationErrors.unsafe_skce_eval_targets(
    κtargets::Kernel, p::AbstractMixtureModel, y, p̃::AbstractMixtureModel, ỹ
)
    p_components = components(p)
    p_probs = probs(p)
    p̃_components = components(p̃)
    p̃_probs = probs(p̃)

    probsi = p_probs[1]
    componentsi = p_components[1]
    s =
        probsi *
        p̃_probs[1] *
        CalibrationErrors.unsafe_skce_eval_targets(
            κtargets, componentsi, y, p̃_components[1], ỹ
        )
    for j in 2:length(p̃_components)
        s +=
            probsi *
            p̃_probs[j] *
            CalibrationErrors.unsafe_skce_eval_targets(
                κtargets, componentsi, y, p̃_components[j], ỹ
            )
    end

    for i in 2:length(p_components)
        probsi = p_probs[i]
        componentsi = p_components[i]

        for j in 2:length(p̃_components)
            s +=
                probsi *
                p̃_probs[j] *
                CalibrationErrors.unsafe_skce_eval_targets(
                    κtargets, componentsi, y, p̃_components[j], ỹ
                )
        end
    end

    return s
end


text:
"""
utility functions helping managing result list
"""
module ResultListUtils
using CUDA,..MetaDataUtils, ..CUDAAtomicUtils
export getResLinIndex,allocateResultLists,@addResult

"""
allocate memory on GPU for storing result list
    totalFpCount-  total number of false positives
    TotalFNCount - total number of false negatives

    in the array first 3 entries will be x,y,z than isGold - 1 if it is  related to gold pass dilatations
        , direction from which result was set and the iteration number in which it was covered
    in the second list we have  the UInt32 Ids generated by the  getResLinIndex function
"""
function allocateResultLists(totalFpCount,TotalFNCount)
  return CUDA.zeros(UInt16, (totalFpCount+ TotalFNCount+1),6 )
end#allocateResultList

 """
 adding result to the result list at correct postion using data from metaData - so we get from the metadata offset and result counter
 metadata - 4 dimensional array holding metaData
 xMeta,yMeta,zMeta -  x,y,z coordinates of block of intrest in meta Data
 resList - list of result (matrix to be more precise) where we will wrtie the results
 resListIndicies - list of indicies related to results
 x,y,z - coordinates where we found point of intrest
 dir - direction from which dilatation covering this voxel had happened
 queueNumber - what fp or fn queue we are intrested in modyfing now
 metaDataDims - dimensions of metadata array
mainArrDims - dimensions of main array
 isGold - indicated is this a gold dilatation step (then it will evaluate to 1 otherwise 0 )
 """
 macro addResult(metaData ,xMeta,yMeta,zMeta, resList,x,y,z, dir,iterNumb,queueNumber,metaDataDims,mainArrDims ,isGold  )
  return esc(quote

# linearIndex =
  resListPos = ($metaData[($xMeta+1),($yMeta+1),($zMeta+1), (getResOffsetsBeg() +$queueNumber) ]+atomicallyAddToSpot( metaData,(($xMeta+1) + ($yMeta)*$metaDataDims[1] + ($zMeta)*$metaDataDims[1]*$metaDataDims[2] + (getNewCountersBeg()+$queueNumber-1)*$metaDataDims[1]*$metaDataDims[2]*$metaDataDims[3]),UInt32(1) ))+1
# qn = $queueNumber
# xm = $xMeta
# ym = $yMeta
# zm = $zMeta
# xx = $x
# yy= $y
# zz=$z
# dd= $dir

#  CUDA.@cuprint "\n resListPos $(resListPos) queueNumber $(qn)  xMeta $(xm) yMeta $(ym)  zMeta $(zm) x $(xx) y $(yy) z $(zz)  dir $(Int64(dd))  \n "

@inbounds $resList[ resListPos, 1]=$x
@inbounds $resList[ resListPos, 2]=$y
@inbounds $resList[ resListPos, 3]=$z
@inbounds $resList[ resListPos, 4]= $isGold
@inbounds $resList[ resListPos, 5]= $dir
@inbounds $resList[ resListPos, 6]= $iterNumb
# @inbounds $resListIndicies[resListPos]=getResLinIndex($x,$y,$z,$isGold,$mainArrDims)
# CUDA.@cuprint "\n linIndex $(getResLinIndex(x,y,z,isGold,mainArrDims))  \n "
#addResHelper(resListIndicies,resListPos,x)
 end)#quote
 end#addResult

 """
 giver the result row that holds data about covered point and in what iteration, from what direction and in what pass it was covered
 resRow - array where first 3 entries are x,y,z positions then is gold,
 """
  function getResLinIndex(x,y,z,isGold,mainArrDims)::UInt32
     # last one is in order to differentiate between gold pass and other pass ...
     return x+ y*mainArrDims[1]+ z* mainArrDims[1]*mainArrDims[2]+ isGold*mainArrDims[1]*mainArrDims[2]*mainArrDims[3]
   end#getResLinIndex

end#ResultListUtils


text:
module LicenseCheck

using licensecheck_jll: licensecheck_jll

export licensecheck, is_osi_approved
export find_licenses, find_license
export find_licenses_by_bruteforce, find_licenses_by_list,
       find_licenses_by_list_intersection

include("OSI_LICENSES.jl")
include("find_licenses.jl")

"""
    licensecheck(text::String) -> @NamedTuple{licenses_found::Vector{String}, license_file_percent_covered::Float64}

Returns a vector of the names of licenses (more specifically, the SPDX 3.10 license identifiers) matched in `text` and the percent of the text covered by these matches.

The full list of license IDs is located at [https://github.com/google/licensecheck/blob/v0.3.1/licenses/README.md](https://github.com/google/licensecheck/blob/v0.3.1/licenses/README.md), and includes the SDPX 3.10 licenses, [https://github.com/spdx/license-list-data/tree/v3.10/text](https://github.com/spdx/license-list-data/tree/v3.10/text).

This provides some of the functionality of `licensecheck.Scan` in the original Go library ([https://github.com/google/licensecheck](https://github.com/google/licensecheck)).

## Example

```julia
julia> using LicenseCheck

julia> text = read(joinpath(pkgdir(LicenseCheck), "LICENSE"), String);

julia> licensecheck(text)
(licenses_found = ["MIT"], license_file_percent_covered = 98.82352941176471)
```
"""
function licensecheck(text::String)
    arr, dims, license_file_percent_covered = ccall((:License,
                                                     licensecheck_jll.licensecheck),
                                                    Tuple{Ptr{Ptr{UInt8}},Cint,Float64},
                                                    (Cstring,), text)
    return (; licenses_found=unsafe_string.(unsafe_wrap(Array, arr, dims; own=true)),
            license_file_percent_covered=license_file_percent_covered)
end

"""
    is_osi_approved(spdx_identifier::String) -> Bool
    is_osi_approved(nt::NamedTuple) -> Bool

Checks if a [SDPX 3.10 identifier](https://spdx.dev/ids/) is
[OSI approved](https://opensource.org/licenses) by checking against an autogenerated list,
or if a `NamedTuple` with a `licenses_found` key contains only OSI-approved SDPX 3.10 identifiers.

!!! warning

    Note that [`licensecheck`](@ref) understands a larger number of licenses than
    SDPX 3.10, and can output results which are not SDPX 3.10 identifiers.
    Such licenses will yield `is_osi_approved(name) == false` regardless of
    their OSI standing.

## Example

```julia
julia> is_osi_approved("MIT")
true

julia> is_osi_approved(find_license(pkgdir(LicenseCheck)))
true

julia> is_osi_approved(find_license(homedir(), allow_brute=false)) # no license there
false

```
"""
is_osi_approved(spdx_identifier::String) = spdx_identifier ∈ OSI_LICENSES
function is_osi_approved(nt::NamedTuple)
    return !isempty(nt.licenses_found) && all(is_osi_approved, nt.licenses_found)
end
is_osi_approved(::Nothing) = false # so that it can always be used with `find_license`

end # module


text:
using Test

include("testtypevers.jl")
include("testutilapi.jl")
include("testwellknown.jl")
include("testcodec.jl")
include("services/testsvc.jl")
include("testprotoc.jl")


text:
#####
##### Competing Poisson processes
#####

export CompetingPoisson

struct CompetingPoisson{T,D,E}
    total_rate::T
    index_distribution::D
    events::E
end

"""
$(SIGNATURES)

The first arrival from competing Poisson point processes with the given `rates`.

`rand` returns a `NamedTuple{(:duration,:event)}`, where `duration` is the duration until
the first event and `event` is drawn from `events`.

Varios properties are supported, see [`propertynames`](@ref).

# Note

The distribution of `duration` is `Exponential(1/sum(rates))`, so this process can also be
thought of as competing exponentials.
"""
function CompetingPoisson(rates::AbstractVector; events = axis(rates, 1))
    total_rate = sum(rates)
    CompetingPoisson(total_rate, Categorical(rates ./ total_rate), events)
end

function Base.rand(rng::AbstractRNG, sampler::SamplerTrivial{<:CompetingPoisson})
    @unpack total_rate, index_distribution, events = sampler[]
    τ = rand(rng, Exponential(1 / total_rate))
    i = rand(rng, index_distribution)
    (duration = τ, event = events[i])
end

function Base.propertynames(cp::CompetingPoisson, private = false)
    public = (:total_rate, :rates, :probabilities, :events)
    private ? (:index_distribution, public...) : public
end

function Base.getproperty(cp::CompetingPoisson, key::Symbol)
    if key ≡ :probabilities
        probs(getfield(cp, :index_distribution))
    elseif key ≡ :rates
        cp.probabilities .* cp.total_rate
    else
        getfield(cp, key)
    end
end

function Base.show(io::IO, cp::CompetingPoisson)
    print(io, "Competing Poisson with total rate $(round(cp.total_rate; sigdigits = 3));")
    for (e, p) in zip(cp.events, cp.probabilities)
        print(io, " P($(e))=$(round(p; sigdigits = 3))")
    end
end


text:
## modeled after sort tests in python's sortperf.py
##
## Kevin Squire
##

module SortPerf

export sortperf, sort_plots, view_sort_plots, save_sort_plots, std_sort_tests, sort_median, sort_scale

import Base.Sort: Algorithm, Forward, ReverseOrdering, ord
import Base.Order.Ordering

using SortingAlgorithms
using DataFrames
using Winston
using Color

# rand functions for testing
randstr(n::Int) = [randstring() for i = 1:n]
randint(n::Int) = rand(1:n,n)

randfns = Dict(Int => randint,
               Int32 => x->round(Int32,randint(x)),
               Int64 => x->round(Int64,randint(x)),
               Int128 => x->round(Int128,randint(x)),
               Float32 => x->round(Float32,rand(x)),
               Float64 => rand,
               String => randstr)

std_types = [Int, Float64, String]
sort_algs = [InsertionSort, HeapSort, MergeSort, QuickSort, RadixSort, TimSort] #, TimSortUnstable]

# DataFrame labels
labels = [symbol("test_type"), symbol("sort_alg"), symbol("log_size"), symbol("size"), symbol("*sort"),
          symbol("\\sort"), symbol("/sort"), symbol("3sort"), symbol("+sort"), symbol("~sort"), symbol("=sort"), symbol("!sort")]

# Corresponding descriptions
sort_descr = Dict( '*' => "random",
                   '\\' => "reversed",
                   '/' => "sorted",
                   '3' => "3 exchanges",
                   '+' => "10 appended",
                   '~' => "4 unique items",
                   '=' => "all equal",
                   '!' => "qsort median killer")

# Test algorithm performance on a data vector
function _sortperf(alg::Algorithm, origdata::Vector, order::Ordering=Forward; replicates=3, skip_median_killer=false)
    srand(1)
    times = Dict[]   # Array of Dicts!
    n = length(origdata)
    logn = log2(length(origdata))

    rec = Dict(symbol("test_type") => string(eltype(origdata)),
               symbol("sort_alg") => string(alg)[1:end-5],
               symbol("log_size") => isinteger(logn) ? round(Int, logn) : logn,
               symbol("size") => length(origdata))

    if alg==RadixSort && !isbits(origdata[1])
        println("Skipping RadixSort for non bitstype $(eltype(origdata))")
        reptimes = Dict{Symbol, Float64}()
        for label in labels[5:end]
            reptimes[label] = 0.0
        end
        return times
    end

    for rep = 1:replicates
        print(" $rep")
        reptimes = Dict{Symbol, Float64}()

        ## Random
        data = copy(origdata)
        gc()
        if !issorted(data, order)
            reptimes[symbol("*sort")] = @elapsed sort!(data, alg, order)
            @assert issorted(data, order)
        end

        ## Sorted
        gc()
        reptimes[symbol("/sort")] = @elapsed sort!(data, alg, order)
        @assert issorted(data, order)

        ## Reverse sorted
        reverse!(data)
        @assert issorted(data, ReverseOrdering(order))
        gc()
        reptimes[symbol("\\sort")] = @elapsed sort!(data, alg, order)
        @assert issorted(data, order)

        ## Sorted with 3 exchanges
        for i = 1:3
            n1 = rand(1:n)
            n2 = rand(1:n)
            data[n1], data[n2] = data[n2], data[n1]
        end
        gc()
        reptimes[symbol("3sort")] = @elapsed sort!(data, alg, order)
        @assert issorted(data, order)

        ## Sorted with 10 unsorted values at end
        if length(data) >= 20
            for i = 1:10
                val = splice!(data, i:i)
                append!(data, val)
            end
            gc()
            reptimes[symbol("+sort")] = @elapsed sort!(data, alg, order)
            @assert issorted(data, order)
        end

        ## Random data with 4 unique values
        idxs = Int[]
        for i = 1:4
            idx = rand(1:n)
            while idx in idxs
                idx = rand(1:n)
            end
            push!(idxs, idx)
        end
        vals = data[idxs]
        data4 = vals[rand(1:4, n)]
        gc()
        reptimes[symbol("~sort")] = @elapsed sort!(data4, alg, order)
        @assert issorted(data, order)

        ## All values equal
        data1 = data[fill(rand(1:n), n)]
        gc()
        reptimes[symbol("=sort")] = @elapsed sort!(data1, alg, order)
        @assert issorted(data1, order)

        ## quicksort median killer: first half descending, second half ascending
        if skip_median_killer || alg==QuickSort && !(eltype(data) <: Integer) && length(data) >= 2^18
            print("(median killer skipped) ")
        else
            last = (length(data)>>1)<<1 # make sure data length is even
            qdata = vcat(data[last:-2:2], data[1:2:last])
            gc()
            reptimes[symbol("!sort")] = @elapsed sort!(qdata, alg, order)
            @assert issorted(qdata, order)
        end

        push!(times, merge(rec, reptimes))
    end
    println()

    times
end

# run sortperf over a range of data lengths
function _sortperf(alg::Algorithm, data::Vector, log2range::Range, order::Ordering=Forward; named...)
    lg2range = filter(x -> 2^x <= length(data), collect(log2range))
    times = Dict[]
    for logsize in lg2range
        println("  $logsize")
        size = 2^logsize
        if (alg == InsertionSort && logsize >= 14)
            println("Skipped")
            continue
        end

        append!(times, _sortperf(alg, data[1:size], order; named...))
    end

    times
end

# Generate a random array and run sortperf
function _sortperf(alg::Algorithm, T::Type, size::Int, args...; named...)
    println("Testing $T...")
    _sortperf(alg, randfns[T](size), args...; named...)
end

# Generate a random array and run sortperf over different ranges of that array
function _sortperf(alg::Algorithm, T::Type, log2range::Range, args...; named...)
    println("Testing $T...")
    _sortperf(alg, randfns[T](2^last(log2range)), log2range, args...; named...)
end

# Run sortperf for a number of types
function _sortperf(alg::Algorithm, types::Vector{DataType}, args...; named...)
    times = Dict[]
    for T in types
        append!(times, _sortperf(alg, T, args...; named...))
    end
    times
end

# Run sortperf on a number of algorithms
function _sortperf(algs::Vector{Algorithm}, args...; named...)
    times = Dict[]
    for alg in algs
        println("\n$alg\n")
        append!(times, _sortperf(alg, args...; named...))
    end
    times
end

# Returns a DataFrame version of sortperf output
sortperf(args...; named...) = DataFrame(_sortperf(args...; named...), labels)

# Get median sort timings
sort_median(df::DataFrame) = groupby(df, [:log_size, :size, :sort_alg, :test_type]) |> median

sort_scale(df::DataFrame, base_sort) = by(df, ["log_size", "size", "test_type"],
                                          x->(row = find(x["sort_alg"] .== base_sort);
                                              ht = size(x,1);
                                              hcat(x[:,3:3],x[:,5:end]./vcat(rep(x[row,5:end],ht)...))))

# Create sort plots
function sort_plots(df, base_sort, cols = ["*sort_median", "\\sort_median", "/sort_median", "3sort_median",
                                           "+sort_median", "~sort_median", "=sort_median", "!sort_median"])
    plots = PlotContainer[]
    sort_algs = unique(df["sort_alg"])
    dc = distinguishable_colors(12)[[1,3,4,5,6,7,8,9,10,11,12]]
    colors = Dict([symbol(alg) for alg in sort_algs], Uint32[convert(RGB24, r) for r in dc])
    linestyles = repmat(["solid", "dotted", "dotdashed"]', 6, 1)[:]
    for test_type_df in groupby(df, "test_type")
        test = test_type_df[1, "test_type"]

        for col in cols
            if length(DataFrames.removeNA(df[col])) == 0; continue; end

            seq = split(col, '_')[1]
            seq_name = sort_descr[seq[1]]

            ts = "$(test)_$(seq_name)" |> x -> replace(x, r"[ \[]", "_") |> x -> replace(x, '\]', "")

            println("Plotting $ts")
            plt = FramedPlot()
            setattr(plt, "xlabel", "log n")
            setattr(plt, "ylabel", "time (relative to $base_sort")
            setattr(plt, "title", "Sort Comparison ($test, $seq_name)")
            #setattr(plt.y1, "log", true )
            setattr(plt, "yrange", (0, 1.35))
            setattr(plt, "width", 768)

            curves = Any[]
            for (i, sort_alg_df_all) in enumerate(groupby(test_type_df, "sort_alg"))
                sort_alg_df = sort_alg_df_all[complete_cases(sort_alg_df_all[:,[col]]),:]
                sort_alg = sort_alg_df[1, "sort_alg"]
                c = colors[sort_alg]
                s = linestyles[i]
                pts = Points(sort_alg_df["log_size"], sort_alg_df[col], "color", c, "symboltype", "circle")
                cv = Curve(sort_alg_df["log_size"], sort_alg_df[col], "color", c, "linestyle", s)
                setattr(cv, "label", replace(sort_alg, "_", "\\_"))
                push!(curves, cv)
                add(plt, pts)
                if length(sort_alg_df["log_size"]) > 1
                    add(plt, cv)
                end
            end

            add(plt, Legend(.1, .9, curves))

            push!(plots, plt)
        end
    end

    plots
end

# Save plots
save_sort_plots(plots, pdffile="sortperf.pdf") = file(plots, pdffile)
#save_sort_plots(df::DataFrame, pdffile="sortperf.pdf") =
#    save_sort_plots(sort_plots(sort_median(df), ["*sort_median", "\\sort_median", "/sort_median", "3sort_median",
#                                                 "+sort_median", "~sort_median", "=sort_median", "!sort_median"]),
#                               pdffile)
save_sort_plots(df::DataFrame, base_sort, pdffile="sortperf.pdf") =
    save_sort_plots(sort_plots(sort_scale(sort_median(df), base_sort), base_sort), pdffile)

function view_sort_plots(plots)
    for p in plots
        Winston.display(p)
    end
end

# Test standard sort functions
function std_sort_tests(;sort_algs=SortPerf.sort_algs, types=SortPerf.std_types, range=6:20, replicates=3,
                        lt::Function=isless, by::Function=identity, rev::Bool=false, order::Ordering=Forward,
                        save::Bool=false, prefix="sortperf", skip_median_killer=false, base_sort="QuickSort")
    sort_times = sortperf(sort_algs, types, range, ord(lt,by,rev,order);
                           replicates=replicates, skip_median_killer=skip_median_killer)

    if save
        pdffile = prefix*".pdf"
        tsvfile = prefix*".tsv"
        save_sort_plots(sort_times, base_sort, pdffile)
        writetable(tsvfile, sort_times)
    end

    sort_times
end

end # module SortPerf


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    add_electric_tariff_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")

Adds the ElectricTariff results to the dictionary passed back from `run_reopt` using the solved model `m` and the `REoptInputs` for node `_n`.
Note: the node number is an empty string if evaluating a single `Site`.

ElectricTariff results:
- `lifecycle_energy_cost_after_tax` lifecycle cost of energy from the grid in present value, after tax
- `year_one_energy_cost_before_tax` cost of energy from the grid over the first year, before considering tax benefits
- `lifecycle_demand_cost_after_tax` lifecycle cost of power from the grid in present value, after tax
- `year_one_demand_cost_before_tax` cost of power from the grid over the first year, before considering tax benefits
- `lifecycle_fixed_cost_after_tax` lifecycle fixed cost in present value, after tax
- `year_one_fixed_cost_before_tax` fixed cost over the first year, before considering tax benefits
- `lifecycle_min_charge_adder_after_tax` lifecycle minimum charge in present value, after tax
- `year_one_min_charge_adder_before_tax` minimum charge over the first year, before considering tax benefits
- `year_one_bill_before_tax` sum of `year_one_energy_cost_before_tax`, `year_one_demand_cost_before_tax`, `year_one_fixed_cost_before_tax`, `year_one_min_charge_adder_before_tax`, and `year_one_coincident_peak_cost_before_tax`
- `lifecycle_export_benefit_after_tax` lifecycle export credits in present value, after tax
- `year_one_export_benefit_before_tax` export credits over the first year, before considering tax benefits
- `lifecycle_coincident_peak_cost_after_tax` lifecycle coincident peak charge in present value, after tax
- `year_one_coincident_peak_cost_before_tax` coincident peak charge over the first year
"""
function add_electric_tariff_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
    r = Dict{String, Any}()
    m[Symbol("Year1UtilityEnergy"*_n)] = p.hours_per_time_step *
        sum(m[Symbol("dvGridPurchase"*_n)][ts, tier] for ts in p.time_steps, tier in 1:p.s.electric_tariff.n_energy_tiers)

    r["lifecycle_energy_cost_after_tax"] = round(value(m[Symbol("TotalEnergyChargesUtil"*_n)]) * (1 - p.s.financial.offtaker_tax_pct), digits=2)
    r["year_one_energy_cost_before_tax"] = round(value(m[Symbol("TotalEnergyChargesUtil"*_n)]) / p.pwf_e, digits=2)

    r["lifecycle_demand_cost_after_tax"] = round(value(m[Symbol("TotalDemandCharges"*_n)]) * (1 - p.s.financial.offtaker_tax_pct), digits=2)
    r["year_one_demand_cost_before_tax"] = round(value(m[Symbol("TotalDemandCharges"*_n)]) / p.pwf_e, digits=2)

    r["lifecycle_fixed_cost_after_tax"] = round(m[Symbol("TotalFixedCharges"*_n)] * (1 - p.s.financial.offtaker_tax_pct), digits=2)
    r["year_one_fixed_cost_before_tax"] = round(m[Symbol("TotalFixedCharges"*_n)] / p.pwf_e, digits=0)

    r["lifecycle_min_charge_adder_after_tax"] = round(value(m[Symbol("MinChargeAdder"*_n)]) * (1 - p.s.financial.offtaker_tax_pct), digits=2)
    r["year_one_min_charge_adder_before_tax"] = round(value(m[Symbol("MinChargeAdder"*_n)]) / p.pwf_e, digits=2)

    r["lifecycle_export_benefit_after_tax"] = -1 * round(value(m[Symbol("TotalExportBenefit"*_n)]) * (1 - p.s.financial.offtaker_tax_pct), digits=2)
    r["year_one_export_benefit_before_tax"] = -1 * round(value(m[Symbol("TotalExportBenefit"*_n)]) / p.pwf_e, digits=0)

    r["lifecycle_coincident_peak_cost_after_tax"] = round(value(m[Symbol("TotalCPCharges"*_n)]) * (1 - p.s.financial.offtaker_tax_pct), digits=2)
    r["year_one_coincident_peak_cost_before_tax"] = round(r["lifecycle_coincident_peak_cost_after_tax"] / p.pwf_e, digits=2)

    r["year_one_bill_before_tax"] = r["year_one_energy_cost_before_tax"] + r["year_one_demand_cost_before_tax"] +
                                    r["year_one_fixed_cost_before_tax"]  + r["year_one_min_charge_adder_before_tax"] + r["year_one_coincident_peak_cost_before_tax"]

    d["ElectricTariff"] = r
    nothing
end

function add_electric_tariff_results(m::JuMP.AbstractModel, p::MPCInputs, d::Dict; _n="")
    r = Dict{String, Any}()
    m[Symbol("energy_purchased"*_n)] = p.hours_per_time_step *
        sum(m[Symbol("dvGridPurchase"*_n)][ts] for ts in p.time_steps)

    r["energy_cost"] = round(value(m[Symbol("TotalEnergyChargesUtil"*_n)]), digits=2)

    r["demand_cost"] = round(value(m[Symbol("TotalDemandCharges"*_n)]), digits=2)

    r["export_benefit"] = -1 * round(value(m[Symbol("TotalExportBenefit"*_n)]), digits=0)

    d["ElectricTariff"] = r
    nothing
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
rad(x, y, z) = sqrt(x^2 + y^2 + z^2)
lat(x, y, z) = asin(z / rad(x, y, z)) # ϕ ∈ [-π/2, π/2]
lon(x, y, z) = atan(y, x) # λ ∈ [-π, π)

r̂ⁿᵒʳᵐ(x, y, z) = norm([x, y, z]) ≈ 0 ? 1 : norm([x, y, z])^(-1)
ϕ̂ⁿᵒʳᵐ(x, y, z) =
    norm([x, y, 0]) ≈ 0 ? 1 : (norm([x, y, z]) * norm([x, y, 0]))^(-1)
λ̂ⁿᵒʳᵐ(x, y, z) = norm([x, y, 0]) ≈ 0 ? 1 : norm([x, y, 0])^(-1)

r̂(x, y, z) = r̂ⁿᵒʳᵐ(x, y, z) * @SVector([x, y, z])
ϕ̂(x, y, z) = ϕ̂ⁿᵒʳᵐ(x, y, z) * @SVector [x * z, y * z, -(x^2 + y^2)]
λ̂(x, y, z) = λ̂ⁿᵒʳᵐ(x, y, z) * @SVector [-y, x, 0]

rfunc(p, x...) = ρuʳᵃᵈ(p, lon(x...), lat(x...), rad(x...)) * r̂(x...)
ϕfunc(p, x...) = ρuˡᵃᵗ(p, lon(x...), lat(x...), rad(x...)) * ϕ̂(x...)
λfunc(p, x...) = ρuˡᵒⁿ(p, lon(x...), lat(x...), rad(x...)) * λ̂(x...)

ρu⃗(p, x...) = rfunc(p, x...) + ϕfunc(p, x...) + λfunc(p, x...)

function vector_field_representation(
    u⃗,
    grid::DiscontinuousSpectralElementGrid,
    rep::AbstractRepresentation,
)
    n_ijk, _, n_e = size(u⃗)
    uᴿ = copy(u⃗)
    v⃗ =
        VectorField(data = (u⃗[:, 1, :], u⃗[:, 2, :], u⃗[:, 3, :]), grid = grid)
    for ijk in 1:n_ijk, e in 1:n_e, s in 1:3
        uᴿ[ijk, s, e] = v⃗(rep)[ijk, e, verbose = false][s]
    end
    return uᴿ
end

vector_field_representation(simulation::Simulation, rep) =
    vector_field_representation(simulation.state.ρu, simulation.model.grid, rep)


text:
"""
    SeisGain(d, t; <keyword arguments>)

Gain a group of traces.

# Arguments
* `d::Array{Real,2}`: two dimensional data.

# Keyword arguments
* `dt::Real=0.002`: sampling interval in secs.
* `kind::AbstractString="time"`: if kind="time", gain = t.^a . * exp(-bt);
                         if kind="agc", automatic gain control is applied.
* `param::Vector{Real}=[2.0,0.0]`: if kind="time", param = [a,b];
                                   if kind="agc", param = [agc_gate]
* `norm::Int=0`: `norm=0` no normalization; `norm=1` normalize each trace by
                  amplitude; `norm=2` normalize each trace by rms value/

# Output
* d::Array{Real, 2}`: gained two dimensional data.

# Example
```julia
julia> d, extent = SeisHypEvents();
       dout = SeisGain(d, kind="agc", param=[0.05]);
       SeisPlot([d dout], extent);
```

Credits: Juan I. Sabbione, Aaron Staton, Mauricio D. Sacchi, 2016

"""
function SeisGain{Td<:Real,Tp<:Real
                  }(d::Array{Td,2}; dt::Real=0.004, kind::AbstractString="time",
                    param::Vector{Tp}=[2.0,0.0], norm::Int=0)

    nt = size(d,1)
    nx = size(d[:,:],2)
    dout = zeros(d)

    if kind == "time"   # Geometrical spreading-like gain

        a = param[1]
        b = param[2]
        t = collect(0:1:nt-1)*dt
        tgain = (t.^a).*exp.(b.*t)
        for k = 1:nx
            dout[:,k] = d[:,k].*tgain
        end

    end

    if kind=="agc"   # AGC

        L = floor(Int,round(param[1]/(2dt)))
        h = triang(2*L+1)

        for k = 1:nx
            aux =  d[:,k]
            e = aux.^2
            rms = sqrt(abs.(conv(e,h)[L+1:nt+L]))
            epsi = 1.e-10*maximum(rms)
            op = rms./(rms.^2+epsi)
            dout[:,k] = d[:,k].*op
        end
    end

    if norm==1     # Normalize by amplitude

        for k = 1:nx
            aux =  d[:,k]
            amax = maximum(abs.(aux))
            dout[:,k] = dout[:,k]/amax
        end

    end

    if norm==2;    # Normalize by rms

        for k = 1:nx
            aux = d[:,k];
            amax = sqrt(sum(aux.^2)/nt);
            dout[:,k] = dout[:,k]/amax;
        end

    end

    return dout

end

function triang(n::Integer)
    [1 - abs((k - (n-1)/2))/(n/2) for k=0:(n-1)]
end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
module Hidden
    export is_sum_of_cubes

    function is_cubic(n::Int)::Bool
        n == mapreduce(x -> x^3, +, digits(n))
    end

    function format_return(cubics::Vector{Int})::String
        if length(cubics) > 0
            string( join(map(x->string(x," "),cubics)), sum(cubics), " Lucky" )
        else
            "Unlucky"
        end
    end

    function is_sum_of_cubes(s::AbstractString)::String
        cubics::Vector{Int} = []
        index = 1
        while index <= length(s)
            m = match(r"\d{1,3}", s, index)
            if m !== nothing
                number = parse(Int,m.match)
                if is_cubic(number)
                    push!(cubics,number)
                end
                index = m.offset + length(m.match)
            else
                index = length(s)+1
            end
        end
        format_return(cubics)
    end

end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
###########################################################
#
#   GenericIncidenceList{V, E, VList, IncList}
#
#   V :         vertex type
#   E :         edge type
#   VList :     the type of vertex list
#   IncList:    the incidence list
#
###########################################################

mutable struct GenericIncidenceList{V, E, VList, IncList} <: AbstractGraph{V, E}
    is_directed::Bool
    vertices::VList
    nedges::Int
    inclist::IncList
end

const SimpleIncidenceList = GenericIncidenceList{Int, IEdge, UnitRange{Int}, Vector{Vector{IEdge}}}
const IncidenceList{V,E} = GenericIncidenceList{V, E, Vector{V}, Vector{Vector{E}}}

@graph_implements GenericIncidenceList vertex_list vertex_map edge_map adjacency_list incidence_list

# construction

simple_inclist(nv::Integer; is_directed::Bool=true) =
    SimpleIncidenceList(is_directed, intrange(nv), 0, multivecs(IEdge, nv))

inclist(vs::Vector{V}, ::Type{E}; is_directed::Bool = true) where {V,E} =
    IncidenceList{V,E}(is_directed, vs, 0, multivecs(E, length(vs)))

inclist(::Type{V}, ::Type{E}; is_directed::Bool = true) where {V,E} = inclist(V[], E; is_directed=is_directed)
inclist(vs::Vector{V}; is_directed::Bool = true) where {V} = inclist(vs, Edge{V}; is_directed=is_directed)
inclist(::Type{V}; is_directed::Bool = true) where {V} = inclist(V[], Edge{V}; is_directed=is_directed)

# First constructors on Dict Inc List version (reusing GenericIncidenceList container and functions, few dispatch changes required)
const IncidenceDict{V,E} = GenericIncidenceList{V, E, Dict{Int,V}, Dict{Int,Vector{E}}}
incdict(vs::Dict{Int,V}, ::Type{E}; is_directed::Bool = true) where {V,E} =
    IncidenceDict{V,E}(is_directed, vs, 0, Dict{Int, E}())
incdict(::Type{V}; is_directed::Bool = true) where {V} = incdict(Dict{Int,V}(), Edge{V}; is_directed=is_directed)

# required interfaces

is_directed(g::GenericIncidenceList) = g.is_directed

num_vertices(g::GenericIncidenceList) = length(g.vertices)
# vertices(g::GenericIncidenceList) = g.vertices

# dictionary enables version
vertices_specific(a::UnitRange{Int}) = a
vertices_specific(a::Vector{V}) where {V} = a
vertices_specific(d::Dict{Int,V}) where {V} = collect(values(d))
vertices(g::GenericIncidenceList) = vertices_specific(g.vertices)

num_edges(g::GenericIncidenceList) = g.nedges

edge_index(e::E, g::GenericIncidenceList{V,E}) where {V,E} = edge_index(e)

out_edges(v::V, g::GenericIncidenceList{V}) where {V} = g.inclist[vertex_index(v, g)]
out_degree(v::V, g::GenericIncidenceList{V}) where {V} = length(out_edges(v, g))
out_neighbors(v::V, g::GenericIncidenceList{V}) where {V} = TargetIterator(g, g.inclist[vertex_index(v, g)])

"""
Find neighbors connected by directed edge towards `vert`.
"""
function in_neighbors(vert::V, gr::GenericIncidenceList{V, Edge{V}, Vector{V}}) where {V}
  inclist = gr.inclist
  targid = vert.index
  inlist = V[]
  for edgelist in inclist
    for ed in edgelist
      if ed.target.index == targid
        push!(inlist, ed.source)
      end
    end
  end
  return inlist
end
function in_neighbors(vert::V, gr::GenericIncidenceList{V, Edge{V}, Dict{Int, V}}) where {V}
  inclist = gr.inclist
  targid = vert.index
  inlist = V[]
  for (key,edgelist) in inclist
    for ed in edgelist
      if ed.target.index == targid
        push!(inlist, ed.source)
      end
    end
  end
  return inlist
end

# mutation

function add_vertex!(vertices::Vector{V}, inclist::Vector{E}, v::V) where {V,E}
    push!(vertices, v)
    push!(inclist, Array{E}(undef, 0))
    v
end
function add_vertex!(vertices::Dict{Int, V}, inclist::Dict{Int,E}, v::V) where {V,E}
  if haskey(vertices, v.index)
    error("Already have index $(v.index) in g")
  end
  vertices[v.index] = v
  inclist[v.index] = Array{E}(undef, 0)
  v
end
function add_vertex!(g::GenericIncidenceList{V,E}, v::V) where {V,E}
  add_vertex!(g.vertices, g.inclist, v)
end

function delete_vertex!(vertices::Dict{Int, V}, inclist::Dict{Int, E}, v::V, outnei) where {V,E}
  nedges = 0
  # delete all connected edges
  for vid in union(map(x->x.index, outnei), v.index)
    count = 0
    keeplist = E()
    for ed in inclist[vid]
      count += 1
      # want to delete any source or dest edge pointing to v (not whole inclist)
      if ed.source.index == v.index || ed.target.index == v.index
        # this edge must be deleted from Array
        nedges += 1
      else
        push!(keeplist, ed)
      end
    end
    inclist[vid] = keeplist
  end
  delete!(inclist, v.index)
  delete!(vertices, v.index)
  return nedges
end

function delete_vertex!(v::V, g::GenericIncidenceList{V,E}) where {V,E}
  # find list of vertices that may have edges to v
  possv = collect(out_neighbors(v, g))

  # delete the vertex
  ned = delete_vertex!(g.vertices, g.inclist, v, possv)

  g.nedges -= (is_directed(g) ? ned : round(Int,ned/2))

  return nothing
end

add_vertex!(g::GenericIncidenceList, x) = add_vertex!(g, make_vertex(g, x))

function add_edge!(g::GenericIncidenceList{V,E}, u::V, v::V, e::E) where {V,E}
    # add an edge between (u, v)
    ui::Int = vertex_index(u, g)
    push!(g.inclist[ui], e)
    g.nedges += 1

    if !g.is_directed
        vi::Int = vertex_index(v, g)
        push!(g.inclist[vi], revedge(e))
    end
end

add_edge!(g::GenericIncidenceList{V,E}, e::E) where {V,E} = add_edge!(g, source(e, g), target(e, g), e)
add_edge!(g::GenericIncidenceList{V, E}, u::V, v::V) where {V,E} = add_edge!(g, u, v, make_edge(g, u, v))


text:
hello = "Hello, World!"
println(hello)

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
function read_vec(fpath; builtin = false)
    if builtin
        f = open(fpath, "r")
        lines = readlines(f)
        close(f)
        map(x -> float(x), lines)
    else
        readdlm(fpath)
    end
end

function write_vec(outpath, vec; msg = true)
  f = open(outpath, "w")
  for elem in vec
    write(f, string(elem) * "\n")
  end
  close(f)
  if msg
    println("Write vector into $outpath")
  end
end

function writedlm_header(fpath, matrix, header::AbstractString)
    f = open(fpath, "w")
    write(f, header * "\n")
    writedlm(f, matrix, ",")
    close(f)
end

function join_two_columns(in1, in2, out)
  f1 = open(in1, "r")
  f2 = open(in2, "r")
  f3 = open(out, "w")
  lines1 = readlines(f1)
  lines2 = readlines(f2)
  println(length(lines1) - length(lines2))
  if length(lines1) == length(lines2)
    for i = 1:length(lines1)
        write(f3, string(replace(lines1[i], "\n", "")) * "," *
              string( replace(lines2[i],"\n","")) * "\n")
    end
  end
  close(f1)
  close(f2)
  close(f3)
end
# mpath = "D:/DATASETS/OUT_SAC16Data/bus1000/w1000/jl/"
# join_two_columns(mpath*"best_columns.txt", mpath*"best_sigmas.txt", mpath*"clsgm.txt")

function correct_path(fpath::AbstractString)
  print(fpath[end])
  if fpath[end] != '/'
    fpath = fpath * "/"
  end
  fpath
end

function write_items(fpath, vec, header)
    f=open(fpath,"w")
    write(f, header * "\n")
    for e in vec
        str=join(e, ",")*"\n"
        write(f,str)
    end
    close(f)
end
function write_items(fpath, vec)
    f=open(fpath,"w")
    for e in vec
        str=join(e, ",")*"\n"
        write(f,str)
    end
    close(f)
end
#write_items("D:/ttt.txt", [(1,2,'a'),(5,7,'v')])
#write_items("D:/ttt1.txt", [(1,2,'a'),(5,7,'v')],"dwdw,de,de,d")

function rms(v)
  n = length(v)
  sqrt(sum(v.*v)/n)
end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
struct StepSize1D1 <: StepSize
end

Base.:*(ss::StepSize1D1, k::Float64) = ss

struct StepSize1Dj <: StepSize
    dq::Float64
end

Base.:*(ss::StepSize1Dj, k::Float64) = StepSize1Dj(ss.dq * k)

struct State1D{P} <: State{P}
    qs::Vector{Float64}
end

State1D(P::Int, xi_init::Float64) = State1D{P}(xi_init * ones(P))

Base.copy(st::State1D{P}) where {P} = State1D{P}(copy(st.qs))

path_xi(st::State1D) = st.qs

function density_diff(sys::PimcSystem, st1::State1D{P}, st2::State1D{P}, j::Int) where{P}
    pre_density = (sys.pot(st2.qs[j]) - sys.pot(st1.qs[j])) / P
    k = sys.mass * P / (2 * hbar^2 * sys.beta^2)

    for j_p in [bead_prev(P, j), bead_next(P, j)]
        link1 = (st1.qs[j] - st1.qs[j_p])^2
        link2 = (st2.qs[j] - st2.qs[j_p])^2
        pre_density += k * (link2 - link1)
    end

    pre_density
end

step_1(::PimcSystem, st::State1D, ss::StepSize1D1) = true, st

function step_j(sys::PimcSystem, st::State1D, j::Int, ss::StepSize1Dj)
    st_new = copy(st)
    st_new.qs[j] += ss.dq * randn()

    pre_density = density_diff(sys, st, st_new, j)

    if rand() < exp(-sys.beta * pre_density)
        true, st_new
    else
        false, st
    end
end

function deriv_estimator_pre(::State1D)
    dlogJ = 0.0
    chvar = [1.0]

    dlogJ, chvar
end

function deriv_estimator_1(sys::PimcSystem, chvar::Vector{Float64}, st::State1D, j::Int)
    dot(chvar, sys.force(st.qs[j]))
end

function deriv_estimator_2(chvar::Vector{Float64}, st::State1D{P}) where {P}
    dot(chvar, 2 .* st.qs[1] .- st.qs[2] .- st.qs[P])
end


text:
# TODO: make it type-stable?
# type parameters can be (nested) Symbol, Tuple, Function, and NametedTuple
# so we can embed all hooks and req prototypes into the type parameter
struct RouteNode
    name::String
    hook::Dict{String, Vector{Function}}
    children::Vector{RouteNode}
end

RouteNode(name) = RouteNode(name, Dict(x => Function[] for x in [map(uppercase∘string, methods); "HOOK"]), RouteNode[])

function find_node(node, seg, i=1)::RouteNode
    i > length(seg) && return node
    token = seg[i]
    child = findfirst(x->startswith(x.name, ':') || x.name == token, node.children)
    if child == nothing
        push!(node.children, RouteNode(token))
        child = length(node.children)
    end
    find_node(node.children[child], seg, i+1)
end

function (node::RouteNode)(req, res)
    seg, route = HTTP.URIs.splitpath(req.uri.path), PDict(:node=>node)

    function _continue_route(node, k)
        if k > length(seg)
            method = req.method
            isempty(node.hook[method]) && return res.route_status = 405
            return _continue_handle(node.hook[method], 1)
        end

        token = seg[k]
        child = findfirst(x->startswith(x.name, ':') || x.name == token, node.children)
        if child == nothing
            return res.route_status = 404
        end
        node = node.children[child]
        route.node = node
        startswith(node.name, ':') && setproperty!(route, Symbol(SubString(node.name, 2)), token)
        return _continue_hook(node, node.hook["HOOK"], 1, k+1)
    end

    function _continue_hook(node, hooks, i, k)
        i > length(hooks) && return _continue_route(node, k)
        route.next = () -> _continue_hook(node, hooks, i+1, k)
        hooks[i](req, res, route)
    end

    function _continue_handle(handlers, i)
        i > length(handlers) && return @warn "the last handler calls `next`" path='/'*join(seg, '/')
        route.next = () -> _continue_handle(handlers, i+1)
        handlers[i](req, res, route)
    end

    _continue_hook(node, node.hook["HOOK"], 1, 1)
end

function build_routing_tree(rules)
    root = RouteNode("")
    push!(root.hook["HOOK"], default) # attach the default hook, which is a part of the basic Restful functionality

    for rule in rules
        method, path = rule
        node = find_node(root, HTTP.URIs.splitpath(path))
        push!(node.hook[method], rule[3:end]...)
    end

    root
end


text:
# Test file for linkage.jl

# These tests compare the results with the R interface of fastcluster,
# NOT with another implementation of hierarchical clustering. Hence,
# they should be taken as a validity of the interface, NOT the package
# as such.

using Test
using RDatasets
using Distances
using DataFrames
using RCall

using Fastcluster

# this function assumes that clusters are indexed 1:k
function are_clusters_equal(c1::Vector{T}, c2::Vector{T}) where {T<:Integer}
    k = maximum(c1)
    if k!=maximum(c2)
        # different number of clusters
        return false
    end
    out = zeros(T, k,k)
    for i = 1:k
        for j = 1:k
            out[i,j] = sum((c1 .== i) .& (c2 .== j))
        end
    end
    return all(sum(out.>=1,dims=1) .== 1) & all(sum(out.>=1,dims=2) .== 1)
end

# Start the testing. Happy testing.

df = dataset("datasets", "iris")
#
# answers = Dict{Symbol,Vector{Int64}}()
#
# answers[:single] = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
# answers[:complete] = [1, 2, 2, 2, 1, 1, 2, 1, 2, 2, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 2, 2, 1, 1, 1, 2, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 2, 2, 1, 1, 2, 1, 2, 1, 1, 3, 1, 3, 2, 1, 2, 1, 2, 1, 2, 2, 1, 1, 1, 2, 3, 2, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 2, 1, 3, 1, 2, 2, 2, 1, 2, 2, 2, 2, 2, 1, 2, 2, 1, 2, 3, 1, 1, 3, 2, 3, 1, 3, 1, 1, 1, 2, 2, 1, 1, 3, 3, 1, 3, 2, 3, 1, 3, 3, 1, 1, 1, 3, 3, 3, 1, 1, 1, 3, 1, 1, 1, 3, 3, 3, 2, 3, 3, 1, 1, 1, 1, 1]
# answers[:average] = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 1, 2, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 1, 2, 2, 2, 3, 2, 2, 3, 1, 3, 2, 3, 2, 2, 2, 2, 2, 2, 2, 3, 3, 2, 2, 2, 3, 2, 2, 3, 2, 2, 2, 3, 3, 3, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
# answers[:weighted] = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 2, 1, 2, 1, 2, 1, 1, 2, 2, 2, 1, 2, 1, 1, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 2, 1, 2, 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1, 3, 2, 2, 3, 1, 3, 2, 3, 2, 2, 2, 1, 1, 2, 2, 3, 3, 2, 2, 1, 3, 2, 2, 3, 2, 2, 2, 3, 3, 3, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2]
# answers[:centroid] = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 1, 2, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 1, 2, 2, 2, 3, 2, 2, 3, 1, 3, 2, 3, 2, 2, 2, 2, 2, 2, 2, 3, 3, 2, 2, 2, 3, 2, 2, 3, 2, 2, 2, 3, 3, 3, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
# answers[:median] = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 2, 1, 2, 1, 2, 1, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 1, 2, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1, 3, 2, 2, 3, 1, 3, 2, 3, 2, 2, 2, 1, 1, 2, 2, 3, 3, 2, 2, 1, 3, 2, 2, 3, 2, 1, 2, 3, 3, 3, 2, 2, 1, 3, 2, 2, 1, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 1]
# answers[:ward] = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 2, 3, 2, 1, 2, 1, 1, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 1, 3, 3, 3, 3, 1, 3, 2, 3, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 3, 3, 2, 2, 2, 2, 3, 2, 3, 2, 2, 2, 2, 3, 3, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 3, 2, 2, 2, 3, 2, 2, 2, 3, 2, 2, 3]
#
# for meth in [:single, :complete, :average, :weighted, :centroid, :median]
#     println("Testing interface with method=:$(string(meth))...")
#     for k in [3]
#         # run once with R interface, and once with Fastcluster.jl interface
#         points = convert(Array{Float64,2},df[:,[:SepalWidth, :SepalLength]])
#         d = pairwise(Euclidean(), points, dims=1)
#         m,h = Fastcluster.linkage(d, meth)
#         c_fastcluster = Fastcluster.cutree(m,(length(m)>>1)+1,k)
#         #@show c_fastcluster
#         @test are_clusters_equal(c_fastcluster, answers[meth])
#     end
#     println("Ok.")
# end
#
# # for
# for meth in [:ward]
#     println("Testing interface with method=:$(string(meth))...")
#     for k in [3]
#         # run once with R interface, and once with Fastcluster.jl interface
#         points = convert(Array{Float64,2},df[:,[:SepalWidth, :SepalLength]])
#         d = pairwise(SqEuclidean(), points, dims=1)
#         m,h = Fastcluster.linkage(d, meth)
#         c_fastcluster = Fastcluster.cutree(m,(length(m)>>1)+1,k)
#         d2 = pairwise(Euclidean(), points, dims=1)
#         #@show c_fastcluster
#         @test are_clusters_equal(c_fastcluster, answers[meth])
#     end
#     println("Ok.")
# end

R".libPaths(c(.libPaths(), temp <- tempdir()))
install.packages(\"fastcluster\", lib=temp, repos='http://cran.us.r-project.org')"
#R"install.packages(\"fastcluster\", lib = Sys.getenv(\"JL_PKG\"))"

# function to wrap R interface:
function cluster_and_cut_R(d::Array{T,2}, method::Symbol, nclusters::Int64) where {T<:Real}

    @rput d

    mthidx = Dict(:single => 0,
              :complete  => 1,
              :average   => 2,
              :weighted  => 3,
              :ward      => 4,
              :centroid  => 5,
              :median    => 6 )

    # for R's fastcluster, method must be: "single", "complete",
    # "average", "mcquitty", "ward.D",
    # "ward.D2", "centroid" or "median".

    R"library('fastcluster')"

    if method == :single
        R"clusters <- hclust(as.dist(d), \"single\")"
    elseif method == :complete
        R"clusters <- hclust(as.dist(d), \"complete\")"
    elseif method == :average
        R"clusters <- hclust(as.dist(d), \"average\")"
    elseif method == :weighted
        R"clusters <- hclust(as.dist(d), \"mcquitty\")"
    elseif method == :ward
        R"clusters <- hclust(as.dist(d), \"ward.D2\")"
    elseif method == :ward1
        R"clusters <- hclust(as.dist(d), \"ward.D\")"
    elseif method == :centroid
        R"clusters <- hclust(as.dist(d), \"centroid\")"
    elseif method == :median
        R"clusters <- hclust(as.dist(d), \"median\")"
    else
        error("Invalid or unsupported method.")
    end

    @rput nclusters
    R"clusterCut <- cutree(clusters, k = nclusters)"
    @rget clusterCut

    return clusterCut

end

for meth in [:single, :complete, :average, :weighted, :centroid, :median]
    println("Testing interface with method=:$(string(meth))...")
    for k in [1 , 3, 5]
        # run once with R interface, and once with Fastcluster.jl interface
        points = convert(Array{Float64,2},df[:,[:SepalWidth, :SepalLength]])
        d = pairwise(Euclidean(), points, dims=1)
        m,h = Fastcluster.linkage(d, meth)
        c_fastcluster = Fastcluster.cutree(m,(length(m)>>1)+1,k)
        r_fastcluster = cluster_and_cut_R(d, meth, k)
        @test are_clusters_equal(c_fastcluster, r_fastcluster)
    end
    println("Ok.")
end

# for
for meth in [:ward]
    println("Testing interface with method=:$(string(meth))...")
    for k in [1 , 3, 5]
        # run once with R interface, and once with Fastcluster.jl interface
        points = convert(Array{Float64,2},df[:,[:SepalWidth, :SepalLength]])
        d = pairwise(SqEuclidean(), points, dims=1)
        m,h = Fastcluster.linkage(d, meth)
        c_fastcluster = Fastcluster.cutree(m,(length(m)>>1)+1,k)
        d2 = pairwise(Euclidean(), points, dims=1)
        r_fastcluster = cluster_and_cut_R(d2, meth, k)
        @test are_clusters_equal(c_fastcluster, r_fastcluster)
    end
    println("Ok.")
end


text:
# to run from the nearly-neutral/src folder:
# julia -L InfAlleles.jl run.jl ../../NearlyNeutralPaper/results21/Richness_figure3/ngens200K/nn_neut_Nmulist_fm45_N25_6400
#=
Estimates of theta = N_e mu based on Watterson "The homozygosity test of neutrality" (1978).

allele_freqs is the vector of allele frequecies of the infinite alleles model population.
=#

export watterson_homozygosity, p_homozygosity, watterson_theta, IQV, w_homoz

function watterson_homozygosity( allele_freqs::Population )
  if length(allele_freqs) == 0
    return 0.0   # Not sure that this is valid
  end
  n = sum(allele_freqs)
  sum_a_sq = 0
  for a in allele_freqs
    sum_a_sq += a^2
  end
  Float64(sum_a_sq)/Float64(n^2)
end
#=
function watterson_theta( allele_freqs::Config )
  1.0/watterson_homozygosity( allele_freqs ) - 1.0
end
=#

function watterson_theta( allele_freqs::Population )
  1.0/watterson_homozygosity( allele_freqs ) - 1.0
end

@doc """ function IQV()
Allan Wilcox's "index of quantitative variation" defined by formula (7)
  of Mark E. Madsen's "Neutral Cultural Transmission in Time Averaged
Archaeological Assemblages" arXiv:1024.2043v2 (2012).
A measure of "evenness".  IQV of a population with a single trait is 0,
while IQV of a population with an even distribution of traits is close to 1.0.
"""
function IQV( allele_freqs::Population )
  k = length(allele_freqs)
  if k == 1
    return 0.0
  end
  Fk = Float64(k)
  wh = watterson_homozygosity( allele_freqs )
  (1.0-wh)*Fk/(Fk-1.0)
end

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
#!/usr/bin/julia
using AssociationsGUI
folder = joinpath(Pkg.dir("AssociationsGUI"), "test", "videofolder")
main(folder)