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

using Markdown
using InteractiveUtils

# ╔═╡ 61dc86f0-238f-11eb-054f-41a45b78e94d
begin
	using  MLJ
	using  MLJLinearModels
	using  CSV
	import RDatasets: dataset
    import DataFrames: describe, select, Not, rename!,head,DataFrame, nename!
	import Plots: plot,plot!,scatter,scatter!,histogram,histogram!
end

# ╔═╡ ff5bd480-29c0-11eb-3112-cb35c17acc60
md" ### Libs"

# ╔═╡ a2efe1c0-29c8-11eb-287d-49b90e597876
md"### Sobre o Problema

url: https://archive.ics.uci.edu/ml/datasets/banknote+authentication

Data Set Information:

Data were extracted from images that were taken from genuine and forged banknote-like specimens. For digitization, an industrial camera usually used for print inspection was used. The final images have 400x 400 pixels. Due to the object lens and distance to the investigated object gray-scale pictures with a resolution of about 660 dpi were gained. Wavelet Transform tool were used to extract features from images.

Attribute Information:

1. variance of Wavelet Transformed image (continuous)
2. skewness of Wavelet Transformed image (continuous)
3. curtosis of Wavelet Transformed image (continuous)
4. entropy of image (continuous)
5. class (integer)
"

# ╔═╡ 8e226310-29c0-11eb-1093-f95e53c029b6
md"### Carga de Dados"

# ╔═╡ 114c91a0-2a90-11eb-365b-4d360bb51ebe
md" É interessante sempre mostrarmos os primeiros exemplos dos dados para ter uma ideia dos valores. Entretanto, como visto no curso, há várias análises possíveis para se fazer como identificar nulos, outliers, ... Tudo isto, é bom ser tratado."

# ╔═╡ 1c623cf0-2d17-11eb-2f20-ad2b49b221de
begin
	url = "https://archive.ics.uci.edu/ml/machine-learning-databases/00267/data_banknote_authentication.txt"
	tmp = download(url)
end

# ╔═╡ 049db9d0-2d19-11eb-24f1-0ba111942cb8
begin
	banknote = DataFrame(CSV.File(tmp; header=false))
	rename!(banknote, Dict(:Column1 => "Variance", :Column2 => "Skewness",:Column3 => "Curtosis",:Column4 => "Entropy",:Column5 => "Target"))
end

# ╔═╡ 9a810b20-29c0-11eb-1ca3-d155108c6e52
md"### Estatísticas dos Dados"

# ╔═╡ 91187080-2a8f-11eb-3001-1973d885e3f0
md"É sempre bom tiramos umas estatísticas dos dados para entender melhor o comportamento das variáveis do nosso problema. Qual a variância de cada variável? Elas estão na mesma ordem de grandeza? Máixmo e minímo ?"

# ╔═╡ 2033e842-294f-11eb-2b38-1969d7440b88
describe(banknote, :mean, :std, :min, :max, :eltype)

# ╔═╡ 2f06c4a0-294f-11eb-2cb1-2da452f98b0f
begin
	data = coerce(banknote, autotype(banknote, :discrete_to_continuous));
	data = coerce(data,:Target => OrderedFactor)#Multiclass)
end

# ╔═╡ d2a9dfe0-29c0-11eb-36a9-e7e3c1d17717
md"### Selecionando Dados"

# ╔═╡ 49b99b50-2a90-11eb-2b54-2704b914334d
md"É nesta fase que você vai escolher seu target e vai selecionar os atributos para trabalhar."

# ╔═╡ 47a80140-294f-11eb-2880-ef7cfb0a8997
begin
	y = data.Target;
    X = select(data, Not(:Target));
end

# ╔═╡ 429fc960-2d25-11eb-3c06-614daa2fc78f
schema(X)

# ╔═╡ e3d83ea0-29a3-11eb-10f5-b5f94102e1f1
md"### Verificando Algoritmos Disponíveis"

# ╔═╡ 60083560-2a90-11eb-373f-3b8a1ba538a0
md"Essa é uma característica bem legal do MLJ. Como nosso dados possuem atributos numéricos e o target numérico, ele me mostra todos os algoritmos que conseguem rodar para este tipo de dado."

# ╔═╡ 868b73b0-29a9-11eb-0ef8-971446a3c919
#models(matching(X, y))
begin
	task(model) = matching(model, X, y) && model.prediction_type == :probabilistic
	models(task)
end

# ╔═╡ 87cedb30-2a90-11eb-26dd-6162d49edd88
md"Abaixo realizamos uma carga no nosso regressor linear."

# ╔═╡ 9abff260-2a90-11eb-026f-59a85ec6b64c
md"O MLJ nos permite montar um pipeline. Pipeline é uma forma de montar uma sequência de operações poupando-lhe código. No nosso caso, queremos padronizar X e y antes de realizarmos o treino do regressor. Podemos fazer isso em uma única linha."

# ╔═╡ fcde2920-29a6-11eb-37f7-4befad3ee3d2
predictor = @pipeline  Standardizer LogisticClassifier(lambda = 0.5)

# ╔═╡ 323d2620-29c0-11eb-258a-4703b4f5a0d2
md"### Partição dos Dados"

# ╔═╡ e57688a2-2a90-11eb-2272-3d9898bf4375
md"Realizamos aqui o que conhecemos como hold-out. Particionamos em treino e teste, que é a forma mais simples de se fazer uma validação. Sabemos que existem outras no curso. No caso, deixei 70% para treino e 30% para teste. Como podem notar, a função gera índices qu correspondem às linhas dos exemplos. Ele não tira cópia dos dados."

# ╔═╡ fca39122-29a6-11eb-3b25-c7f0df7fc68f
train, test  = partition(eachindex(y), 0.7, shuffle=true);

# ╔═╡ 287ae7e0-2a91-11eb-0d5b-b138443bdc29
md"Um conceito interessante do MLJ é o machine. Esta função faz uma mapeamento do preditor ao dado. Quando aplicamos machine, estamos gerando um nó num grafo de computação e um estado associado. Existem algmas vantagens neste paradigma como: associar um mesmo preditor a vários dados distintos que podem ser processados em diferentes recursos; evitar re execução de um mesmo preditor e um mesmo dado."

# ╔═╡ fc474050-29a6-11eb-1fbd-87a88d38525a
pred_machine = machine(predictor, X, y);

# ╔═╡ 4449f8d0-29b5-11eb-3f83-c352d034ecb2
md"### Treino do Classificador"

# ╔═╡ 91918630-2a91-11eb-1113-43b18086f208
md" Aqui temos enfim, o treino do nosso regressor para o data X, apenas no conjunto de treino. Isto porque estamos passando os índices do exemplos."

# ╔═╡ 03caa2a0-2d2a-11eb-0742-413f150b9056
fit!(pred_machine, rows=train);

# ╔═╡ 50f36580-29b5-11eb-06f2-cf4b9e142979
md"### Análise da Predição"

# ╔═╡ ef445730-2a91-11eb-1552-c3e7d6e5ff34
md"Aqui analisamos cada predição e cada target. Este tipo de análise é mais apropriada quando nosso conjunto de dados possui uma relação de ordem temporal entre os exemplos (que não é o caso). Entretanto, é uma análise boa a se fazer. Podemos notar que há um atendência de subestimar os valores mais altos do que valores mais baixos. Você pode analisar o que pode estar gerando dificldade nisto analisando alguns exemplos individualmente."

# ╔═╡ 36bd3310-2dbd-11eb-0e46-bf6c0b1e80a9
md"Aqui pegamos as classes e probabilidades"

# ╔═╡ 152f7100-29a7-11eb-2300-735ca72b15c3
preds = predict(pred_machine, rows=test)

# ╔═╡ 275bab42-2dbd-11eb-1bd7-a1c0cb2965a7
md"Pegando probabilidades da predição por classe"

# ╔═╡ 1336f910-2d9c-11eb-222b-ade11f5b7b51
begin
	probs_fraude  = broadcast(pdf,preds,1.0);
	probs_nfraude = broadcast(pdf,preds,0.0);
end

# ╔═╡ d4536cb0-2d9b-11eb-16f7-993d0065c683
#mode.(ŷ)
ŷ = predict_mode(pred_machine,rows = test)

# ╔═╡ 35b079a0-2d9f-11eb-0fa6-1bb075367008
measures(matching(y))

# ╔═╡ dbb70eb2-2dc0-11eb-3da4-754b2e0de392
md"#### Métricas para o limiar padrão (0.5)"

# ╔═╡ 32c66e80-2d9e-11eb-220e-093e86fdc72f
cm = confusion_matrix(ŷ,y[test])

# ╔═╡ 3593e090-2e05-11eb-0557-13c977ff89b2
false_positive_rate(cm),false_negative_rate(cm)

# ╔═╡ 463a4470-2e05-11eb-0cac-2f56d62bfa20
true_positive_rate(cm),true_negative_rate(cm) # recall

# ╔═╡ 1fedd0a0-2e07-11eb-1f26-8575ef999663
positive_predictive_value(cm),negative_predictive_value(cm) # precision

# ╔═╡ 6a6e0550-2d9e-11eb-1136-119a2997dcae
accuracy(cm)

# ╔═╡ 0af080a0-2e05-11eb-0069-8d908435f94b
balanced_accuracy(ŷ,y[test])

# ╔═╡ c8288560-2dc8-11eb-14f6-19ae67fb2357
begin
	f1 = FScore{1.0}()
	f1(cm)
end

# ╔═╡ 00989c50-2dc9-11eb-31fc-0d8ab39b5a92
matthews_correlation(cm)

# ╔═╡ 77646ed0-2dc5-11eb-35c5-4b1582878184
md"#### Métricas para limiares distintos"

# ╔═╡ f3292f10-2d9d-11eb-1e76-8de80adce2e6
tprs, fprs, ts = roc_curve(preds,y[test]);

# ╔═╡ 96c67790-2e07-11eb-0000-19536e55469e
ts[180]

# ╔═╡ 903e2090-2dc5-11eb-2855-bfc16c6145b6
md" Área embaixo da curva"

# ╔═╡ 58f2c540-2d9e-11eb-3a14-0931bcd2cc5d
auc(preds,y[test])

# ╔═╡ c29c1380-29b5-11eb-254c-2d026e15f903
plot(tprs,fprs,label="LG Predictions",ylabel="Verdadeiros Positivos",xlabel="Falsos Positivos")

# ╔═╡ 89d06380-29bf-11eb-0d5b-4533c94745d5
md"#### Distribuição das probabilidades das predições"

# ╔═╡ 8b3ac8f0-29bf-11eb-274e-f7b07eb35e7c
histogram(probs_fraude,bins=10,label="Fraude")

# ╔═╡ 1e35c150-29de-11eb-0069-cb75699b6d61
histogram(probs_nfraude,bins=10,label="Não Fraude")

# ╔═╡ 76b4c570-29b5-11eb-32d5-75da0001a894
md"### Parâmetros encontrados"

# ╔═╡ 46a9fc40-2a8e-11eb-301f-077f9a6d50fe
md" Muitas vezes é interessante analisarmos a relevância dos coeficientes por váriosmotivos. Isto porque possibilita descobrir quais são as variáveis que ajudam mais num problema, ajudando inclusive na tomada de decisão de um ser humano. Adicionalmente, pode ser usado para compactar o modelo (ambientes com poucos recursos), eliminando atributos pouco relevantes. Finalmente, pode-se utilizar para verificar se seu modelo usa alguma variável que seja discrimnatória em termos de sociedade (ex. idade para tomada de crédito). "

# ╔═╡ 8d50cc90-294f-11eb-012b-3b973cbd1f1e
begin
	fp = fitted_params(pred_machine)
	coefs = fp.fitted_params_given_machine[fp.machines[2]]
	#fp.machines
end

# ╔═╡ e102cfc0-295c-11eb-2c5b-834000ca2077
#begin
#	res = ŷ .- y[test]
#	plot(res)
#	scatter!(res)
#end

# ╔═╡ 08b31090-2e01-11eb-27a7-c3758492386e
begin
	names = [string(e[1]) for e in coefs[2] ]
	vals  = [e[2] for e in coefs[2] ]
	plot(names,vals,linetype = :bar,label="Coefs")
end

# ╔═╡ Cell order:
# ╟─ff5bd480-29c0-11eb-3112-cb35c17acc60
# ╠═61dc86f0-238f-11eb-054f-41a45b78e94d
# ╟─a2efe1c0-29c8-11eb-287d-49b90e597876
# ╟─8e226310-29c0-11eb-1093-f95e53c029b6
# ╟─114c91a0-2a90-11eb-365b-4d360bb51ebe
# ╠═1c623cf0-2d17-11eb-2f20-ad2b49b221de
# ╠═049db9d0-2d19-11eb-24f1-0ba111942cb8
# ╟─9a810b20-29c0-11eb-1ca3-d155108c6e52
# ╟─91187080-2a8f-11eb-3001-1973d885e3f0
# ╠═2033e842-294f-11eb-2b38-1969d7440b88
# ╠═2f06c4a0-294f-11eb-2cb1-2da452f98b0f
# ╟─d2a9dfe0-29c0-11eb-36a9-e7e3c1d17717
# ╟─49b99b50-2a90-11eb-2b54-2704b914334d
# ╠═47a80140-294f-11eb-2880-ef7cfb0a8997
# ╠═429fc960-2d25-11eb-3c06-614daa2fc78f
# ╟─e3d83ea0-29a3-11eb-10f5-b5f94102e1f1
# ╟─60083560-2a90-11eb-373f-3b8a1ba538a0
# ╠═868b73b0-29a9-11eb-0ef8-971446a3c919
# ╟─87cedb30-2a90-11eb-26dd-6162d49edd88
# ╟─9abff260-2a90-11eb-026f-59a85ec6b64c
# ╠═fcde2920-29a6-11eb-37f7-4befad3ee3d2
# ╟─323d2620-29c0-11eb-258a-4703b4f5a0d2
# ╟─e57688a2-2a90-11eb-2272-3d9898bf4375
# ╠═fca39122-29a6-11eb-3b25-c7f0df7fc68f
# ╠═287ae7e0-2a91-11eb-0d5b-b138443bdc29
# ╠═fc474050-29a6-11eb-1fbd-87a88d38525a
# ╟─4449f8d0-29b5-11eb-3f83-c352d034ecb2
# ╟─91918630-2a91-11eb-1113-43b18086f208
# ╠═03caa2a0-2d2a-11eb-0742-413f150b9056
# ╟─50f36580-29b5-11eb-06f2-cf4b9e142979
# ╟─ef445730-2a91-11eb-1552-c3e7d6e5ff34
# ╟─36bd3310-2dbd-11eb-0e46-bf6c0b1e80a9
# ╠═152f7100-29a7-11eb-2300-735ca72b15c3
# ╟─275bab42-2dbd-11eb-1bd7-a1c0cb2965a7
# ╠═1336f910-2d9c-11eb-222b-ade11f5b7b51
# ╠═d4536cb0-2d9b-11eb-16f7-993d0065c683
# ╠═35b079a0-2d9f-11eb-0fa6-1bb075367008
# ╟─dbb70eb2-2dc0-11eb-3da4-754b2e0de392
# ╠═32c66e80-2d9e-11eb-220e-093e86fdc72f
# ╠═3593e090-2e05-11eb-0557-13c977ff89b2
# ╠═463a4470-2e05-11eb-0cac-2f56d62bfa20
# ╠═1fedd0a0-2e07-11eb-1f26-8575ef999663
# ╠═6a6e0550-2d9e-11eb-1136-119a2997dcae
# ╠═0af080a0-2e05-11eb-0069-8d908435f94b
# ╠═c8288560-2dc8-11eb-14f6-19ae67fb2357
# ╠═00989c50-2dc9-11eb-31fc-0d8ab39b5a92
# ╠═77646ed0-2dc5-11eb-35c5-4b1582878184
# ╠═f3292f10-2d9d-11eb-1e76-8de80adce2e6
# ╠═96c67790-2e07-11eb-0000-19536e55469e
# ╟─903e2090-2dc5-11eb-2855-bfc16c6145b6
# ╠═58f2c540-2d9e-11eb-3a14-0931bcd2cc5d
# ╟─c29c1380-29b5-11eb-254c-2d026e15f903
# ╟─89d06380-29bf-11eb-0d5b-4533c94745d5
# ╟─8b3ac8f0-29bf-11eb-274e-f7b07eb35e7c
# ╟─1e35c150-29de-11eb-0069-cb75699b6d61
# ╟─76b4c570-29b5-11eb-32d5-75da0001a894
# ╟─46a9fc40-2a8e-11eb-301f-077f9a6d50fe
# ╠═8d50cc90-294f-11eb-012b-3b973cbd1f1e
# ╟─e102cfc0-295c-11eb-2c5b-834000ca2077
# ╟─08b31090-2e01-11eb-27a7-c3758492386e

max_stars_repo_name: lalvim/machinelearning

max_issues_repo_name: lalvim/machinelearning

max_forks_repo_name: lalvim/machinelearning


text:
export Resource, InstrumentationInfo

"""
    Resource(;attributes=StaticAttrs(), schema_url="")

Quoted from [the specification](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/overview.md#resources):

> Resource captures information about the entity for which telemetry is recorded. For example, metrics exposed by a Kubernetes container can be linked to a resource that specifies the cluster, namespace, pod, and container name.
>
> Resource may capture an entire hierarchy of entity identification. It may describe the host in the cloud and specific container or an application running in the process.
"""
Base.@kwdef struct Resource{A<:StaticAttrs}
    attributes::A = StaticAttrs()
    schema_url::String = ""
end

#####

"""
    InstrumentationInfo(;name="Main", version=v"0.0.1-dev")

Usually used in an instrumentation package.
"""
Base.@kwdef struct InstrumentationInfo
    name::String = "Main"
    version::VersionNumber = v"0.0.1-dev"
    schema_url::String = ""
end

ext: jl


text:
include("swarm.jl")

function PsoAlgorithm(swarm, targetValue, maxIteration)
    globalBestParticle = swarm[1];
    globalBestParticle = calculateGlobalBestParticle(swarm,targetValue,globalBestParticle);
    for i in 1:maxIteration
        if globalBestParticle.currentValue == targetValue
            printSwarm(swarm)
            println("SOLUTION : ");
            printParticle(globalBestParticle);
            println("Solution found after $i iteration");
            return

        end
        updateSwarmVelocity!(swarm,globalBestParticle,targetValue);
        updateSwarmPosition!(swarm,targetValue);
        globalBestParticle = calculateGlobalBestParticle(swarm,targetValue,globalBestParticle);
    end
    printSwarm(swarm);
    print("GLOBAL SOLUTION : ");
    printParticle(globalBestParticle);
end

function calculateGlobalBestParticle(swarm,targetValue, globalBestParticle)
    for i in 1:length(swarm)
        if abs(swarm[i].currentValue - targetValue) < abs(globalBestParticle.currentValue - targetValue)
            globalBestParticle = deepcopy(swarm[i]);
        end
    end
    return globalBestParticle;
end

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
#currently, matrix solve operations are on hold.

@doc """
  SigmoidNumbers.get_unscaled_replacement_row!(rr, M, row, cache, quire)

  takes the matrix M and specifies a good "replacement row" for it.  You should
  also supply a cache vector which is used to store magic values, it should be
  the length of rr.
"""
function get_unscaled_replacement_row!(rr::Vector{T}, M::Matrix{T}, row, cached_coefficients::Vector{T}, quire) where {T}
  l = length(rr)

  @assert(size(M) == (l, l))
  @assert(length(cached_coefficients) == l)
  @assert(row > 1)

  #generate a quire for use in the function, and a variable to hold intermediates.
  res = zero(T)

  #first column is special.  simply set the magic value cache to be the negative
  #of the row number, then zero it out.
  cached_coefficients[1] = - M[row, 1]
  rr[1] = zero(T)

  for column = (2:row - 1)
    set!(quire, M[row, column]) #we'll subtract values from this as a tally.

    #scan down the vrows to obtain the coefficients for value reduction.
    for vrow = 1:(column - 1)
      res = fdp!(quire, cached_coefficients[vrow], M[vrow, column])
    end

    #finally, set the early column in the row to be zero.
    cached_coefficients[column] = -res
    rr[column]    = zero(T)
  end

  for column = row:l
    set!(quire, M[row, column])
    #scan down the vrows to obtain reduce the coeffients again.
    for vrow = 1:(row -1)
      res = fdp!(quire, cached_coefficients[vrow], M[vrow, column])
    end
    #set the column in the replacement row.
    rr[column] = res
  end

  #don't bother returning anything, since otherwise it would be confusing.
  nothing
end

function diminish_result_vector!(v::Vector{T}, row, cached_coefficients::Vector{T}, quire) where {T}
  #alter the vector in the same fashion.
  set!(quire, v[row])
  res = zero(T)
  for trow = 1:(row - 1)
    res = fdp!(quire, v[trow], cached_coefficients[trow])
  end
  v[row] = res
end

function rescale_row!(M::Matrix{T}, row) where {T}
  value = M[row, row]
  for col = (row + 1):size(M,2)
    M[row,col] /= value
  end
  M[row, row] = one(T)
  value
end

@doc """
  SigmoidNumbers.row_echelon!(M, v)

  converts (M, v) matrix/vector into row-echelon form, with partial pivoting
  Uses quires.
"""
function row_echelon!(M::Matrix{T}, v::Vector{T}, quire = Quire(T)) where {T}
  eqs = length(v)
  @assert size(M) == (eqs, eqs)
  #convert to row-echelon form.

  #create a temporary array that stores intermediate row values.
  t_row = zeros(T, eqs)
  #create a temporary array that stores the prime row values.
  cached_coefficients = zeros(T, eqs)
  cached_replacement_row = zeros(T, eqs)
  #a usable result variable for obtaining a quire result.
  res = zero(T)

  for row = 1:eqs           #go down each row.
    #partial pivoting.
    if row < eqs
      #find the maximum entry among things in this column.
      swap_index = indmax(abs.(M[row:eqs,row])) + row - 1
      if swap_index != row
        cached_replacement_row = M[row, :]
        M[row,:] = M[swap_index, :]
        M[swap_index, :] = cached_replacement_row

        cache = v[row]
        v[row] = v[swap_index]
        v[swap_index] = cache
      end
    end

    if row > 1
      get_unscaled_replacement_row!(cached_replacement_row, M, row, cached_coefficients, quire)

      #replace the row of the matrix.
      M[row,:] = cached_replacement_row

      diminish_result_vector!(v, row, cached_coefficients, quire)
    end

    scale = rescale_row!(M, row)
    v[row] /= scale

  end

  (M, v)
end

@doc """
  SigmoidNumbers.solve_row_echelon!(result, M, v, quire)

  solves a system of linear equations that's in row-echelon form.
"""
function solve_row_echelon!(result::Vector{T}, M::Matrix{T}, v::Vector{T}, quire = Quire(T)) where {T}
  l = length(v)
  @assert(length(result) == l)
  @assert(size(M) == (l,l))
  res = zero(T)
  result[l] = v[l]
  for row = l-1:-1:1
    set!(quire, v[row])
    for col = l:-1:(row+1)
      res = fdp!(quire, M[row, col], -result[col])
    end
    result[row] = res
  end
end

function solve(M::Matrix{T}, v::Vector{T}) where {T}
  quire = Quire(T)
  (Mre, vre) = row_echelon!(copy(M), copy(v), quire)
  r = Vector{T}(length(v))
  solve_row_echelon!(r, Mre, vre, quire)
  r
end

function solve_with_refine(M::Matrix{T}, v::Vector{T}) where {T}
    r = M \ v
    r1 = refine(r, M, v)
    r2 = refine(r, M, v)
end

@doc """
  find_residual(M, r, v)
  calculates v - M * r, using exact dot products.
"""
function find_residuals(M::Matrix{T}, r::Vector{T}, v::Vector{T}, quire = Quire(T)) where {T}
  res = zeros(T, length(v))
  for dim = 1:length(v)
    set!(quire, v[dim])
    for (m, ri) in zip(M[dim, :], r)
      res[dim] = fdp!(quire, m, -ri)
    end
  end
  res
end

@doc """
  refine(r, M, v)

  calculates M * r and finds the residuals with respect to the solution v,
  and then solves those residuals and applies it back to r.
"""
function refine(r::Vector{T}, M::Matrix{T}, v::Vector{T}, quire = Quire(T)) where {T}
  #first, find the residual.
  residual = find_residuals(M, r, v, quire)
  #then solve the residual, and augment r.
  deltas = solve(M, residual)
  r + deltas
end

################################################################################
#for testing purposes, random_exact_row(n) and random_exact_matrix(n) routines.

function obliterate_lsb(pvalue::Posit{N,ES}) where {N,ES}
  (s, exp, frc) = posit_components(pvalue)
  if frc == zero(UInt64)
    #bump up an exponent.
    Posit{N,ES}(s, exp + 1, frc)
  else
    #find the trailing zero, xor it out.
    Posit{N,ES}(s, exp, frc ⊻ (one(UInt64) << trailing_zeros(frc)))
  end
end

@doc """
  SigmoidNumbers.random_exact_row(n)
  creates a row that is summable to an exact number.
"""
function random_exact_row(::Type{Posit{N,ES}}, count, randomizer::Function = rand, q::Quire = Quire(Posit{N,ES})) where {N,ES}
  #first generate a <count> number of random values.
  current_row_vector = Posit{N,ES}.(randomizer(count))

  #also write a function that gets updates the lsb.
  function get_row_lsb()
    zero!(q)
    exact_sum(current_row_vector, q)
    find_lsb(q)
  end

  old_lsb = get_row_lsb()

  while true

    (old_lsb == find_lsb(Posit{N,ES}(q))) && return current_row_vector

    #next twiddle the values
    (_, idx) = findmin(map(find_lsb, current_row_vector))

    #obliterate the lsb in the indexed value
    current_row_vector[idx] = obliterate_lsb(current_row_vector[idx])

    new_lsb = get_row_lsb()

    #eject and redo in the case that our bitflipping is not trivially soluble.  The
    #less than should not strictly be necessary but it's there to prevent
    (new_lsb <= old_lsb) && return random_exact_row(Posit{N,ES}, count, randomizer, q)

    old_lsb = new_lsb
  end
end

function random_exact_matrix(::Type{Posit{N,ES}}, count, randomizer::Function = rand, q::Quire = Quire(Posit{N,ES})) where {N,ES}
  M = zeros(Posit{N,ES}, count, count)
  for row = 1:count
    M[row, :] = random_exact_row(Posit{N,ES}, count, randomizer)
  end
  M
end

exact_rowsum(M::Matrix{Posit{N,ES}}) where {N,ES} = [exact_sum(M[row, :]) for row = 1:size(M,1)]

export random_exact_matrix, exact_rowsum

Base.big(x::Posit{N,ES}) where {N,ES} = big(Float64(x))

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
@recipe function f(model::StateSpaceModel, forec::Forecast)
    if !isunivariate(model)
        error("This plot recipe currently works for univariate models only.")
    end
    dists = Vector{Distribution}(undef, length(forec.expected_value))
    for d in eachindex(dists)
        dists[d] = Normal(forec.expected_value[d][1], sqrt(forec.covariance[d][1]))
    end
    n = length(model.system.y)
    forec_idx = collect(n + 1: n + length(forec.expected_value))
    fillrange_color = :steelblue
    # Plot the series
    @series begin
        seriestype := :path
        seriescolor := "black"
        return model.system.y
    end
    # Plot the forecast
    @series begin
        seriescolor := "blue"
        return forec_idx, mean.(dists)
    end
    # Plot the prediction interval
    @series begin
        seriestype := :path
        seriescolor := fillrange_color
        fillcolor := fillrange_color
        fillalpha := 0.5
        fillrange := quantile.(dists, 0.975)
        label := ""
        return forec_idx, quantile.(dists, 0.025)
    end
    @series begin
        seriestype := :path
        seriescolor := fillrange_color
        fillcolor := fillrange_color
        fillalpha := 0.5
        fillrange := quantile.(dists, 0.90)
        label := ""
        return forec_idx, quantile.(dists, 0.10)
    end
end

RecipesBase.@recipe function f(model::StateSpaceModel, scenarios::Array{<:AbstractFloat, 3})
    if !isunivariate(model)
        error("This plot recipe currently works for univariate models only.")
    end
    n = length(model.system.y)
    n_scen = size(scenarios, 3)
    scenarios_idx = collect(n + 1: n + size(scenarios, 1))
    # Plot the series
    @series begin
        seriestype := :path
        seriescolor := "black"
        return model.system.y
    end
    # Plot the scenarios
    labels = ["" for s in 1:n_scen - 1]
    pushfirst!(labels, "scenarios")
    @series begin
        seriestype := :path
        seriescolor := "grey"
        linewidth := 0.2
        label := permutedims(labels)
        return scenarios_idx, scenarios[:, 1, :]
    end
end


text:
function RichardsonExpansion(arg0::CR3BPSystem, arg1::LagrangianPoints)
    return RichardsonExpansion((CR3BPSystem, LagrangianPoints), arg0, arg1)
end

function compute_halo_first_guess(obj::RichardsonExpansion, arg0::jdouble, arg1::LibrationOrbitFamily, arg2::jdouble, arg3::jdouble)
    return jcall(obj, "computeHaloFirstGuess", PVCoordinates, (jdouble, LibrationOrbitFamily, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function compute_lyapunov_first_guess(obj::RichardsonExpansion, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(obj, "computeLyapunovFirstGuess", PVCoordinates, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function get_cr3bp_system(obj::RichardsonExpansion)
    return jcall(obj, "getCr3bpSystem", CR3BPSystem, ())
end

function get_halo_orbital_period(obj::RichardsonExpansion, arg0::jdouble)
    return jcall(obj, "getHaloOrbitalPeriod", jdouble, (jdouble,), arg0)
end

function get_lagrangian_point(obj::RichardsonExpansion)
    return jcall(obj, "getLagrangianPoint", LagrangianPoints, ())
end

function get_lyapunov_orbital_period(obj::RichardsonExpansion, arg0::jdouble)
    return jcall(obj, "getLyapunovOrbitalPeriod", jdouble, (jdouble,), arg0)
end


text:
# Is z within an isapprox of any x[i] ?
# Assumes that x is sorted from lowest to highest
function nearby(z::T, x::Vector{T}) where {T}
    top = length(x)
    bot = 1
    mid = div(top, 2)
    while bot < mid < top
        if z > x[mid]
            bot = mid
        elseif z < x[mid]
            top = mid
        else
            return mid
        end
        mid = bot + div(top-bot, 2)
    end
    isapprox(x[bot], z) && return bot
    isapprox(x[top], z) && return top
    return -1   # should error
end

# For general location points with precalculated weights
"""
    bary(z, f, x, w)

evaluate f(z)

# Arguments
- `z::Float64`:         the point at which to evaluate f
- `f::Vector{Float64}`: vector of function values at x
- `x:Vector{Float64}`:  vector of eval locations of f (sorted)
- `w:Vector{Float64}`:  weights for locations x
"""
function bary(z::T, f::Vector{T}, x::Vector{T}, w::Vector{T}) where {T}
    # assert(length(f) == length(x) == length(w))
    num = zero(T)
    den = zero(T)
    @inbounds for j = 1:length(f)
        t = w[j] / (z - x[j])
        num += t * f[j]
        den += t
    end
    fz = num / den
    fz = isnan(fz) ? f[nearby(z, x)] : fz
end

function bary(z, f)
    n = length(f)
    x = chebpts(n)
    num = den = zero(z)
    t = 1 / 2(z - x[1])
    num += t * f[1]
    den += t
    sgn = -1
    @inbounds for j = 2:n-1
        t = sgn / (z - x[j])
        num += t*f[j]
        den += t
        sgn = -sgn
    end
    t = sgn / 2(z - x[n])
    num += t * f[n]
    den += t
    fz = num / den
    fz = isnan(fz) ? f[nearby(z, x)] : fz
end


text:
using OceanParameterizations, Plots

𝒟 = data("strong_wind", reconstruct_fluxes=false)
𝒟_reconstructed = data("strong_wind", reconstruct_fluxes=true)
z = 𝒟_reconstructed.uw.z
t = 𝒟_reconstructed.t
Nt = length(𝒟.t)
output_gif_directory = "TestReconstructFluxes"
animate_gif((𝒟_reconstructed.uw.coarse, 𝒟.uw.coarse[:,1:Nt-1]), z, t, "uw",
            x_label=["∫(-du/dt + fv)dz", "truth"],
            filename="uw_reconstructed",
            directory=output_gif_directory)
animate_gif((𝒟_reconstructed.vw.coarse, 𝒟.vw.coarse[:,1:Nt-1]), z, t, "vw",
            x_label=["∫(-dv/dt - fu)dz", "truth"],
            filename="vw_reconstructed",
            directory=output_gif_directory)
animate_gif((𝒟_reconstructed.wT.coarse, 𝒟.wT.coarse[:,1:Nt-1]), z, t, "wT",
            x_label=["∫(-dw/dt)dz", "truth"],
            filename="wT_reconstructed",
            directory=output_gif_directory)


text:
using Plots
x = 1:10; y = rand(10); # These are the plotting data
plot(x, y);

println("HELLO")

function f(z)
    a = sum(z.^2,dims=2)
    return a
end

println(f([1 2; 3 4]))

l = [1 2; 3 4]
println(l[1,1])
println(l[2,1])
println(l[1,2])

for a in 1:10
end

bbox = [0 0; 2 2]
h0 = 1
x = bbox[1,1]:h0:bbox[2,1]
y = bbox[1,2]:h0:bbox[2,2]
println(x)
println(y)

# Funkcja meshgrid znana z Matlaba jest niedostępna w Julii
# Tworzoymy jej podstawową implementację
# Obliczamy długość podanych w parametrach funkcji wektorów
# Reprezentują one wymiarny naszej siatki
# Nastepnie funkcja repeat tworzy wektor skladajacy sie z
# Elemenntów pierwotnego wektora, ale kazdy element zostal powtórzon????????????y m ilość razy
# Co finalnie daje nam długi wektor wspolrzednych x punktów naszej siatki
# Funkcja reshape łamie wektor do macierzy o wymiarach m na na
# Ta sama prodecura powtórzona jest dla punktów y tylko w tym przypadku
# Funkcja repeat powtarza cały wektor y n razy, na koniec łamie go do macierzy
# Zwracana jest para macierzy gx, gy

function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T
    m, n = length(vy), length(vx)

    gx = reshape(repeat(vx, inner = m, outer = 1), m, n)
    gy = reshape(repeat(vy, inner = 1, outer = n), m, n)

    return gx, gy
end

println(repeat([1 2; 3 4], inner=(2,1), outer=(1, 3)))
println(repeat([1 2; 3 4], inner=(2,3)))
#println(meshgrid(x))

p =

max_stars_repo_name: jstarczewski/julia-fem

max_issues_repo_name: jstarczewski/julia-fem

max_forks_repo_name: jstarczewski/julia-fem


text:
using Test
using YaoExperiment: yaorepl_parse, yaorepl_trans,  YaoREPLState

@testset "yaorepl_parse" begin
    @test yaorepl_parse("e op reg") == :(expect(op, reg))
    @test yaorepl_parse("e' op reg") == :(expect'(op, reg))
    @test yaorepl_parse("m! reg") == :(measure!(reg))
    @test yaorepl_parse("cr(x)") == :(cr(x))
    @test yaorepl_parse("3=>X") == :(3=>X)
    @test yaorepl_parse("reg |> X") == :(reg |>X)
end

@testset "yaorepl_trans" begin
    s = YaoREPLState()
    @test yaorepl_trans(:(expect(op, reg)), s) == :(expect(op, reg))
    @test yaorepl_trans(:(expect'(op, reg)), s) == :(expect'(op, reg))
    @test yaorepl_trans(:(measure!(reg)), s) == :(measure!(reg))
    @test yaorepl_trans(:(help()), s) == :(help())
    @test_throws ErrorException yaorepl_trans(:(code()), s) == :(code())
    @test_throws ErrorException yaorepl_trans(:(cr(x)), s)
    @test yaorepl_trans(:(reg = zero_state(5)), s) == :(reg = zero_state(5))
    @test yaorepl_trans(:(reg = zero_state(5, nbatch=10)), s) == :(reg = zero_state(5, nbatch=10))
    @test yaorepl_trans(:(reg |>X), s) == :(apply!(reg, X))
    @test yaorepl_trans(:(op = begin nqubits=3; 2=>X; end), s) == :(op = chain(3, [put(3,2=>X)]))
end


text:
export InformedEditAllocationMove

struct InformedEditAllocationMove <: InvMcmcMove
    ν::Int
    P::Matrix{Float64}   #
    counts::Vector{Int}  # Track acceptance
    function InformedEditAllocationMove(;ν::Int=3)
        new(
            ν,
            # zeros(Int,ν),
            # [0],
            [0,0]
        )
    end
end

Base.show(io::IO, x::EditAllocationMove) = print(io, "EditAllocationMove(ν=$(x.ν))")


text:
using MTIWrapper
using Test
using CSV

cwd = dirname(@__FILE__)

mti_dir = "$cwd/tmp"
MTIWrapper.install_web_api(mti_dir)
test_df = CSV.read("$cwd/test_df.txt")

@testset "Title and Abstract" begin
    ti_abs_file = "$cwd/ti-abs.txt"
    output_file = "$cwd/ti-abs-results.txt"

    rm(ti_abs_file, force=true)
    rm(output_file, force=true)

    MTIWrapper.title_abstracts_to_request_file(test_df, ti_abs_file)
    MTIWrapper.mti_batch_query(mti_dir, ENV["NCBI_EMAIL"], ENV["UMLS_USER"], ENV["UMLS_PSSWD"],
            ti_abs_file, output_file)
    header = ["uid","term", "cui", "score", "type", "misc", "location", "path"]
    df = CSV.read(output_file, header = header, datarow=3, delim='|')

    println(df)
    @test size(df) == (34,8)
end

@testset "Abstract" begin
    abs_file = "$cwd/abs.txt"
    output_file = "$cwd/abs-results.txt"

    rm(abs_file, force=true)
    rm(output_file, force=true)

    MTIWrapper.abstracts_to_request_file(test_df, abs_file)
    MTIWrapper.mti_batch_query(mti_dir, ENV["NCBI_EMAIL"], ENV["UMLS_USER"], ENV["UMLS_PSSWD"],
            abs_file, output_file)
    header = ["uid","term", "cui", "score", "type", "misc", "location", "path"]
    df = CSV.read(output_file, header = header, datarow=3, delim='|')

    println(df)
    @test size(df) == (28,8)
end


text:
import StringUtils

function is_leap_year(year::Int64)
    if rem(year, 400) == 0
        return true
    elseif rem(year, 100) == 0
        return false
    elseif rem(year, 4) == 0
        return true
    else
        return false
    end
end

function days_in_month(month::Int64, year::Int64)
    if month == 4 || month == 6 || month == 9 || month == 11
        return 30
    elseif month == 2
        if is_leap_year(year)
            29
        else
            28
        end
    else
        31
    end
end

function main()
    day = 1
    num = 0
    for year in 1900:2000
        for month in 1:12
            day += days_in_month(month, year)
            day = rem(day, 7)
            if day == 0 # Sunday on first of the month
                if year == 1900
                    num += (month == 12) ? 1 : 0
                elseif year == 2000
                    num += (month == 12) ? 0 : 1
                else
                    num += 1
                end
            end
        end
    end
    @printf("%d Sundays fell on the first of the month during the twentieth century (1 Jan 1901 to 31 Dec 2000)\n", num)
end

main()


text:
export BrickletOneWireSearchBusLowLevel
struct BrickletOneWireSearchBusLowLevel
    identifier_length::Integer
    identifier_chunk_offset::Integer
    identifier_chunk_data::Vector{Integer}
    status::Integer
end

export BrickletOneWireRead
struct BrickletOneWireRead
    data::Integer
    status::Integer
end

export BrickletOneWireSPITFPErrorCount
struct BrickletOneWireSPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletOneWireIdentity
struct BrickletOneWireIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletOneWireSearchBus
struct BrickletOneWireSearchBus
    identifier::Vector{Integer}
    status::Integer
end

export BrickletOneWire
"""
Communicates with up 64 1-Wire devices
"""
mutable struct BrickletOneWire <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletOneWire(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_one_wire")
        deviceInternal = package.BrickletOneWire(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export search_bus_low_level
"""
    $(SIGNATURES)

Returns a list of up to 64 identifiers of the connected 1-Wire devices.
Each identifier is 64-bit and consists of 8-bit family code, 48-bit ID and
8-bit CRC.

To get these identifiers the Bricklet runs the
`SEARCH ROM algorithm <https://www.maximintegrated.com/en/app-notes/index.mvp/id/187>`__,
as defined by Maxim.
"""
function search_bus_low_level(device::BrickletOneWire)
    return device.deviceInternal.search_bus_low_level()
end

export reset_bus
"""
    $(SIGNATURES)

Resets the bus with the 1-Wire reset operation.
"""
function reset_bus(device::BrickletOneWire)
    return device.deviceInternal.reset_bus()
end

export write
"""
    $(SIGNATURES)

Writes a byte of data to the 1-Wire bus.
"""
function write(device::BrickletOneWire, data)
    return device.deviceInternal.write(data)
end

export read
"""
    $(SIGNATURES)

Reads a byte of data from the 1-Wire bus.
"""
function read(device::BrickletOneWire)
    return device.deviceInternal.read()
end

export write_command
"""
    $(SIGNATURES)

Writes a command to the 1-Wire device with the given identifier. You can obtain
the identifier by calling :func:`Search Bus`. The MATCH ROM operation is used to
write the command.

If you only have one device connected or want to broadcast to all devices
you can set the identifier to 0. In this case the SKIP ROM operation is used to
write the command.
"""
function write_command(device::BrickletOneWire, identifier, command)
    return device.deviceInternal.write_command(identifier, command)
end

export set_communication_led_config
"""
    $(SIGNATURES)

Sets the communication LED configuration. By default the LED shows 1-wire
communication traffic by flickering.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is off.
"""
function set_communication_led_config(device::BrickletOneWire, config)
    device.deviceInternal.set_communication_led_config(config)
end

export get_communication_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Communication LED Config`
"""
function get_communication_led_config(device::BrickletOneWire)
    return device.deviceInternal.get_communication_led_config()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletOneWire)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletOneWire, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletOneWire)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletOneWire, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletOneWire, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletOneWire, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletOneWire)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletOneWire)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletOneWire)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletOneWire, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletOneWire)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletOneWire)
    return device.deviceInternal.get_identity()
end

export search_bus
"""
Returns a list of up to 64 identifiers of the connected 1-Wire devices.
Each identifier is 64-bit and consists of 8-bit family code, 48-bit ID and
8-bit CRC.

To get these identifiers the Bricklet runs the
`SEARCH ROM algorithm <https://www.maximintegrated.com/en/app-notes/index.mvp/id/187>`__,
as defined by Maxim.
"""
function search_bus(device::BrickletOneWire)

    lock(device.stream_lock) do
        ret = search_bus_low_level(device, )
            identifier_length = ret.identifier_length
        identifier_out_of_sync = ret.identifier_chunk_offset != 0

        identifier_data = ret.identifier_chunk_data

        while !identifier_out_of_sync && length(identifier_data) < identifier_length
            ret = search_bus_low_level(device, )
                identifier_length = ret.identifier_length
            identifier_out_of_sync = ret.identifier_chunk_offset != length(identifier_data)
            identifier_data += ret.identifier_chunk_data
        end

        if identifier_out_of_sync # discard remaining stream to bring it back in-sync
            while ret.identifier_chunk_offset + 7 < identifier_length
                ret = search_bus_low_level(device, )
                    identifier_length = ret.identifier_length
            end

            throw(TinkerforgeStreamOutOfSyncError("Identifier stream is out-of-sync"))
        end
    end

        return SearchBus(identifier_data[:identifier_length], ret.status)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct AppsPermissionsApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_appsPermissionsInfo(_api::AppsPermissionsApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/apps.permissions.info", ["slackAuth"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Returns list of permissions this app has on a team.
Return: Dict{String, Any}
"""
function appsPermissionsInfo(_api::AppsPermissionsApi; _mediaType=nothing)
    _ctx = _swaggerinternal_appsPermissionsInfo(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function appsPermissionsInfo(_api::AppsPermissionsApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_appsPermissionsInfo(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_appsPermissionsRequest(_api::AppsPermissionsApi, scopes::String, trigger_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/apps.permissions.request", ["slackAuth"])
    Swagger.set_param(_ctx.query, "scopes", scopes)  # type String
    Swagger.set_param(_ctx.query, "trigger_id", trigger_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Allows an app to request additional scopes
Param: scopes::String (required)
Param: trigger_id::String (required)
Return: Dict{String, Any}
"""
function appsPermissionsRequest(_api::AppsPermissionsApi, scopes::String, trigger_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_appsPermissionsRequest(_api, scopes, trigger_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function appsPermissionsRequest(_api::AppsPermissionsApi, response_stream::Channel, scopes::String, trigger_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_appsPermissionsRequest(_api, scopes, trigger_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export appsPermissionsInfo, appsPermissionsRequest


text:
import Base.Test: @testset, @test, @test_throws
import StaticArrays: SVector

using MultipleScattering

concen_particles1 = [
    Particle(Acoustic(2.0,2.0,2),Sphere((0.0,0.0),2.0)),
    Particle(Acoustic(0.2,0.2,2),Sphere((0.0,0.0),1.0))
]
concen_particles2 = [
    Particle(Acoustic(2.0,1.0,2),Sphere((6.0,3.0),1.)),
    Particle(Acoustic(0.4,0.8,2),Sphere((6.0,3.0),1.4))
]
particle = Particle(Acoustic(1.5,1.5,2),Sphere((-6.0,-3.0),0.8))

ps = [CapsuleParticle(concen_particles2...), CapsuleParticle(concen_particles1...), particle]
# plot(ps)

medium = Acoustic(0.8, 0.5 + 0.0im,2)

ps = [CapsuleParticle(concen_particles1...)]
source = plane_source(medium, SVector(-3.0,0.0), SVector(1.0,0.0), 1.)
sim = FrequencySimulation(ps, source)

using Plots; #pyplot()

ω  = 1.4
bounds = Rectangle([-3.,-2.],[3.,2.]);
simres = run(sim,bounds,[ω]; res=40, basis_order=8)

plot(simres,ω,seriestype=:contour)
plot!(sim)

ωs = 0.:0.01:2.0
# simres = run(sim,bounds,ωs; res=30, basis_order =8);
simres = run(sim,bounds,ωs; res=15, basis_order = 26);

ts = 0.:0.4:50.
timres = TimeSimulationResult(simres; t_vec = ts)
#timres = run(sim,bounds,0.:0.02:2.0; ts = ts, res=20)

t=10.0
anim = @animate for t in ts
    # plot(timres,t,seriestype=:contour, clim=(-0.8,0.8), c=:balance)
    plot(timres,t,seriestype=:contour, clim=(0.0,0.8), c=:balance)
    plot!(sim)
end
#
gif(anim,"time_capsule_balance.gif", fps = 4)
# gif(anim,"time_capsule.gif", fps = 4)

concen = CapsuleParticle(
    Particle(Acoustic(2.0,2.0,2),Sphere((0.0,0.0),2.0)),
    Particle(Acoustic(0.2,0.2,2),Sphere((0.0,0.0),1.0))
)

n=12
ps = [Particle(Acoustic(0.1,0.1,2),Sphere((4*cos(θ),4*sin(θ)),0.5)) for θ in LinRange(0.,2pi,n+1)[1:n] ]
ps = [concen; ps]
plot(ps)

ω = 1.4
source_position = SVector(0.0,-3.0)
amplitude = 1.0
source = point_source(medium, source_position, amplitude)
sim = FrequencySimulation(ps, source)
sim_source = FrequencySimulation(source)

ω = 0.01
# checkout: ω = 0.63, strange resonance
bounds = Rectangle([-5.,-4.5],[5.,5.]);
simres = run(sim,bounds,[ω]; res=45, basis_order=3)

plot(simres,ω,seriestype=:contour)
plot!(sim)

ωs = 0.:0.0005:3.0

# simres = run(sim,bounds,ωs; res=57, min_basis_order = 4, basis_order = 14);
simres = run(sim,bounds,ωs; res=57, min_basis_order = 5, basis_order = 14);

ts = 0.:0.4:50.
timres = TimeSimulationResult(simres; t_vec = ts)
#timres = run(sim,bounds,0.:0.02:2.0; ts = ts, res=20)

timres2 = TimeSimulationResult(timres.field .+ 0.008, timres.x, timres.t)

t=0.0
anim = @animate for t in ts
    # plot(timres,t,seriestype=:contour, clim=(-0.8,0.8), c=:balance)
    plot(timres2,t,seriestype=:contour, clim=(-0.05,0.15), c=:balance)
    plot!(sim)
end

gif(anim,"capsulse_dance.gif", fps = 4)

ω = 0.65
anim = @animate for ω in 0.:0.01:2.
    # plot(timres,t,seriestype=:contour, clim=(-0.8,0.8), c=:balance)
    plot(simres,ω,seriestype=:contour, clim=(-1.1,1.1), c=:balance)
    plot!(sim)
end

gif(anim,"capsulse_dance_modes.gif", fps = 4)


text:
export ==, !=, hash, isless

#############################
## Degree
#############################
function isless(deg1::Degree, deg2::Degree)
  return (deg1.explvar<deg2.explvar) || (deg1.explvar==deg2.explvar && deg1.conjvar<deg2.conjvar)
end

#############################
## Variable
#############################
function isless(a::Variable, b::Variable)
  return isless(a.name, b.name)
end

#############################
## Exponent
#############################
function isless(exp1::Exponent, exp2::Exponent)
    ## Handle one exponent, empty structures are not handled by following general sort alg.
    if length(exp1) == 0        # if exp1 == 1
        return length(exp2) > 0
    end

    if length(exp2) == 0        # exp1!=1, if exp2 == 1
        return false
    end

    state1 = start(exp1)
    state2 = start(exp2)
    while !done(exp1, state1) && !done(exp2, state2)
        (i1, state1) = next(exp1, state1)
        (i2, state2) = next(exp2, state2)

        if isequal(i1, i2)
            continue
        else
            return isless(i1,i2)
        end
        if done(exp1, state1) && done(exp2, state2)
            return false
        elseif done(exp1, state1) && !done(exp2, state2)
            return true
        elseif !done(exp1, state1) && done(exp2, state2)
            return false
        end
    end

    return false
end

"""
    isless_degree(exp1::Exponent, exp2::Exponent)

    Order sorting elements on their sum of degrees at a first level.
    Test show performance is comparable to previous sorting function.
"""
function isless_degree(exp1::Exponent, exp2::Exponent)
    # First order level: sum of degrees
    exp1_deg = 0
    state1 = start(exp1)
    while !done(exp1, state1)
        (i1, state1) = next(exp1, state1)
        exp1_deg += i1[2].explvar + i1[2].conjvar
    end

    exp2_deg = 0
    state2 = start(exp2)
    while !done(exp2, state2)
        (i2, state2) = next(exp2, state2)
        exp2_deg += i2[2].explvar + i2[2].conjvar
    end

    if exp1_deg < exp2_deg
        return true
    elseif exp1_deg > exp2_deg
        return false
    end

    # Second order level: sort with variables and degrees
    state1 = start(exp1)
    state2 = start(exp2)
    while !done(exp1, state1) && !done(exp2, state2)
        (i1, state1) = next(exp1, state1)
        (i2, state2) = next(exp2, state2)

        if isequal(i1, i2)
            continue
        else
            return isless(i1,i2)
        end
        if done(exp1, state1) && done(exp2, state2)
            return false
        elseif done(exp1, state1) && !done(exp2, state2)
            return true
        elseif !done(exp1, state1) && done(exp2, state2)
            return false
        end
    end

    return false
end

#############################
## Polynomial
#############################
# NOTE: order on polynomials ?


text:
include("Network.jl")
#Sets
set_generators = 1:12
set_nodes = 1:24
set_nodes_ref=1
set_nodes_noref=2:24
set_nodes_g= [1,2,7,13,15,15,16,18,21,22,23,23]
set_nodes_d= [1,2,3,4,5,6,7,8,9,10,13,14,15,16,18,19,20]
set_demands = 1:17

#Parameters
generators=[152	  13.32
            152	  13.32
            350 	20.7
            591	  20.93
            60	  26.11
            155	  10.52
            155	  10.52
            400	  6.02
            400  	5.47
            300  	0
            310	  10.52
            350	  10.89]

demands=[67.48	   32
         60.38	   35
         111.88	   40
         46.17	   33
         44.40	   32
         85.24	   30
         78.14	   34
         106.55	   29
         108.33	   30
         120.76	   34
         165.15	   31
         120.76	   35
         197.12	   37
         62.15	   32
         207.77	   35
         113.65	   32
         79.91	  37]

capacity=1
varcost=2
load_node=1
utility_node=2

network=hcat(From,To,LineCapacity,Susceptance)

start_node=network[:,1]
end_node=network[:,2]
F_max=network[:,3]
B=network[:,4]

n_node=24
n_link=68

links = Array{Tuple{Int64, Int64},2}(undef, n_link, 1)

for i=1:n_link
  links[i] = (start_node[i], end_node[i])
end

links_rev = Array{Tuple{Int64, Int64},2}(undef, n_link, 1)

for i=1:n_link
  links_rev[i] = (end_node[i], start_node[i])
end

F_max_dict=Dict()
B_dict=Dict()

for i=1:n_link
  F_max_dict[(start_node[i], end_node[i])] = F_max[i]
  B_dict[(start_node[i], end_node[i])] = B[i]
end

MapG=[(1,1),(2,2),(3,7),(4,13),(5,15),(6,15),(7,16),(8,18),(9,21),(10,22),(11,23),(12,23)]
MapD=[(1,1),(2,2),(3,3),(4,4),(5,5),(6,6),(7,7),(8,8),(9,9),(10,10),(11,13),(12,14),(13,15),(14,16),(15,18),(16,19),(17,20)]


text:
mutable struct State
  rng::AbstractRNG
  parameters::Parameters
  float::Vector{Float32}
  integer::Vector{Int32}
  boolean::Vector{Bool}
  name::Vector{Symbol}
  code::Vector{Any}
  exec::Vector{Any}
  instructions::Dict{Symbol, Function}
  definitions::Dict{Symbol, Any}
  flag_quote_name::Bool
  flag_quote_code::Bool

  State() = new(MersenneTwister(0), Parameters(), Float32[], Int32[], Bool[],
    Any[], Any[], Any[], Dict{Symbol, Function}(), Dict{Symbol, Any}(), false, false)
end

function pp_stacks(s::State)
  println("FLOAT: ($(join(s.float, " ")))")
  println("INTEGER: ($(join(s.integer, " ")))")
  println("BOOLEAN: ($(join(s.boolean, " ")))")
  println("NAME: ($(join(s.name, " ")))")
  println("CODE: ($(join(s.code, " ")))")
  println("EXEC: ($(join(s.exec, " ")))")
end

# Enables a registered Push instruction, specified by its name,
# by inserting it into the instruction table of a given state.
enable_instruction!(s::State, i::String) = enable_instruction!(s, Symbol(i))
enable_instruction!(s::State, i::Symbol) = s.instructions[i] = fetch_instruction(i)

# Returns a list with all enabled instructions and definitions within
# a given state.
list_instructions(s::State) =
  vcat(collect(keys(s.instructions)), collect(keys(s.definitions)))


text:
using Test

@time begin
@testset "PotentialUQ.jl" begin
    # include("test_mcmc.jl")
    include("test_hmc.jl")

    # include("snap_test_MAP.jl")
    # include("snap_test_sample.jl")
    # include("snap_GaN_sample.jl")
end
end


text:
using Test, MLToolkit
using Distributions: Normal
using MLToolkit.DistributionsX: test_stat
using Statistics: mean, var
using Flux: gpu

@testset "Noise" begin
    dim1, dim2 = 10, 5
    n_samples, atol = 5_000, 0.01

    @testset "UniformNoise" begin
        noise = UniformNoise(dim1) |> gpu
        @test size(noise) == (dim1,)
        @test size(rand(noise)) == size(noise)
        @test size(rand(noise, dim2)) == (size(noise)..., dim2)
        @test size(rand(UniformNoise(dim1, dim2))) == (dim1, dim2)

        x = test_stat(mean, noise, n_samples, atol * prod(size(noise)))
        test_stat(var, noise, x, atol * prod(size(noise)))

        @test logpdf(noise, x) ≈ log.(ones(dim1, n_samples) / 2)
    end

    @testset "GaussianNoise" begin
        noise = GaussianNoise(dim1) |> gpu
        @test size(noise) == (dim1,)
        @test size(rand(noise)) == size(noise)
        @test size(rand(noise, dim2)) == (size(noise)..., dim2)
        @test size(rand(GaussianNoise(dim1, dim2))) == (dim1, dim2)

        x = test_stat(mean, noise, n_samples, atol * prod(size(noise)))
        test_stat(var, noise, x, atol * prod(size(noise)))

        @test logpdf(noise, x) ≈ map(xi -> logpdf(Normal(0, 1), xi), x)
    end
end


text:
using JuliaEuler
using Dates

@testset "Examples 1-24" begin

    @testset "Ex1" begin
        @test ex1(10) == 23
        @test ex1(16) == 60
        @test ex1(1_000) == 233_168
    end

    @testset "Ex2" begin
        @test ex2(90) == 44
        @test ex2(4_000_000) == 4613732
    end

    @testset "Ex3" begin
        @test ex3(13_195) == 29
        @test ex3(600851475143) == 6857
    end

    @testset "Ex4" begin
        @test ex4(2) == 9009
        @test ex4(3) == 906609
    end

    @testset "Ex5" begin
        @test ex5(10) == 2520
        @test ex5(20) == 232792560
    end

    @testset "Ex6" begin
        @test ex6(10) == 2640
        @test ex6(100) == 25164150
    end

    @testset "Ex7" begin
        @test ex7(6) == 13
        @test ex7(10_001) == 104743
    end

    @testset "Ex8" begin
        @test ex8(4) == 5832
        @test ex8(13) == 23514624000
    end

    @testset "Ex9" begin
        @test ex9(12) == 3*4*5
        @test ex9(1000) == 31875000
    end

    @testset "Ex10" begin
        @test ex10(10) == 17
        @test ex10(2_000_000) == 142913828922
    end

    @testset "Ex11" begin
        @test ex11(2) == 99*97
        @test ex11(4) == 70600674
    end

    @testset "Ex12" begin
        @test ex12(5) == 28
        @test ex12(500) == 76576500
    end

    @testset "Ex13" begin
        @test ex13(1) == "3710728753"
        @test ex13(100) == "5537376230"
    end

    @testset "Ex14" begin
        @test ex14(8) == 7
        @test ex14(999_999) == 837799
    end

    @testset "Ex15" begin
        @test ex15(3) == 20
        @test ex15(20) == 137846528820
    end

    @testset "Ex16" begin
        @test ex16(15) == 26
        @test ex16(1000) == 1366
    end

    @testset "Ex17" begin
        @test ex17() == 21124
    end

    @testset "Ex18" begin
        @test ex18(3) == 75+64+82
        @test ex18(15) == 1074
    end

    @testset "Ex19" begin
        @test ex19(Date(2000,1,1), Date(2000,12,1)) == 1
        @test ex19(Date(1901,1,1), Date(2000,12,1)) == 171
    end

    @testset "Ex20" begin
        @test ex20(10) == 27
        @test ex20(100) == 648
    end

    @testset "Ex21" begin
        @test ex21(284) == 504
        @test ex21(10000) == 31626
    end

    @testset "Ex22" begin
        @test ex22() == 871198282
    end

    @testset "Ex23" begin
        @test ex23(24) == sum(1:23)
        @test ex23(28123) == 4179871
    end

     @testset "Ex24" begin
        @test [ex24(split("012",""),i) for i in 1:6] == [12,21,102,120,201,210]
        @test ex24(split("0123456789",""), 1_000_000) == 2783915460
    end
end

max_stars_repo_name: binnisb/JuliaEuler

max_issues_repo_name: binnisb/JuliaEuler

max_forks_repo_name: binnisb/JuliaEuler


text:
import Optim.optimize
import Base: axpy!, scale!

export frank_wolfe_sketched

# type FrankWolfeParams<:OptParams
# 	maxiters::Int
# 	abstol::AbstractFloat
# 	stepsize::StepSizeRule
# end
# stop(p::FrankWolfeParams, eps) = eps < p.abstol ? true : false
# FrankWolfeParams() = FrankWolfeParams(50, 1e-10, DecreasingStepSize(2,1))

function frank_wolfe_sketched(z::AbstractArray, # starting point
	objective::Function,
	grad_objective::Function,
	constraint::Function, # evaluates constraint function
	alpha::Number, # bound on constraint function
	min_lin_st_constraint::Function, # solves min_{c(x)<=alpha} g \dot x (as a function of alpha and g)
	sketch::AbstractSketch = IdentitySketch(x),
	params::FrankWolfeParams = FrankWolfeParams(),
	ch::ConvergenceHistory = ConvergenceHistory("fw_sketched");
	LB = -Inf,
	UB = Inf,
	verbose = false)

	# initialize
	if verbose @printf("%10s%12s%12s%12s%12s%12s\n", "iteration", "UB", "LB", "gap", "rel gap", "time") end

	t = time()
	t0 = copy(t)
	for k=1:params.maxiters

		# function and gradient evaluation
		objval = objective(z)
		G = grad_objective(z)
		# G = A'*g, where g is the gradient wrt the compressed variable
		# so G.factors[1] is the compression operator A
		A = G.factors[1]
		g = vec(G.factors[2])

		# solve fenchel problem
		# Delta is solution to min_{c(X)<=alpha} grad(X_old) \dot X
		Delta, linearized_obj = min_lin_st_constraint(G, alpha)
		dz = A*Delta

		# check stopping condition
		UB = min(UB, objval)
 	   # <G, Delta - X> = <A'g, Delta - X> = <g, A*Delta - A*X> = <g, dz - z>
		LB = max(LB, objval + linearized_obj - dot(g, z))
		old_t, t = t, time()
		if verbose && k%10==0
			cx = constraint(z)
			@printf("%10d%12.4e%12.4e%12.4e%12.4e%12.4e\n", k, UB, LB, UB - LB, (UB - LB)/min(abs(LB), abs(UB)), t - t0)
		end
		update_ch!(ch, t - old_t; obj = UB, dual_obj = LB)
		stop(params, UB, LB) && break
		t = time()

		# choose step size with stepsize rule
		# eg a = 1/k
		a = step(params.stepsize, objective, z, z-dz;
			     objval=objval,
			     normgradsq=UB-LB)

		# take step
		# z = (1-a)*z + a*dz
		scale!(z, 1-a); axpy!(a, dz, z)

		cgd_update!(sketch, Delta, a)
	end
	return reconstruct(sketch)
end

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
module BenchSet

using BenchmarkTools
using Transducers

const SUITE = BenchmarkGroup()

itr = (x => y for x in 1:1000 for y in 1:x);
SUITE["xf"] = @benchmarkable Dict(eduction($itr))
SUITE["base"] = @benchmarkable Dict($itr)

end  # module
BenchSet.SUITE


text:
using FinElt
using FinElt.PlanarPoisson

mesh = read_msh_file("keyhole.msh")
#h, sr = shape_params(mesh)
essential_bc = [ "Gamma" ]
f(x) = 4.0
vp = VariationalProblem(mesh, essential_bc)
add_bilin_form!(vp,  "Omega", grad_dot_grad!)
add_lin_functnl!(vp, "Omega", source_times_func!, f)

A, b = assembled_linear_system(vp)

ufree = A \ b

u = complete_soln(ufree, vp)
write_pos_file("keyhole.pos", "keyhole.msh") do fid
    save_warp_nodal_scalar_field(u, "u", fid)
end


text:
#GOALS
#=
1 Get the range of the dataset: dataCountryYearsNum() -> #looks at the whole dataset of score files and reads the header to return a dictionary for the year->countrynumber, min year, max year
2 Check the parameters: paramCheck() -> all the variable values supplied are valid
3 Simulate the null unbiases score distribution in a simulation: scoreSimDist -> Given the full range of years and interval/windowSize return the low2high score accumulation for each interval in a di=ctionary: windowDist[string(yr,"-",yr+windowSize)]
4 using the scoreSim dictionary of ranked list of sampled scores, return the threshold for each window as a dictionary -> windowConf[string(yr,"-",yr+windowSize)] = confalpha
5
=#

#need GraphViz to visualize the results
#using GraphViz
using Statistics

#csvDir2AdjList.jl is the file with functions to call upon the real data to get the aggregates
include("csvDir2AdjList.jl")
#include("biasGraphView.jl")
#include("biasAnalysis.jl")

#>>MAIN<<
function biasesESC(startYr = 1980, endYr = 1990, windowSize = 5, alpha = 0.05)

    #load data and get the dictionary for the country num per year
    #looks at the whole dataset of score files and reads the header to return a dictionary for the year->countrynumber,
    #min year, max year
    countryYearsNum, yrMin, yrMax = dataCountryYearsNum()

    #check params
    paramCheck(startYr, endYr, windowSize, yrMin, yrMax)

    #simulate scores to create a distribution1
    #return the low2high score accumulation for each interval in a dictionary: windowDist[string(yr,"-",yr+windowSize)]
    windowDist = scoreSimDist(startYr, endYr, windowSize, countryYearsNum)
    #println(windowDist)

    #first the upper
    tailSide = "upper"
    #get confidence intervals for the lower or upper end
    #using the scoresim, get a dictionary for the threshold of significances:
    #dictionary -> windowConf[string(yr,"-",yr+windowSize)] = confalpha
    windowConf = windowConfValues(startYr, endYr, windowSize, windowDist, tailSide, alpha)
    println("finished windowConf")
    #Now we must obtain the averages for each country from to country (we need to have the CSV data read)
    #return the dictionary of the time windows and the aggregate adjacency list of scores, country names, average scores aggregates for each window, and the thresholds of when the average surpasses the alpha value for significance. keys: "countries" "thresholdSigAdjList" "avgScoreAggregateAdjList" "scoreAggregateAdjList"
    winAggDictUpper = windowsDictThresholdsAdjList(windowConf, startYr, endYr, windowSize)
    println(winAggDictUpper)
    println("finished Agg Dict Upper")
    #To reduce computations later on, we now add to the dictionary a key to the total threshold surpassings as a count for the set of windows (1 per window significance) and get the total country name list for the full time period each window covers
    winAggDictUpper = dictTotalThresholdsAdjListWindowCount(winAggDictUpper,startYr,endYr,windowSize)

    println("finished upper")

    #same as above but now for the lower end of the distribution
    tailSide = "lower"
    windowConfLower = windowConfValues(startYr, endYr, windowSize, windowDist, tailSide, alpha)
    winAggDictLower = windowsDictThresholdsAdjList(windowConfLower, startYr, endYr, windowSize)
    winAggDictLower = dictTotalThresholdsAdjListWindowCount(winAggDictLower,startYr,endYr,windowSize)

    winAggDictUpper["alpha"] = alpha
    winAggDictLower["alpha"] = alpha
    winAggDictUpper["windowSize"] = windowSize
    winAggDictLower["windowSize"] = windowSize
    winAggDictUpper["side"] = "Upper"
    winAggDictLower["side"] = "Lower"

    println("finished lower")

    #produce the GraphViz visualizations for this dataset; from ("biasGraphView.jl")
#    graphAvoid(winAggDictUpper,winAggDictLower)
#    println("finished graphing")
#    #produce scatter plots for the aggregate edge types for countries; from ("biasAnalysis.jl")
#    analyzeBiases(winAggDictUpper,winAggDictLower)
#    println("finished scatter plots")
    return [winAggDictUpper,winAggDictLower]
end

#return the year windows' confidence interval values specified
#using the scoreSim dictionary of ranked list of sampled scores, return the threshold for each window as a dictionary -> windowConf[string(yr,"-",yr+windowSize)] = confalpha
function windowConfValues(startYr, endYr, windowSize, windowDist, tailSide, alpha)

    if(tailSide == "upper" || tailSide == "right")
        alpha = 1 - alpha
    end

    windowConf = Dict() #hold the distribution of the scores
    if(tailSide == "upper"  || tailSide == "right")
        windowConf["tailSide"] = "upper"
    else
        windowConf["tailSide"] = "lower"
    end

    yr = startYr
    while( (yr+windowSize) <= endYr )
        #each window is a null dist unique due to the voting schemes
        #simulate the scores for the dist
        distTmp = windowDist[string(yr,"-",yr+windowSize)]
        sampleSize = length(distTmp)
        confIndAlpha =  max(1,floor(Int,alpha*sampleSize))
        confalpha = distTmp[confIndAlpha]

        windowConf[string(yr,"-",yr+windowSize)] = confalpha
        yr = yr + windowSize
    end

    return windowConf

end

#return the year windows of distributions for scores
#Given the full range of years and interval/windowSize return the low2high score accumulation for each interval in a dictionary: windowDist[string(yr,"-",yr+windowSize)]
function scoreSimDist(startYr, endYr, windowSize, countryYearsNum)

    #Generate NULL distribution for each set of years in the windows
    windowDist = Dict() #hold the distribution of the scores
    yr = startYr
    while( (yr+windowSize) <= endYr )
        #each window is a null dist unique due to the voting schemes
        #simulate the scores for the dist
        distTmp = scoreSim(yr,yr+windowSize,countryYearsNum)
        windowDist[string(yr,"-",yr+windowSize)] = distTmp
        yr = yr + windowSize
    end

    return windowDist
end

#simulate the score distribution for each window span provided
function scoreSim(startYr,endYr,countryYearsNum)
    AVG_SIMULATION = []
    iterNum = 3000
    for ii = 1:iterNum
        ONE_SIMULATION = []
        for yr = startYr:endYr
            NUM = countryYearsNum[yr]
	    if(yr >= 1975 || yr == 1963 || yr == 1962)
                score = Allocated(yr,NUM);

            elseif( (1964<=yr<= 1966) || yr==1974 || (1967<=yr<=1970) || (1957<=yr<=1961))
                score = Sequential(yr,NUM)
            elseif(1971<=yr<=1973)
                score = Rated(yr,NUM)
            else
                score = Allocated(-1,NUM)
            end
            append!(ONE_SIMULATION,score);

        end
        avgSim = mean(ONE_SIMULATION);
        append!(AVG_SIMULATION,avgSim);
    end
    sortedAVG_SIMULATION = sort(AVG_SIMULATION,rev=false)#returns low to high

    return sortedAVG_SIMULATION
end

#here each country can receive a set of scores with consecutive points awarded
#in sequence for that year it has an equal chance of receiving each score
function Sequential(yr,NUM)
    SCORES1 = [5,3,1]
    SCORES2 = ones(Int,1,10)
    score = 0
    if(1964 <= yr <= 1966)
        for ii=1:length(SCORES1)
            position = ceil.(rand(1,1)*NUM)
            if((Int.(position))[1] == 1)
                score = SCORES1[ii] + score
            end
        end
    elseif(yr==1974 || (1967<=yr<=1970) || (1957<=yr<=1961))
        for ii=1:length(SCORES2)
            position = ceil.(rand(1,1)*NUM)
            if((Int.(position))[1] == 1)
                score = SCORES2[ii] + score
            end
        end
    end
    return score

end

function Allocated(yr,NUM)
    SCORES1 = [3,2,1]
    SCORES2 = [5,4,3,2,1]
    SCORES3 = [12,10,8,7,6,5,4,3,2,1]
    position = ceil.(rand(1,1)*NUM)
    if(yr >= 1975 && yr <= 2016)
        SCORES = SCORES3
    elseif(yr == 1962)
        SCORES = SCORES1
    elseif(yr == 1963)
        SCORES = SCORES2
    else
        SCORES = SCORES3
    end
    if(position[1] <= length(SCORES))
        score = SCORES[Int.(position)]
    else
        score = 0
    end
    return score
end

function Rated(yr,NUM)
    SCORES1 = [5,4,3,2,1]
    if(1971<=yr<=1973)
        X1 = SCORES1[rand(1:end)]
        X2 = SCORES1[rand(1:end)]
    else
        return -1
    end
    score = X1 + X2
    return score
end

#get the years of the data provided and country number in each year
#looks at the whole dataset of score files and reads the header to return a dictionary for the year->countrynumber, min year, max year
function dataCountryYearsNum()
    countryYearsNum = Dict{Integer,Integer}()
    resultsFile = readdir("./dataTables/")
    yrMin = 100000
    yrMax = -1
    for rf in resultsFile
        fileTmp = open(string("./dataTables/",rf))
        linesTmp = readlines(fileTmp)  #readfile lines
        yrTmp = parse(Int,((split(rf,"."))[1]))
        countryNumTmp = length(split(linesTmp[1],",")) - 1
        countryYearsNum[yrTmp] = countryNumTmp
        close(fileTmp)
        if(yrTmp < yrMin)
            yrMin = yrTmp
        end
        if(yrTmp > yrMax)
            yrMax = yrTmp
        end
    end

    return countryYearsNum, yrMin, yrMax
end

#fn to accept the parameters and check for the validity
function paramCheck(startYr, endYr, windowSize, yrMin, yrMax)
    if(startYr >= endYr)
        println("the start year needs to be before the end year")
        quit()
    end
    if((startYr+4) >= endYr)
        println("the start year needs to be at least 4 years prior the end year")
        quit()
    end
    #sanity check input years
    if( endYr < startYr || startYr < yrMin || endYr > yrMax )
        print(string("year range improperly set, for the analysis end year must be greater than start and the smallest year is $(yrMin) and largest $(yrMax) with smallest first"))
        quit()
    end
    if( (startYr+windowSize) > endYr)
        print("not enough years between start and end for analysis due to window size")
        quit()
    end

end

max_stars_repo_name: mantzaris/eurovisionBias

max_issues_repo_name: mantzaris/eurovisionBias

max_forks_repo_name: mantzaris/eurovisionBias


text:
# using Printf

"Parameters for right-hand side function of variational partitioned Runge-Kutta methods."
mutable struct ParametersVPRKdegenerate{DT, TT, ET <: IODE{DT,TT}, D, S} <: AbstractParametersVPRK{DT,TT,ET,D,S}
    equ::ET
    tab::TableauVPRK{TT}
    Δt::TT
    cache::NonlinearFunctionCacheVPRK{DT}

    t::TT
    q::Vector{DT}
    p::Vector{DT}
end

function ParametersVPRKdegenerate(equ::ET, tab::TableauVPRK{TT}, Δt::TT,
                                  cache::NonlinearFunctionCacheVPRK{DT}) where {DT, TT, ET <: IODE{DT,TT}}
    q = zeros(DT, equ.d)
    p = zeros(DT, equ.d)

    ParametersVPRKdegenerate{DT, TT, ET, equ.d, tab.s}(equ, tab, Δt, cache, 0, q, p)
end

"Compute solution of degenerate symplectic partitioned Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST},
                params::ParametersVPRKdegenerate{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    # create temporary vectors
    q̅ = zeros(ST,D)
    p̅ = zeros(ST,D)

    quote
        @assert length(x) == length(b)

        compute_solution!(x, $q̅, $p̅, params)

        # compute b = - [q̅-q-BV]
        for k in 1:div(D,2)
            b[k] = params.q[k] - $q̅[k]
            for i in 1:S
                b[k] += params.Δt * params.tab.q.b[i] * params.cache.V[k,i]
            end
        end

        # compute b = - [p̅-p-BF]
        for k in 1:div(D,2)
            b[div(D,2)+k] = params.p[k] - $p̅[k]
            for i in 1:S
                b[div(D,2)+k] += params.Δt * params.tab.p.b[i] * params.cache.F[k,i]
            end
        end
    end
end

function compute_solution!(
                x::Vector{ST}, q̅::Vector{ST}, p̅::Vector{ST},
                params::ParametersVPRKdegenerate{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    @assert length(x) == length(q̅) == length(p̅)

    # copy x to q
    q̅ .= x

    # compute p̅ = ϑ(q)
    params.equ.ϑ(params.t + params.Δt, q̅, p̅)

end

"Variational partitioned Runge-Kutta integrator."
mutable struct IntegratorVPRKdegenerate{DT, TT,
                SPT <: ParametersVPRK{DT,TT},
                PPT <: ParametersVPRKdegenerate{DT,TT},
                SST <: NonlinearSolver{DT},
                STP <: NonlinearSolver{DT},
                IT  <: InitialGuessPODE{DT,TT}} <: AbstractIntegratorVPRK{DT,TT}

    sparams::SPT
    pparams::PPT

    solver::SST
    projector::STP
    iguess::IT

    cache::NonlinearFunctionCacheVPRK{DT}

    q::Vector{Vector{TwicePrecision{DT}}}
    p::Vector{Vector{TwicePrecision{DT}}}
end

function IntegratorVPRKdegenerate(equation::ET, tableau::TableauVPRK{TT}, Δt::TT) where {DT, TT, ET <: IODE{DT,TT}}
    D = equation.d
    M = equation.n
    S = tableau.s

    # create cache for internal stage vectors and update vectors
    cache = NonlinearFunctionCacheVPRK{DT}(D,S)

    # create solver params
    sparams = ParametersVPRK(equation, tableau, Δt)

    # create projector params
    pparams = ParametersVPRKdegenerate(equation, tableau, Δt, cache)

    # create nonlinear solver
    solver = create_nonlinear_solver(DT, D*S, sparams)

    # create projector
    projector = create_nonlinear_solver(DT, D, pparams)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create solution vectors
    q = create_solution_vector(DT, D, M)
    p = create_solution_vector(DT, D, M)

    # create integrator
    IntegratorVPRKdegenerate{DT, TT, typeof(sparams), typeof(pparams), typeof(solver), typeof(projector), typeof(iguess)}(
                sparams, pparams, solver, projector, iguess, cache, q, p)
end

function initial_guess!(int::IntegratorVPRKdegenerate, m::Int)
    for i in 1:int.sparams.tab.s
        evaluate!(int.iguess, m, int.cache.y, int.cache.z, int.cache.v, int.sparams.tab.q.c[i], int.sparams.tab.p.c[i])
        for k in 1:int.sparams.equ.d
            int.solver.x[int.sparams.equ.d*(i-1)+k] = int.cache.v[k]
        end
    end
end

"Integrate ODE with variational partitioned Runge-Kutta integrator."
function integrate_step!(int::IntegratorVPRKdegenerate{DT,TT}, sol::SolutionPDAE{DT,TT}, m::Int, n::Int) where {DT,TT}
        # check if m and n are compatible with solution dimensions
        check_solution_dimension_asserts(sol, m, n)

        # set time and solution for nonlinear solver
        int.sparams.t = sol.t[0] + (n-1)*int.sparams.Δt
        int.sparams.q .= int.q[m]
        int.sparams.p .= int.p[m]

        # set time and solution for solution solver
        int.pparams.t = sol.t[0] + (n-1)*int.pparams.Δt
        int.pparams.q .= int.q[m]
        int.pparams.p .= int.p[m]

        # compute initial guess
        initial_guess!(int, m)

        # call nonlinear solver
        solve!(int.solver)

        # print solver status
        print_solver_status(int.solver.status, int.solver.params, n)

        # check if solution contains NaNs or error bounds are violated
        check_solver_status(int.solver.status, int.solver.params, n)

        # compute vector fields at internal stages
        compute_stages!(int.solver.x, int.cache.Q, int.cache.V, int.cache.P, int.cache.F, int.sparams)

        # compute internal stage solution
        update_solution!(int.q[m], int.cache.V, int.sparams.tab.q.b, int.sparams.tab.q.b̂, int.sparams.Δt)
        update_solution!(int.p[m], int.cache.F, int.sparams.tab.p.b, int.sparams.tab.p.b̂, int.sparams.Δt)

        # set initial guess for solution
        int.projector.x .= int.q[m]

        # call solution solver
        solve!(int.projector)

        # print solver status
        print_solver_status(int.projector.status, int.projector.params, n)

        # check if solution contains NaNs or error bounds are violated
        check_solver_status(int.projector.status, int.projector.params, n)

        # copy solution to integrator
        compute_solution!(int.projector.x, int.cache.q̅, int.cache.p̅, int.pparams)
        int.q[m] .= int.cache.q̅
        int.p[m] .= int.cache.p̅

        # copy solution to initial guess for next time step
        update!(int.iguess, m, sol.t[0] + n*int.sparams.Δt, int.q[m], int.p[m])

        # take care of periodic solutions
        cut_periodic_solution!(int.q[m], int.sparams.equ.periodicity)

        # copy to solution
        copy_solution!(sol, int.q[m], int.p[m], n, m)
end


text:
LShapeEigenfunction(domain::LShape; stride::Int = 1) =
    LShapeEigenfunction(domain, stride, arb[])

function Base.show(io::IO, u::LShapeEigenfunction)
    println(io, "Eigenfunction for the L-shaped domain")
    if !haskey(io, :compact) || !io[:compact]
        println(io, "domain: $(u.domain)")
        print(io, "number of set coefficients: $(length(u.coefficients))")
    end
end

active_boundaries(::LShape, ::LShapeEigenfunction) = 1:4

function (u::LShapeEigenfunction)(
    (r, θ)::Tuple{T,T},
    λ::arb,
    ks::UnitRange{Int};
    boundary = nothing,
) where {T<:Union{arb,arb_series}}
    rsqrtλ = r * sqrt(λ)
    res = similar(ks, T)
    for i in eachindex(ks)
        k = 1 + (ks[i] - 1) * u.stride

        ν = u.domain.parent(2k // 3)

        res[i] = bessel_j(ν, rsqrtλ) * sin(ν * θ)
    end

    return res
end

function (u::LShapeEigenfunction)(
    (r, θ)::Tuple{T,T},
    λ::arb;
    boundary = nothing,
) where {T<:Union{arb,arb_series}}
    coeffs = coefficients(u)
    return sum(coeffs .* u((r, θ), λ, 1:length(coeffs), boundary = boundary))
end

function norm(domain::LShape, u::LShapeEigenfunction, λ::arb)
    θ_integral = 3 // 4 * domain.parent(π)

    CC = ComplexField(domain.parent.prec)
    # The integrals goes from zero to the lower bound for θ. However
    # the function has a branch cut at zero and Arb has problem
    # handling this. We therefore integrate a small distance away from
    # zero.
    a = CC(1e-1)
    b = CC(1)

    r_integral = domain.parent(0)
    for k = 1:min(4, length(u.coefficients))
        ν::arb = domain.parent(2k // 3)
        c2 = u.coefficients[k]^2

        f = r -> begin
            r = real(r)
            if contains_negative(r)
                return CC(NaN)
            end

            CC(c2 * r * bessel_j(ν, sqrt(λ) * r)^2)
        end
        y = Nemo.integrate(CC, f, a, b)
        r_integral += real(y)
    end

    return sqrt(θ_integral * r_integral)
end


text:
"""
    RegularCartesianGrid{FT, TX, TY, TZ, R} <: AbstractGrid{FT, TX, TY, TZ}

A Cartesian grid with with constant grid spacings `Δx`, `Δy`, and `Δz` between cell centers
and cell faces, elements of type `FT`, topology `{TX, TY, TZ}`, and coordinate ranges
of type `R`.
"""
struct RegularCartesianGrid{FT, TX, TY, TZ, R} <: AbstractGrid{FT, TX, TY, TZ}
    # Number of grid points in (x,y,z).
    Nx :: Int
    Ny :: Int
    Nz :: Int
    # Halo size in (x,y,z).
    Hx :: Int
    Hy :: Int
    Hz :: Int
    # Domain size [m].
    Lx :: FT
    Ly :: FT
    Lz :: FT
    # Grid spacing [m].
    Δx :: FT
    Δy :: FT
    Δz :: FT
    # Range of coordinates at the centers of the cells.
    xC :: R
    yC :: R
    zC :: R
    # Range of grid coordinates at the faces of the cells.
    xF :: R
    yF :: R
    zF :: R
end

"""
    RegularCartesianGrid([FT=Float64]; size,
                         extent = nothing, x = nothing, y = nothing, z = nothing,
                         topology = (Periodic, Periodic, Bounded), halo = (1, 1, 1))

Creates a `RegularCartesianGrid` with `size = (Nx, Ny, Nz)` grid points.

Keyword arguments
=================

- `size` (required): A tuple prescribing the number of grid points in non-`Flat` directions.
                     `size` is a 3-tuple for 3D models, a 2-tuple for 2D models, and either a
                     scalar or 1-tuple for 1D models.

- `topology`: A 3-tuple `(Tx, Ty, Tz)` specifying the topology of the domain.
              `Tx`, `Ty`, and `Tz` specify whether the `x`-, `y`-, and `z` directions are
              `Periodic`, `Bounded`, or `Flat`. The topology `Flat` indicates that a model does
              not vary in that directions so that derivatives and interpolation are zero.
              The default is `topology=(Periodic, Periodic, Bounded)`.

- `extent`: A tuple prescribing the physical extent of the grid in non-`Flat` directions.
            The origin for three-dimensional domains is the oceanic default `(0, 0, -Lz)`.

- `x`, `y`, and `z`: Each of `x, y, z` are 2-tuples that specify the end points of the domain
                     in their respect directions. Scalar values may be used in `Flat` directions.

*Note*: _Either_ `extent`, or all of `x`, `y`, and `z` must be specified.

- `halo`: A tuple of integers that specifies the size of the halo region of cells surrounding
          the physical interior for each non-`Flat` direction.

The physical extent of the domain can be specified via `x`, `y`, and `z` keyword arguments
indicating the left and right endpoints of each dimensions, e.g. `x=(-π, π)` or via
the `extent` argument, e.g. `extent=(Lx, Ly, Lz)` which specifies the extent of each dimension
in which case 0 ≤ x ≤ Lx, 0 ≤ y ≤ Ly, and -Lz ≤ z ≤ 0.

A grid topology may be specified via a tuple assigning one of `Periodic`, `Bounded, and `Flat`
to each dimension. By default, a horizontally periodic grid topology `(Periodic, Periodic, Bounded)`
is assumed.

Constants are stored using floating point values of type `FT`. By default this is `Float64`.
Make sure to specify the desired `FT` if not using `Float64`.

Grid properties
===============

- `(Nx, Ny, Nz)::Int`: Number of physical points in the (x, y, z)-direction

- `(Hx, Hy, Hz)::Int`: Number of halo points in the (x, y, z)-direction

- `(Lx, Ly, Lz)::FT`: Physical extent of the grid in the (x, y, z)-direction

- `(Δx, Δy, Δz)::FT`: Cell width in the (x, y, z)-direction

- `(xC, yC, zC)`: (x, y, z) coordinates of cell centers.

- `(xF, yF, zF)`: (x, y, z) coordinates of cell faces.

Examples
========

* A default grid with Float64 type:

```jldoctest
julia> using Oceananigans

julia> grid = RegularCartesianGrid(size=(32, 32, 32), extent=(1, 2, 3))
RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}
                   domain: x ∈ [0.0, 1.0], y ∈ [0.0, 2.0], z ∈ [-3.0, 0.0]
                 topology: (Periodic, Periodic, Bounded)
  resolution (Nx, Ny, Nz): (32, 32, 32)
   halo size (Hx, Hy, Hz): (1, 1, 1)
grid spacing (Δx, Δy, Δz): (0.03125, 0.0625, 0.09375)
```

* A default grid with Float32 type:

```jldoctest
julia> using Oceananigans

julia> grid = RegularCartesianGrid(Float32; size=(32, 32, 16), x=(0, 8), y=(-10, 10), z=(-π, π))
RegularCartesianGrid{Float32, Periodic, Periodic, Bounded}
                   domain: x ∈ [0.0, 8.0], y ∈ [-10.0, 10.0], z ∈ [-3.1415927, 3.1415927]
                 topology: (Periodic, Periodic, Bounded)
  resolution (Nx, Ny, Nz): (32, 32, 16)
   halo size (Hx, Hy, Hz): (1, 1, 1)
grid spacing (Δx, Δy, Δz): (0.25f0, 0.625f0, 0.3926991f0)
```

* A two-dimenisional, horizontally-periodic grid:

```jldoctest
julia> using Oceananigans

julia> grid = RegularCartesianGrid(size=(32, 32), extent=(2π, 4π), topology=(Periodic, Periodic, Flat))
RegularCartesianGrid{Float64, Periodic, Periodic, Flat}
                   domain: x ∈ [0.0, 6.283185307179586], y ∈ [0.0, 12.566370614359172], z ∈ [0.0, 0.0]
                 topology: (Periodic, Periodic, Flat)
  resolution (Nx, Ny, Nz): (32, 32, 1)
   halo size (Hx, Hy, Hz): (1, 1, 0)
grid spacing (Δx, Δy, Δz): (0.19634954084936207, 0.39269908169872414, 0.0)
```

* A one-dimensional "column" grid:

```jldoctest
julia> using Oceananigans

julia> grid = RegularCartesianGrid(size=256, z=(-128, 0), topology=(Flat, Flat, Bounded))
RegularCartesianGrid{Float64, Flat, Flat, Bounded}
                   domain: x ∈ [0.0, 0.0], y ∈ [0.0, 0.0], z ∈ [-128.0, 0.0]
                 topology: (Flat, Flat, Bounded)
  resolution (Nx, Ny, Nz): (1, 1, 256)
   halo size (Hx, Hy, Hz): (0, 0, 1)
grid spacing (Δx, Δy, Δz): (0.0, 0.0, 0.5)
```
"""
function RegularCartesianGrid(FT=Float64;
                                  size,
                                     x = nothing, y = nothing, z = nothing,
                                extent = nothing,
                              topology = (Periodic, Periodic, Bounded),
                                  halo = nothing
                              )

    TX, TY, TZ = validate_topology(topology)
    size = validate_size(TX, TY, TZ, size)
    halo = validate_halo(TX, TY, TZ, halo)
    Lx, Ly, Lz, x, y, z = validate_regular_grid_domain(TX, TY, TZ, FT, extent, x, y, z)

    # Unpacking
    Nx, Ny, Nz = N = size
    Hx, Hy, Hz = H = halo
                 L = (Lx, Ly, Lz)
    Δx, Δy, Δz = Δ = L ./ N
                X₁ = (x[1], y[1], z[1])

    # Face-node limits in x, y, z
    xF₋, yF₋, zF₋ = XF₋ = @. X₁ - H * Δ
    xF₊, yF₊, zF₊ = XF₊ = @. XF₋ + total_extent(topology, halo, Δ, L)

    # Cell-node limits in x, y, z
    xC₋, yC₋, zC₋ = XC₋ = @. XF₋ + Δ / 2
    xC₊, yC₊, zC₊ = XC₊ = @. XC₋ + L + Δ * (2H - 1)

    TFx, TFy, TFz = total_length.(Face, topology, N, H)
    TCx, TCy, TCz = total_length.(Cell, topology, N, H)

    # Include halo points in coordinate arrays
    xF = range(xF₋, xF₊; length = TFx)
    yF = range(yF₋, yF₊; length = TFy)
    zF = range(zF₋, zF₊; length = TFz)

    xC = range(xC₋, xC₊; length = TCx)
    yC = range(yC₋, yC₊; length = TCy)
    zC = range(zC₋, zC₊; length = TCz)

    # Offset.
    xC = OffsetArray(xC, -Hx)
    yC = OffsetArray(yC, -Hy)
    zC = OffsetArray(zC, -Hz)

    xF = OffsetArray(xF, -Hx)
    yF = OffsetArray(yF, -Hy)
    zF = OffsetArray(zF, -Hz)

    return RegularCartesianGrid{FT, TX, TY, TZ, typeof(xC)}(
        Nx, Ny, Nz, Hx, Hy, Hz, Lx, Ly, Lz, Δx, Δy, Δz, xC, yC, zC, xF, yF, zF)
end

short_show(grid::RegularCartesianGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} =
    "RegularCartesianGrid{$FT, $TX, $TY, $TZ}(Nx=$(grid.Nx), Ny=$(grid.Ny), Nz=$(grid.Nz))"

function domain_string(grid)
    xₗ, xᵣ = x_domain(grid)
    yₗ, yᵣ = y_domain(grid)
    zₗ, zᵣ = z_domain(grid)
    return "x ∈ [$xₗ, $xᵣ], y ∈ [$yₗ, $yᵣ], z ∈ [$zₗ, $zᵣ]"
end

function show(io::IO, g::RegularCartesianGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ}
    print(io, "RegularCartesianGrid{$FT, $TX, $TY, $TZ}\n",
              "                   domain: $(domain_string(g))\n",
              "                 topology: ", (TX, TY, TZ), '\n',
              "  resolution (Nx, Ny, Nz): ", (g.Nx, g.Ny, g.Nz), '\n',
              "   halo size (Hx, Hy, Hz): ", (g.Hx, g.Hy, g.Hz), '\n',
              "grid spacing (Δx, Δy, Δz): ", (g.Δx, g.Δy, g.Δz))
end

"""
    with_halo(new_halo, old_grid::RegularCartesianGrid)

Returns a new `RegularCartesianGrid` with the same properties as
`old_grid` but with halos set to `new_halo`.

Note that in contrast to the constructor for `RegularCartesianGrid`,
`new_halo` is expected to be a 3-`Tuple` by `with_halo`. The elements
of `new_halo` corresponding to `Flat` directions are removed (and are
therefore ignored) prior to constructing the new `RegularCartesianGrid`.
"""
function with_halo(new_halo, old_grid::RegularCartesianGrid)

    FT = eltype(old_grid)
    Nx, Ny, Nz = size = (old_grid.Nx, old_grid.Ny, old_grid.Nz)
    topo = topology(old_grid)

    x = x_domain(old_grid)
    y = y_domain(old_grid)
    z = z_domain(old_grid)

    # Remove elements of size and new_halo in Flat directions as expected by grid
    # constructor
    size = pop_flat_elements(size, topo)
    new_halo = pop_flat_elements(new_halo, topo)

    new_grid = RegularCartesianGrid(eltype(old_grid); size=size, x=x, y=y, z=z,
                                    topology=topo, halo=new_halo)

    return new_grid
end


text:
"""
Internal function that wraps Requests.get
"""
function _get_request(uri::String)
    resp = HTTP.get(uri)
    if resp.status != 200
        #desc = STATUS_CODES[status]
        error("Expected status code 200 but received $resp.status")
    end
    return resp
end

"""
Internal function that parses the response
"""
function _parse_response(data, datatype::String)
    if datatype == "csv"
        return readdlm(data.body, ',', header=true)
    elseif datatype == "json"
        return JSON.Parser.parse(String(data.body))
    end
end

"""
Internal function that helps forms the request uri
"""
function _form_uri_tail(outputsize::String, datatype::String)
    o = "outputsize=$outputsize"
    d = "datatype=$datatype"
    a = "apikey=" * ENV["ALPHA_VANTAGE_API_KEY"]
    "&" * o * "&" * d * "&" * a
end

"""
Internal function that helps forms the request uri
"""
function _form_uri_head(func::String)
    uri = "$(alphavantage_api)query?"
    f = "function="* uppercase(func)
    uri * f
end

function _parse_params(params)
    if isempty(params)
        return ""
    else
        args = keys(params)
        values = collect(params)
        return mapreduce(i-> "$(args[i])=$(params[i])&", *, 1:length(params))
    end
end