text:
global max_spectra_to_use = 200
if max_spectra_to_use < 200
   @warn "param.in setting max_spectra_to_use to " * string(max_spectra_to_use)
end

global tophap_ccf_mask_scale_factor=1.6

global fits_target_str
if fits_target_str == "Solar"
   global linelist_for_ccf_filename = "G2.espresso.mas"
   hostname = gethostname()
   if occursin("aci.ics.psu.edu",hostname)
      global ancilary_solar_data_path = "/gpfs/group/ebf11/default/ebf11/neid_solar"
   elseif occursin("nuc8",hostname)  # Eric's home machine :)
      global ancilary_solar_data_path = "/home/eford/Data/SolarSpectra/NEID_solar/"
   end
   global ccf_mid_velocity = 0
   global bjd_first_good = 2458745.1296134139
   global bjd_last_good = 2458745.283
   global df_files
   global df_files_use = df_files |>
      @filter( _.target == fits_target_str ) |>
      @filter(bjd_first_good <= _.bjd < bjd_last_good) |>
      @take(max_spectra_to_use) |>
      DataFrame
end

if fits_target_str == "101501"
   global espresso_mask_filename = "G9.espresso.mas"
   global psu_mask3_filename = "HD101501q90f95n8e=false_mask.csv"
   global psu_mask4_filename = "HD101501q45f95n8e=false_mask.csv"

   global ccf_mid_velocity = -5.0e3
   global df_files
   global df_files_use = df_files |>
      @filter( _.target == fits_target_str ) |>
      @take(max_spectra_to_use) |>
      DataFrame
end

if fits_target_str == "10700"
   global espresso_mask_filename = "G9.espresso.mas"
   global psu_mask3_filename = "HD10700q90f80n7e=false_mask.csv"
   global psu_mask4_filename = "HD10700q50f80n7e=false_mask.csv"
   global ccf_mid_velocity = -16640.0
   global df_files
   global df_files_use = df_files |>
      @filter( _.target == fits_target_str ) |>
      @take(max_spectra_to_use) |>
      DataFrame
end

if fits_target_str == "34411"  #G1?
   global espresso_mask_filename = "G2.espresso.mas"
   global psu_mask3_filename = "HD34411q90f80n0e=false_mask.csv"
   global psu_mask4_filename = "HD34411q65f80n0e=false_mask.csv"
   global ccf_mid_velocity = 66500.0
   global df_files
   global df_files_use = df_files |>
      @filter( _.target == fits_target_str ) |>
      @take(max_spectra_to_use) |>
      DataFrame
end

if fits_target_str == "26965" # K0?
   global espresso_mask_filename = "G9.espresso.mas"
   global psu_mask3_filename = "HD26965q90f95n3e=false_mask.csv"
   global psu_mask4_filename = "HD26965q55f95n3e=false_mask.csv"
   global ccf_mid_velocity = -40320.0
   global df_files
   global df_files_use = df_files |>
      @filter( _.target == fits_target_str ) |>
      @take(max_spectra_to_use) |>
      DataFrame
end


text:
# multinomial idea.

# x is the input.
# function evalsquaredRKHS(   x::Vector{T},
#                             θ::SquaredFunctionsKernelType,
#                             X::Vector{Vector{T}},
#                             c::Vector{T},
#                             𝓐)::T where T <: Real
#
#     # get kernel components.
#     N = length(X)
#     @assert length(c) == N
#
#     f_x_components = collect( c[i]*evalkernel(x, X[i], θ.base_params) for i = 1:N)
#
#     # sum.
#     out = @distributed (+) for u in 𝓐
#         Combinatorics.multinomial(u...)*prod(f_x_components.^u)
#     end
#
#     return out
# end

function evalsquaredRKHS(   x::Vector{T},
                            θ::SquaredFunctionsKernelType,
                            X::Vector{Vector{T}},
                            c::Vector{T},
                            𝓐)::T where T <: Real

    # get kernel components.
    N = length(X)
    @assert length(c) == N

    f_x_components = collect( c[i]*evalkernel(x, X[i], θ.base_params) for i = 1:N)

    # sum.
    out = sum( Combinatorics.multinomial(u...)*prod(f_x_components.^u) for u in 𝓐 )

    return out
end

# assume y is strictly positive.
# to query: collect( evalsquaredRKHS(x, θp, η.X, η.c, 𝓐) for x in xq )
function fitsquarefunc( X::Vector{Vector{T}},
                        y::Vector{T},
                        θ::KT,
                        σ²::T) where {T,KT}
    #
    N = length(X)

    # fit suare root of y.
    y_sqrt = sqrt.(y)
    η = RKHSProblemType( zeros(T,N),
                         X,
                         θ,
                         σ²)
    fitRKHS!(η, y_sqrt)

    # prepare multinomial indices.
    θp = SquaredFunctionsKernelType(θ)
    𝓐 = SymTensorTools.getsubscriptarray(2,N)

    return η, 𝓐, θp
end

ext: jl

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


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

using Markdown
using InteractiveUtils

# ╔═╡ 710bd51c-18d7-11eb-2a05-c7f9f8f34324
begin
	using HypergeometricFunctions
	using LinearAlgebra
	using SparseArrays
	using PyPlot
	using BenchmarkTools
end

# ╔═╡ 51590b9c-3c39-11eb-0a39-6bb9f81b7d68
md"#### Fractional Fourier-Kravchuk evolution

Date: December 11, 2020
"

# ╔═╡ 02f4dd56-18d9-11eb-2d38-379eca82ec3a
begin
	rc("text", usetex = true)
	ion()
end

# ╔═╡ be283892-18d0-11eb-0154-758e466c7a40
# Global constants are fine, I suppose
begin
	const p = 2.0 # For the symmetric Kravchuk polynomial
	const j = 16 # Angular momentum number
	const N = 2*j # Dimension
end;

# ╔═╡ ffe361c6-18d0-11eb-2aa6-af741f44c0f9
# There's a problem in the _₂F₁ who have a NaN value when n = q = N, the bottom right corner of the array
function k_p(n,q)
	if n == q == N # Ad hoc condition added for completeness of the calculation
		return 0
	else
		return HypergeometricFunctions._₂F₁general2(-n, -q, -N, p)
	end
end;

# ╔═╡ a4b42958-3c10-11eb-3efa-b505bbd7df59
begin
	A(n, j) = ((-1)^n / 2^j)
	B(n, j) = binomial(2 * j, n)
end;

# ╔═╡ d3a38eca-3c10-11eb-20a8-6f0267436c21
k_f(n,q) = A(n, j) * sqrt(B(n, j) * B(j + q, j)) * k_p(n, j + q);

# ╔═╡ d86e82a6-39e1-11eb-3e5d-05a1af7f904e
K = reshape([k_f(l,m) for m in -j:j for l in 0:N], (N+1,N+1)); # Kravchuk matrix

# ╔═╡ be9d003e-3c14-11eb-3f4b-4771e5c06bec
id(α,n) = exp(-1im * pi * α * n / 2.0); # Unit roots

# ╔═╡ c9a7e70c-3c17-11eb-3038-7327c598667a
idT = [[id(α,n) for n in 0:N] for α in range(0, stop = 2, length = 315)]; # Vectors for the unit roots of dimensions N×(partition of 2π)

# ╔═╡ c635d95a-3c18-11eb-37b1-5dcccc976e80
dia(α) = Diagonal(idT[α]); # Diagonal matrix of the unit roots vectors, are bounded for α ∈ [1, ..., length of the 2π partition]

# ╔═╡ d8c7b53a-3c19-11eb-2890-759ed15872ff
begin
	pint(q) = sparsevec(Dict(q + j + 1 => 1), N+1); # Vector of zeros with a unit entry at q ∈ (-j, ..., j)
	eigK(n) = K[n,:]; # Eigenfunctions of the Kravchuk oscillator
	disK(z) = [k_f(0,m+z) for m in -j:j]; # Displaced n = 0 eigenstates of the Kravchuk oscillator
end

# ╔═╡ aa913f54-3c34-11eb-267a-6dc645a56600
vK = transpose(K) * normalize(disK(8)); # Transformation of the initial condition to the Kravchuk space

# ╔═╡ 35f1846c-3c2d-11eb-20fa-7954ca2d82a9
vFrT(α) = K * dia(α) * vK; # Propagation function in Kravchuk space and returning to configuration space

# ╔═╡ 28a9a1de-3c20-11eb-09b4-61e0b5d3c9c5
FrK = [abs2.(vFrT(α)) for α in 1:1:length(idT)]; # Propagation

# ╔═╡ 3456ff94-3c35-11eb-0000-ef9e984797bb
# Sadly the resulting matrix is 315 × N, so we need to reaccomodate to be N × 315
begin
	FK = zeros(Float64,(N+1,315))
	for l in 1:N+1
		for m in 1:315
			FK[l,m] = FrK[m][l]
		end
	end
end

# ╔═╡ 91da32e4-3c35-11eb-11ec-3fd798a5d2b2
begin
	fig1 = figure(figsize = (10,5))
	ax1 = gca()
	ax1.imshow(FK, cmap = "hot")
	ax1.set_aspect("auto")

	tight_layout()
end

# ╔═╡ Cell order:
# ╟─51590b9c-3c39-11eb-0a39-6bb9f81b7d68
# ╠═710bd51c-18d7-11eb-2a05-c7f9f8f34324
# ╠═02f4dd56-18d9-11eb-2d38-379eca82ec3a
# ╠═be283892-18d0-11eb-0154-758e466c7a40
# ╠═ffe361c6-18d0-11eb-2aa6-af741f44c0f9
# ╠═a4b42958-3c10-11eb-3efa-b505bbd7df59
# ╠═d3a38eca-3c10-11eb-20a8-6f0267436c21
# ╠═d86e82a6-39e1-11eb-3e5d-05a1af7f904e
# ╠═be9d003e-3c14-11eb-3f4b-4771e5c06bec
# ╠═c9a7e70c-3c17-11eb-3038-7327c598667a
# ╠═c635d95a-3c18-11eb-37b1-5dcccc976e80
# ╠═d8c7b53a-3c19-11eb-2890-759ed15872ff
# ╠═aa913f54-3c34-11eb-267a-6dc645a56600
# ╠═35f1846c-3c2d-11eb-20fa-7954ca2d82a9
# ╠═28a9a1de-3c20-11eb-09b4-61e0b5d3c9c5
# ╠═3456ff94-3c35-11eb-0000-ef9e984797bb
# ╠═91da32e4-3c35-11eb-11ec-3fd798a5d2b2


text:
# extending Base.in to Symbolics.Num types (useful for HyperGraphs of Num type)
Base.in(x::Num, collection::AbstractVector) = any(isequal(x, y) for y in collection)

# returns the most specialised type
get_lowest_subtype(t1::Type, t2::Type) = t1 <: t2 ? t1 : t2
get_lowest_subtype(t1::Type, t2::Type, tn...) = get_lowest_subtype(get_lowest_subtype(t1, t2), tn...)

# helper function to unwrap Num types
get_value(n::Num) = Symbolics.value(n)
get_value(eq::Equation) = Equation(get_value(eq.lhs), get_value(eq.rhs))
get_value(eqs::Vector{Equation}) = get_value.(eqs)

# temporary addition in order to buffer out changes in Symbolics; allows to fall
# back on a default case while the code is still being updated
get_value(x::T) where {T<:Symbolic} = x

# helper that returns the name of a symbolic variable, for a few different types
get_sym_name(s::Term) = s.f.name
get_sym_name(s::Sym) = s.name
get_sym_name(s::Num) = get_sym_name(get_value(s))

# helpers that abstract away the Symbolics syntax
function _check_symvar(var::Num)
    if isa(get_value(var), Sym)
        true
    else
        vartype = typeof(get_value(var)).name.wrapper
        error("$var of type $vartype; expected SymbolicUtils.Sym")
    end
end
_make_subname(symname::Symbol, idx::Int) = Symbol(symname, join(Symbolics.map_subscripts.(idx), "ˏ"))
genvar(symname::Symbol) = first(@variables $symname)
function genvar(symname::Symbol, idx::Int)
    (name = _make_subname(symname, idx); first(@variables $name))
end
function genvar(symname::Symbol, var::Num)
    _check_symvar(var) && first(@variables $symname($var))
end
function genvar(symname::Symbol, idx::Int, var::Num)
    _check_symvar(var) && (name = _make_subname(symname, idx); first(@variables $name($var)))
end

# helper acting as a gateway to genvar. this might later check whether the variable already exists or not; if it does, it would then just evaluate it and not generate it again. may be renamed to use_var when that behaviour is implemented
make_var(symname::Symbol, args...) = genvar(symname, args...)

# helper function that returns the symbolic variable given as input with the given subscript
function subscript_var(var::Symbolic, sub::Int; iv = Symbolics.variable(:t))
    symname = get_sym_name(var)
    if var isa Term
        Symbolics.variable(symname, sub, T = Symbolics.FnType)(iv)
    elseif var isa Sym
        Symbolics.variable(symname, sub)
    end
end
subscript_var(v::Num, sub; iv = Symbolics.variable(:t)) = Num(subscript_var(get_value(v), sub, iv = iv))
subscript_vars(vars::Vector{T}, sub; iv = Symbolics.variable(:t)) where {T<:SymTypes} = [subscript_var(var, sub, iv = iv) for var in vars]

# wrap SymbolicUtils variables in chemical hyperedges and hypergraphs into Symbolics.Num
# useful when converting ModelingToolkit types to chemical hypergraphs and hyperedges and still use the Simulacrum functions that are defined on Num
function Base.convert(::Type{ChemicalHyperEdge{Num}}, che::ChemicalHyperEdge{T}) where {T<:Symbolic}
    ChemicalHyperEdge(SpeciesSet(Num.(objects(che.src)), src_stoich(che)), SpeciesSet(Num.(objects(che.tgt)), tgt_stoich(che)), weight(che))
end
function Base.convert(::Type{ChemicalHyperGraph{Num}}, chx::ChemicalHyperGraph{T}) where {T<:Symbolic}
    ChemicalHyperGraph{Num}(Num.(vertices(chx)), convert.(ChemicalHyperEdge{Num}, hyperedges(chx)))
end

# convert directed edge indices to oriented hyperedges (one source, several targets from aggregated edges)
function edge_idx_to_hyperedge_idx(idx::Vector{Tuple{Int, Int}})

    # group ordinary edge indices according to their source vertex;
    # edge indices in each group have the same source vertex but different target vertex
    groups = [filter(v -> first(v) == i, idx) for i in 1:maximum(first.(idx))]

    # make hyperedge indices by concatenating all targets within each group
    elist = Vector{Tuple{Vector{Int}, Vector{Int}}}(undef, length(groups))
    @inbounds for (i, group) in enumerate(groups)
        elist[i] = ([first(group[1])], last.(group))
    end
    elist

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 96eb9ae0-30bc-11eb-1e22-890708244624
begin
	using RDatasets
	using MLJ
	using Flux
	import Random.seed!; seed!(123)
	using Plots
end

# ╔═╡ 0d3bf19e-3102-11eb-2661-73574bc6552b
md"### O problema

url: https://archive.ics.uci.edu/ml/datasets/Iris/

Data Set Information:

This is perhaps the best known database to be found in the pattern recognition literature. Fisher's paper is a classic in the field and is referenced frequently to this day. (See Duda & Hart, for example.) The data set contains 3 classes of 50 instances each, where each class refers to a type of iris plant. One class is linearly separable from the other 2; the latter are NOT linearly separable from each other.

Predicted attribute: class of iris plant.

This is an exceedingly simple domain.

This data differs from the data presented in Fishers article (identified by Steve Chadwick, spchadwick '@' espeedaz.net ). The 35th sample should be: 4.9,3.1,1.5,0.2,Iris-setosa where the error is in the fourth feature. The 38th sample: 4.9,3.6,1.4,0.1,Iris-setosa where the errors are in the second and third features.

Attribute Information:

1. sepal length in cm
2. sepal width in cm
3. petal length in cm
4. petal width in cm
5. class:
-- Iris Setosa
-- Iris Versicolour
-- Iris Virginica"

# ╔═╡ e72f92ae-30c4-11eb-1940-6d50060c106d
iris = RDatasets.dataset("datasets", "iris");

# ╔═╡ 3dadcf40-3114-11eb-39db-7108c541c711
head(iris)

# ╔═╡ 8a2e0720-3102-11eb-3f44-7d2bc065e1ba
describe(iris,:mean,:std,:min,:max,:eltype)

# ╔═╡ f0538220-30c4-11eb-37ad-a125378947fe
y, X = unpack(iris, ==(:Species), colname -> true, rng=123);

# ╔═╡ 53d0da90-3116-11eb-06b3-a3c6ebc1a281
md"É interessante verificarmos as dimensõs dos nossos dados e também o número de classes. MLJFlux já identifica para você e ajusta a entrada e saída da rede. Entretanto, a arquitetura é simplória. Veremos mais a frente como podemos fazer a nossa arquitetura."

# ╔═╡ b00ae7b0-30e9-11eb-3649-41d9c676166c
size(X),unique(y)

# ╔═╡ fbad2992-3101-11eb-01c3-5f4c0a471575
md"### Criando o classificador"

# ╔═╡ 87101960-30c2-11eb-367f-dfeab5e2cffb
begin
	@load NeuralNetworkClassifier
	clf = NeuralNetworkClassifier()
end

# ╔═╡ a65f7890-30c4-11eb-2c1f-ed94a0ec1243
mach = machine(clf, X, y)

# ╔═╡ ea5995c2-3101-11eb-36a3-25db3da1790b
md"### Particionando os dados

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."

# ╔═╡ 6d1c16a0-30c5-11eb-1e52-6f98c2648e77
train, test = partition(eachindex(y), 0.9, shuffle=true, rng=1234)

# ╔═╡ d572a200-3101-11eb-2190-4ff72e2dc603
md"### Treinando a Rede

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."

# ╔═╡ b827b20e-30c2-11eb-2229-f703a237dcc4
fit!(mach,rows=train)

# ╔═╡ bbef79c0-3115-11eb-3948-f3c4bf0be576
md" Aqui conseguimos ver a arquitetura da rede padrão do MLJFlux. Podemos alterar se quisermos."

# ╔═╡ 365eb290-30e2-11eb-3665-a3d1951a9fcc
fitted_params(mach).chain

# ╔═╡ a7500660-3101-11eb-1430-a5021429fdec
md"### Avaliando a qualidade no conjunto de teste

Aqui coloco algumas métricas já conhecidas de classificação para avaliarmos no conjunto de teste."

# ╔═╡ 87dd9c60-30df-11eb-0dd7-1976c80dd295
ŷ = predict_mode(mach, rows=test)

# ╔═╡ f92ad200-30e1-11eb-0168-35b7013a8ba5
cm = confusion_matrix(ŷ,y[test])

# ╔═╡ 05c4ae4e-30e2-11eb-0018-37d1d82e0cc4
accuracy(cm)

# ╔═╡ 0735f190-30e2-11eb-1b10-974050eb025d
matthews_correlation(cm)

# ╔═╡ 24a34c80-30f3-11eb-3224-9b3dbd304190
md"## Analisando a qualidade da arquitetura

Uma coisa comum na avaliação de redes é a análise da qualidade por época. Os resultados aqui são da predição a cada época. "

# ╔═╡ 384339e0-31ba-11eb-0d12-fbf707d43b44
md"#### Analisando a sensibilidade com relação as épocas

É importante descobrir o número ideal de épocas. Para isso, devemos registrar o erro por época."

# ╔═╡ 60da1820-30e2-11eb-3992-73afe463dc94
begin
	re = range(clf, :epochs, lower=1, upper=550, scale=:log10)
	curve = learning_curve(clf,X,y,
						   range=re,					                             resampling=Holdout(fraction_train=0.9,
			         shuffle=true,rng=123),#CV(nfolds=10,  shuffle=true, rng=123))
		             measure = cross_entropy)

end

# ╔═╡ 4271ca40-31b4-11eb-2d56-2b521a23f1bd
plot(curve.parameter_values,
		   curve.measurements,
		   xlab=curve.parameter_name,
		   xscale=curve.parameter_scale,
		   ylab = "Cross Entropy",
		   label = "Validation Error")

# ╔═╡ 484c9a70-31ba-11eb-0b28-7bb96e703c3a
md"#### Analisando a sensibilidade com relação à taxa de aprendizado

Sabemos que a solução do método de aprendizado depende da taxa de aprendizado.
"

# ╔═╡ 04b46700-31b7-11eb-1702-a7af6b66681c
begin
	rl = range(clf, :lambda, lower=0.001, upper=1)
	curvel = learning_curve(clf,X,y,
						   range=rl,					                             resampling=Holdout(fraction_train=0.9,
			         shuffle=true,rng=123),#CV(nfolds=10,  shuffle=true, rng=123))
		             measure = cross_entropy)

end

# ╔═╡ 1780b820-31b7-11eb-255d-b99cfb673820
plot(curvel.parameter_values,
		   curvel.measurements,
		   xlab=curvel.parameter_name,
		   xscale=curvel.parameter_scale,
		   ylab = "Cross Entropy",
		   label = "Validation Error")

# ╔═╡ 30647c00-31b7-11eb-3eed-3999ee53a4a9
clf.lambda

# ╔═╡ 65213270-31b8-11eb-3682-e51a60edd85d
clf.epochs

# ╔═╡ 92346140-31b5-11eb-0317-8f9a5589b0df
curve

# ╔═╡ c5d00cc0-31b5-11eb-0abf-a5eed9d49f11
clf.epochs=300

# ╔═╡ 5afe85c0-31ba-11eb-309d-1f526cf67d66
md"### Realizando uma validação cruzada

O objetivo aqui é estimar com menor viés possível a qualidade da rede para a tarefa"

# ╔═╡ 8d6ff0d0-31b4-11eb-1860-ada5b3b6b84c
history= evaluate!(
    mach;
    resampling=CV(nfolds=10,  shuffle=true, rng=123),#Holdout(fraction_train=0.7, shuffle=true, rng=123),
    operation=predict_mode,
    measure=[accuracy, misclassification_rate],
    verbosity = 3
)

# ╔═╡ 102dfbb0-31b6-11eb-2019-c520f466d6a6
history

# ╔═╡ 1c0a34c0-30ef-11eb-26b0-b38e8f85eca6
md"### Criando a sua própria arquitetura de rede"

# ╔═╡ 17e44710-3116-11eb-0f52-2f6657afcead
md" Um ponto interessant é podermos modificar nossa arquitetura da rede e verificar se o desempenho melhora ou piora. No MLJFlux é preciso criar uma estrutura que será passada por argumento e nós sobrescrevemos a função build."

# ╔═╡ 6c58bf70-30ed-11eb-2fe4-111b91c5fdd2
begin
	mutable struct MyNetwork <: MLJFlux.Builder
		n1 :: Int
		n2 :: Int
	end

	function MLJFlux.build(nn::MyNetwork, n_in, n_out)
		return Chain(Dense(n_in, nn.n1,relu),
			   Dense(nn.n1, nn.n2,relu),
			   Dense(nn.n2, n_out))
	end
end

# ╔═╡ 3503f1a0-30ef-11eb-0dc0-1bbcdf6478b7
myclf = NeuralNetworkClassifier(builder=MyNetwork(10,10))

# ╔═╡ 71a84080-31ac-11eb-034e-69fde3328410
mymach = machine(myclf, X, y)

# ╔═╡ 06495780-31ba-11eb-3e53-9732577482ef
md"### Analisando a sensibilidade com relação as épocas"

# ╔═╡ 813e9d90-30f4-11eb-361a-c1a55c36f13c
begin
	re1 = range(myclf, :epochs, lower=1, upper=250, scale=:log10)
	curvee1 = learning_curve(myclf,X,y,
						   range=re1,					                             resampling=Holdout(fraction_train=0.9,
			         shuffle=true,rng=123),#CV(nfolds=10,  shuffle=true, rng=123))
		             measure = cross_entropy)
end

# ╔═╡ 27f2a5e0-31b4-11eb-1076-e5851c3482e1
plot(curvee1.parameter_values,
		   curvee1.measurements,
		   xlab=curvee1.parameter_name,
		   xscale=curvee1.parameter_scale,
		   ylab = "Cross Entropy",
	       label = "Error")

# ╔═╡ 15ef1440-31ba-11eb-0a5d-8d9578bf921f
md"### Realizando uma validação cruzada

O objetivo aqui é estimar com menor viés possível a qualidade da rede para a tarefa
"

# ╔═╡ ad02ba40-31b9-11eb-033e-2b06818b1dea
history2 = evaluate!(
    mymach;
    resampling=CV(nfolds=10,  shuffle=true, rng=123),#Holdout(fraction_train=0.7, shuffle=true, rng=123),
    operation=predict_mode,
    measure=[accuracy, misclassification_rate],
    verbosity = 3
)

# ╔═╡ ccafdbc0-31b9-11eb-2f03-07aaf5a96058
history2

# ╔═╡ f5af28a0-31b9-11eb-14e7-5794590639c4
md"### Um simples holdout"

# ╔═╡ 2e6caf2e-31ad-11eb-295a-29d5a6d93bc0
fit!(mymach,rows=train)

# ╔═╡ 9e552180-31ab-11eb-3713-275f29c2ce1d
ŷn = predict_mode(mymach, rows=test)

# ╔═╡ bc93e5ee-31ab-11eb-245e-ed426c077c4e
cmn = confusion_matrix(ŷn,y[test])

# ╔═╡ 4ebcee80-31ad-11eb-3d22-9b957c247948
accuracy(cmn)

# ╔═╡ Cell order:
# ╠═96eb9ae0-30bc-11eb-1e22-890708244624
# ╟─0d3bf19e-3102-11eb-2661-73574bc6552b
# ╠═e72f92ae-30c4-11eb-1940-6d50060c106d
# ╠═3dadcf40-3114-11eb-39db-7108c541c711
# ╠═8a2e0720-3102-11eb-3f44-7d2bc065e1ba
# ╠═f0538220-30c4-11eb-37ad-a125378947fe
# ╟─53d0da90-3116-11eb-06b3-a3c6ebc1a281
# ╠═b00ae7b0-30e9-11eb-3649-41d9c676166c
# ╟─fbad2992-3101-11eb-01c3-5f4c0a471575
# ╠═87101960-30c2-11eb-367f-dfeab5e2cffb
# ╠═a65f7890-30c4-11eb-2c1f-ed94a0ec1243
# ╟─ea5995c2-3101-11eb-36a3-25db3da1790b
# ╠═6d1c16a0-30c5-11eb-1e52-6f98c2648e77
# ╟─d572a200-3101-11eb-2190-4ff72e2dc603
# ╠═b827b20e-30c2-11eb-2229-f703a237dcc4
# ╟─bbef79c0-3115-11eb-3948-f3c4bf0be576
# ╠═365eb290-30e2-11eb-3665-a3d1951a9fcc
# ╟─a7500660-3101-11eb-1430-a5021429fdec
# ╠═87dd9c60-30df-11eb-0dd7-1976c80dd295
# ╠═f92ad200-30e1-11eb-0168-35b7013a8ba5
# ╠═05c4ae4e-30e2-11eb-0018-37d1d82e0cc4
# ╠═0735f190-30e2-11eb-1b10-974050eb025d
# ╠═24a34c80-30f3-11eb-3224-9b3dbd304190
# ╟─384339e0-31ba-11eb-0d12-fbf707d43b44
# ╠═60da1820-30e2-11eb-3992-73afe463dc94
# ╠═4271ca40-31b4-11eb-2d56-2b521a23f1bd
# ╟─484c9a70-31ba-11eb-0b28-7bb96e703c3a
# ╠═04b46700-31b7-11eb-1702-a7af6b66681c
# ╠═1780b820-31b7-11eb-255d-b99cfb673820
# ╠═30647c00-31b7-11eb-3eed-3999ee53a4a9
# ╠═65213270-31b8-11eb-3682-e51a60edd85d
# ╠═92346140-31b5-11eb-0317-8f9a5589b0df
# ╠═c5d00cc0-31b5-11eb-0abf-a5eed9d49f11
# ╟─5afe85c0-31ba-11eb-309d-1f526cf67d66
# ╠═8d6ff0d0-31b4-11eb-1860-ada5b3b6b84c
# ╠═102dfbb0-31b6-11eb-2019-c520f466d6a6
# ╟─1c0a34c0-30ef-11eb-26b0-b38e8f85eca6
# ╟─17e44710-3116-11eb-0f52-2f6657afcead
# ╠═6c58bf70-30ed-11eb-2fe4-111b91c5fdd2
# ╠═3503f1a0-30ef-11eb-0dc0-1bbcdf6478b7
# ╠═71a84080-31ac-11eb-034e-69fde3328410
# ╠═06495780-31ba-11eb-3e53-9732577482ef
# ╠═813e9d90-30f4-11eb-361a-c1a55c36f13c
# ╠═27f2a5e0-31b4-11eb-1076-e5851c3482e1
# ╟─15ef1440-31ba-11eb-0a5d-8d9578bf921f
# ╠═ad02ba40-31b9-11eb-033e-2b06818b1dea
# ╠═ccafdbc0-31b9-11eb-2f03-07aaf5a96058
# ╠═f5af28a0-31b9-11eb-14e7-5794590639c4
# ╠═2e6caf2e-31ad-11eb-295a-29d5a6d93bc0
# ╠═9e552180-31ab-11eb-3713-275f29c2ce1d
# ╠═bc93e5ee-31ab-11eb-245e-ed426c077c4e
# ╠═4ebcee80-31ad-11eb-3d22-9b957c247948

max_stars_repo_name: lalvim/machinelearning

max_issues_repo_name: lalvim/machinelearning

max_forks_repo_name: lalvim/machinelearning


text:
using PkgBenchmark
using Random
using MIToS.Utils
using MIToS.MSA

@benchgroup "Utils" begin
    include(joinpath("Utils", "GeneralUtils.jl"))
end

@benchgroup "MSA" begin
    include(joinpath("MSA", "Residues.jl"))
    include(joinpath("MSA", "Annotations.jl"))
end


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

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Testing

function JuliaFEM.get_mesh(::Type{Val{Symbol("two elements 1.0x0.5 with 0.1 gap in y direction")}})
    mesh = Mesh()
    add_node!(mesh, 1, [0.0, 0.0])
    add_node!(mesh, 2, [1.0, 0.0])
    add_node!(mesh, 3, [1.0, 0.5])
    add_node!(mesh, 4, [0.0, 0.5])
    add_node!(mesh, 5, [0.0, 0.6])
    add_node!(mesh, 6, [1.0, 0.6])
    add_node!(mesh, 7, [1.0, 1.1])
    add_node!(mesh, 8, [0.0, 1.1])
    add_element!(mesh, 1, :Quad4, [1, 2, 3, 4])
    add_element!(mesh, 2, :Quad4, [5, 6, 7, 8])
    add_element!(mesh, 3, :Seg2, [1, 2])
    add_element!(mesh, 4, :Seg2, [7, 8])
    add_element!(mesh, 5, :Seg2, [4, 3])
    add_element!(mesh, 6, :Seg2, [6, 5])
    add_element_to_element_set!(mesh, :LOWER, 1)
    add_element_to_element_set!(mesh, :UPPER, 2)
    add_element_to_element_set!(mesh, :LOWER_BOTTOM, 3)
    add_element_to_element_set!(mesh, :UPPER_TOP, 4)
    add_element_to_element_set!(mesh, :LOWER_TOP, 5)
    add_element_to_element_set!(mesh, :UPPER_BOTTOM, 6)
    return mesh
end

function JuliaFEM.get_model(::Type{Val{Symbol("two element contact")}})

    mesh = get_mesh("two elements 1.0x0.5 with 0.1 gap in y direction")

    upper = Problem(Elasticity, "UPPER", 2)
    upper.properties.formulation = :plane_stress
    upper.elements = create_elements(mesh, "UPPER")
    update!(upper.elements, "youngs modulus", 288.0)
    update!(upper.elements, "poissons ratio", 1/3)

    lower = Problem(Elasticity, "LOWER", 2)
    lower.properties.formulation = :plane_stress
    lower.elements = create_elements(mesh, "LOWER")
    update!(lower.elements, "youngs modulus", 288.0)
    update!(lower.elements, "poissons ratio", 1/3)

    bc_upper = Problem(Dirichlet, "UPPER_TOP", 2, "displacement")
    bc_upper.elements = create_elements(mesh, "UPPER_TOP")
    #update!(bc_upper.elements, "displacement 1", -17/90)
    #update!(bc_upper.elements, "displacement 1", -17/90)
    update!(bc_upper.elements, "displacement 1", -0.2)
    update!(bc_upper.elements, "displacement 2", -0.2)

    bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 2, "displacement")
    bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
    update!(bc_lower.elements, "displacement 1", 0.0)
    update!(bc_lower.elements, "displacement 2", 0.0)

    interface = Problem(Contact, "LOWER_TO_UPPER", 2, "displacement")
    interface.properties.dimension = 1
    interface_slave_elements = create_elements(mesh, "LOWER_TOP")
    interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
    update!(interface_slave_elements, "master elements", interface_master_elements)
    interface.elements = [interface_master_elements; interface_slave_elements]

    solver = Solver(Nonlinear)
    push!(solver, upper, lower, bc_upper, bc_lower, interface)
    return solver

end

@testset "test simple two element contact" begin
    solver = get_model("two element contact")
    solver()
    contact = solver["LOWER_TO_UPPER"]
    master = first(contact.elements)
    slave = last(contact.elements)
    u = master("displacement", [0.0], 0.0)
    la = slave("reaction force", [0.0], 0.0)
    info("u = $u, la = $la")
    @test isapprox(u, [-0.2, -0.15])
    @test isapprox(la, [0.0, -30.375])
    # FIXME
end


text:
function f(x::Int64)
    return x && true
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
import Pkg

Pkg.update()

root_directory = dirname(@__DIR__)

kernelabstractions = Pkg.PackageSpec(path = root_directory)

Pkg.develop(kernelabstractions)
if !(VERSION < v"1.6-")
    rockernels = Pkg.PackageSpec(path = joinpath(root_directory, "lib", "ROCKernels"))
    Pkg.develop(rockernels)

    cudakernels = Pkg.PackageSpec(path = joinpath(root_directory, "lib", "CUDAKernels"))
    Pkg.develop(cudakernels)
end
Pkg.build()
Pkg.precompile()


text:
using PowerModels
using Ipopt
using .PP2PM

function run_powermodels(json_path)
    pm = PP2PM.load_pm_from_json(json_path)
    result = PowerModels.run_dc_opf(pm, Ipopt.IpoptSolver(),
                                    setting = Dict("output" => Dict("branch_flows" => true)))
    return result
end

max_stars_repo_name: gisce/pandapower

max_issues_repo_name: gisce/pandapower

max_forks_repo_name: gisce/pandapower


text:
using CryptoCurrencies
import CryptoCurrencies: unit, name, BTC, LTC, XMR, ZEC

using Test

currencies = [BTC,LTC,XMR,ZEC]
units = [8,8,8,8]
names = ["Bitcoin","Litecoin","Monero","Zcash"]

for (ccy,u,n) in zip(currencies,units,names)
    @test unit(ccy) == u
    @test name(ccy) == n
end


text:
using MinMaxFilter
using Base.Test
using MAT

# #######
#
# 1 dimension
#
# Compare to matlab slowminmaxfilt_algo.m
#
# t = [1:1024] ./1024; d = sin(2*pi*4*t);
# [minval, maxval] = slowminmaxfilt_algo(d, 100)
# dlmwrite('minmax_output.txt',[minval; maxval], 'delimiter', '\t', 'precision', '%.12f')
#
# ######

matlab = readdlm(joinpath(dirname(@__FILE__), "data", "minmax_output.txt"),'\t')

t = (1:1024)./1024
d = sin(2*pi*4*t)
minval, maxval = minmax_filter(d, 100)

@test_approx_eq minval matlab[1,:]
@test_approx_eq maxval matlab[2,:]

# #######
#
# 2 dimension
#
# Compare to matlab minmaxfilt.m
#
# [X,Y,Z] = peaks(100);
# smax = minmaxfilt(Z, 11, 'max')
#
# ######

filen = matopen(joinpath(dirname(@__FILE__), "data", "2d_array.mat"))
A = read(filen, "Z")
A = convert(Array{AbstractFloat}, A)
close(filen)

filen = matopen(joinpath(dirname(@__FILE__), "data", "2d_array_max11.mat"))
max_matlab = read(filen, "smax")
max_matlab = convert(Array{AbstractFloat}, max_matlab)
close(filen)

filen = matopen(joinpath(dirname(@__FILE__), "data", "2d_array_min11.mat"))
min_matlab = read(filen, "smin")
min_matlab = convert(Array{AbstractFloat}, min_matlab)
close(filen)

minval, maxval = minmax_filter(A, [11, 11])

@test_approx_eq maxval max_matlab
@test_approx_eq minval min_matlab

# #######
#
# 3 dimension
#
# Compare to matlab minmaxfilt.m
#
# amax=minmaxfilt(image,5,'max','valid');
#
# ######

filen = matopen(joinpath(dirname(@__FILE__), "data", "3d_array.mat"))
A = read(filen, "image")
A = convert(Array{AbstractFloat}, A)
close(filen)

filen = matopen(joinpath(dirname(@__FILE__), "data", "3d_array_max5.mat"))
max_matlab = read(filen, "amax")
max_matlab = convert(Array{AbstractFloat}, max_matlab)
close(filen)

filen = matopen(joinpath(dirname(@__FILE__), "data", "3d_array_min5.mat"))
min_matlab = read(filen, "amin")
min_matlab = convert(Array{AbstractFloat}, min_matlab)
close(filen)

minval, maxval = minmax_filter(A, [5,5,5])

@test_approx_eq maxval max_matlab
@test_approx_eq minval min_matlab


text:
# This file was generated, do not modify it. # hide
ENV["GKSwstype"] = "nul" # hide

using SimilaritySearch

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# # `ClimateTools.jl` and `ClimatePlots.jl`
#
# [ClimateTools.jl](https://juliaclimate.github.io/ClimateTools.jl/dev/) is a collection of commonly-used tools aimed to ease the typical steps in (1) analyzing climate models outputs from netCDF files that follow CF-conventions and (2) creating climate scenarios. [ClimatePlots.jl](https://juliaclimate.github.io/ClimatePlots.jl/dev/) is the associated plotting library.

# + {"slideshow": {"slide_type": "subslide"}}
if false #set to true if you have not yet installed packages listed below
    using Pkg
    Pkg.add(PackageSpec(name="ClimateTools", rev="master"))

    Pkg.add("PyCall")
    ENV["PYTHON"]=""
    Pkg.build("PyCall")
    Pkg.add(PackageSpec(name="ClimatePlots", rev="master"))

    run(`wget http://esgf-data1.diasjp.net/thredds/fileServer/esg_dataroot/cmip5/output1/MIROC/MIROC5/piControl/day/atmos/day/r1i1p1/v20161012/tas/tas_day_MIROC5_piControl_r1i1p1_20000101-20091231.nc`)
    run(`mv tas_day_MIROC5_piControl_r1i1p1_20000101-20091231.nc ../inputs/`)
end

# + {"slideshow": {"slide_type": "-"}, "cell_type": "markdown"}
# _Note: `tas_day_MIROC5_piControl_*.nc` was downloaded [here](http://esgf-data1.diasjp.net/thredds/fileServer/esg_dataroot/cmip5/output1/MIROC/MIROC5/piControl/day/atmos/day/r1i1p1/v20161012/tas/tas_day_MIROC5_piControl_r1i1p1_20000101-20091231.nc) by selecting `piControl,day,tas,MIROC5` in [the search engine](https://esgf-node.llnl.gov/search/cmip5/)_
#
# ```
# project=CMIP5, model=MIROC5, Atmosphere and Ocean Research Institute (The University of Tokyo),
# experiment=pre-industrial control, time_frequency=day, modeling realm=atmos, ensemble=r1i1p1,
# Description: MIROC5 model output prepared for CMIP5 pre-industrial control
# ```

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## Get Meta-Data From File
#
# _Note: uncomment the final line to display the file meta data_

# + {"slideshow": {"slide_type": "-"}}
using ClimateTools, ClimatePlots

p="../inputs"
fil="$p/tas_day_MIROC5_piControl_r1i1p1_20000101-20091231.nc"
#fil="$p/clt_day_MIROC5_historical_r4i1p1_19500101-19591231.nc"
d=Dataset(fil);

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## Read & Plot A Variable

# + {"slideshow": {"slide_type": "fragment"}}
p1=joinpath(dirname(pathof(ClimateTools)),"../test/data")

#fil1="$p1/orog_fx_GFDL-ESM2G_historicalMisc_r0i0p0.nc"
fil1="$p1/sresa1b_ncar_ccsm3-example.nc"
model1 = load(fil1, "pr", data_units="mm")
contourf(model1, region = "Mollweide");

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## Extract Subset And Plot Region
#
# _Note : see `ClimatePlots.jl/src/maps_definition.jl`_

# + {"slideshow": {"slide_type": "fragment"}}
poly_reg = [[NaN -65 -80 -80 -65 -65];[NaN 42 42 52 52 42]]
model2 = load(fil, "tas", poly=poly_reg)
contourf(model2, region = "Quebec");

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## Inspect Meta Data

# + {"slideshow": {"slide_type": "fragment"}}
#d
#model1
model2
# -
# ## Read & Plot A Variable (bis)

p3="../outputs/nctiles-newfiles/interp"
tst=sum(occursin.("ETAN.nc",readdir(p3)))>0
if tst
    model3 = load("$p3/ETAN.nc", "ETAN")
    contourf(model3, region = "Mollweide")
end

max_stars_repo_name: abelolivelanna/GlobalOceanNotebooks

max_issues_repo_name: abelolivelanna/GlobalOceanNotebooks

max_forks_repo_name: abelolivelanna/GlobalOceanNotebooks


text:
import AbstractPlotting
using WGLMakie
using AbstractPlotting.Colors: color, reducec
using ImageTransformations

@testset "WGLMakie" begin
    function test_handler(session, req)
        scene = AbstractPlotting.scatter([1, 2, 3, 4], resolution=(500,500), color=:red)
        return DOM.div(scene)
    end
    testsession(test_handler) do app
        # Lets not be too porcelainy about this ...
        @test evaljs(app, js"document.querySelector('canvas').style.width") == "500px"
        @test evaljs(app, js"document.querySelector('canvas').style.height") == "500px"
        img = WGLMakie.session2image(app.session)
        ratio = evaljs(app, js"window.devicePixelRatio")
        @info "device pixel ration: $(ratio)"
        if ratio != 1
            img = ImageTransformations.imresize(img, (500, 500))
        end
        img_ref = AbstractPlotting.FileIO.load(joinpath(@__DIR__, "test_reference.png"))
        # AbstractPlotting.FileIO.save(joinpath(@__DIR__, "test_reference.png"), img)
        meancol = AbstractPlotting.mean(color.(img) .- color.(img_ref))
        # This is pretty high... But I don't know how to look into travis atm
        # And seems it's not suuper wrong, so I'll take it for now ;)
        @test (reducec(+, 0.0, meancol) / 3) <= 0.025
    end
end


text:
using Distributed, DataFrames
addprocs(48)

@everywhere  begin
    using Jedi
    using LinearAlgebra
    using Distributions
    using SharedArrays
    emat = 2 * (ones(4, 4) - Matrix{Float64}(I, 4, 4))
    Est(l) = 2*(3/4 * l - 5)
    N = 1000
    f0 = 20/2N
    f = fermi_fitness(f0=f0, E_Star=Est)
end

rho = [0, 0.1, 0.5, 1., 2]
E = SharedArray{Float64, 2}(length(rho), 1000)
L = SharedArray{Float64, 2}(length(rho), 1000)
RHO = SharedArray{Float64, 2}(length(rho), 1000)

@everywhere function run(rho, nu, l)
    pop = driver_trailer_l(N=1000, l_0=l, L=l)
    initiate_rand!(pop, 1)
    for i in 1:50000000
        Jedi.bp_substitution!(pop, emat, f)
        if rand() < rho/N
            driver_mutation!(pop)
        end
        if rand() < nu/N
            Jedi.l_substitution!(pop, emat, f)
        end
    end
    Gamma = get_energy(pop, emat)[1]
    l_arr = pop.l[1]
    return Gamma, l_arr
end

@sync @distributed for j in 1:1000
    for r in 1:length(rho)
        E[r, j], L[r, j] = run(rho[r], 0.01, 10)
        RHO[r, j] = rho[r]
    end
    println("Run $j done.")
end

df = DataFrame(gamma=[(E...)...], l=[(L...)...], rho=[(RHO...)...])
CSV.write("script2_results.csv", df)


text:
export Pool, MaxPool

"""
$(TYPEDEF)

Represents a pooling operation.

`p(x)` is shorthand for [`pool(x, p)`](@ref) when `p` is an instance of `Pool`.

## Fields:
$(FIELDS)
"""
struct Pool{N} <: Layer
    strides::NTuple{N, Int}
    pooling_function::Function
end

function Base.show(io::IO, p::Pool)
    (_, stride_height, stride_width, _) = p.strides
    function_display_name = Dict(
        MIPVerify.maximum => "max",
    )
    print(io,
        "$(function_display_name[p.pooling_function]) pooling with a $(stride_height)x$(stride_width) filter and a stride of ($stride_height, $stride_width)"
    )
end

Base.hash(a::Pool, h::UInt) = hash(a.strides, hash(string(a.pooling_function), hash(:Pool, h)))

"""
$(SIGNATURES)

Convenience function to create a [`Pool`](@ref) struct for max-pooling.
"""
function MaxPool(strides::NTuple{N, Int}) where {N}
    Pool(strides, MIPVerify.maximum)
end

"""
$(SIGNATURES)

For pooling operations on an array where a given element in the output array
corresponds to equal-sized blocks in the input array, returns (for a given
dimension) the index range in the input array corresponding to a particular
index `output_index` in the output array.

Returns an empty array if the `output_index` does not correspond to any input
indices.

# Arguments
* `stride::Integer`: the size of the operating blocks along the active
     dimension.

"""
function getsliceindex(input_array_size::Integer, stride::Integer, output_index::Integer)::Array{Int, 1}
    parent_start_index = (output_index-1)*stride+1
    parent_end_index = min((output_index)*stride, input_array_size)
    if parent_start_index > parent_end_index
        return []
    else
        return parent_start_index:parent_end_index
    end
end

"""
$(SIGNATURES)

For pooling operations on an array, returns a view of the parent array
corresponding to the `output_index` in the output array.
"""
function getpoolview(input_array::AbstractArray{T, N}, strides::NTuple{N, Int}, output_index::NTuple{N, Int})::SubArray{T, N} where {T, N}
    it = zip(size(input_array), strides, output_index)
    input_index_range = map(x -> getsliceindex(x...), it)
    return view(input_array, input_index_range...)
end

"""
$(SIGNATURES)

For pooling operations on an array, returns the expected size of the output
array.
"""
function getoutputsize(input_array::AbstractArray{T, N}, strides::NTuple{N, Int})::NTuple{N, Int} where {T, N}
    output_size = ((x, y) -> round(Int, x/y, RoundUp)).(size(input_array), strides)
    return output_size
end

"""
$(SIGNATURES)

Returns output from applying `f` to subarrays of `input_array`, with the windows
determined by the `strides`.
"""
function poolmap(f::Function, input_array::AbstractArray{T, N}, strides::NTuple{N, Int}) where {T, N}
    output_size = getoutputsize(input_array, strides)
    output_indices = collect(CartesianIndices(output_size))
    return ((I) -> f(getpoolview(input_array, strides, I.I))).(output_indices)
end

"""
$(SIGNATURES)

Computes the result of applying the pooling function `params.pooling_function` to
non-overlapping cells of `input` with sizes specified in `params.strides`.
"""
function pool(
    input::AbstractArray{T, N},
    params::Pool{N}) where {T<:JuMPReal, N}
    poolmap(params.pooling_function, input, params.strides)
end

(p::Pool)(x::Array{<:Real}) = MIPVerify.pool(x, p)
(p::Pool)(x::Array{<:JuMPLinearType}) = (Memento.info(MIPVerify.LOGGER, "Specifying $p ... "); MIPVerify.pool(x, p))


text:
module PrettyPrint
export pprint, pprintln, pformat, is_simple, is_simple_t, is_atom_t, is_static_t, pp_impl, pp_impl_dynamic, pp_impl_static, MaxIndentExpected

struct PPPair{FieldName, FieldContent}
    name :: FieldName
    sep :: String
    content :: FieldContent
end

struct PPMetaSymbol{S}
    s :: S
end

Base.show(io::IO, p::PPMetaSymbol) = Base.print(io, p.s)

is_simple_t(a) = isprimitivetype(a)
is_simple_t(::Type{<:Number}) = true
is_simple_t(::Type{<:AbstractString}) = true
is_simple_t(::Type{<:AbstractChar}) = true
is_simple_t(::Type{<:Symbol}) = true
is_simple_t(::Type{<:Complex}) = true
is_simple_t(::Type{<:PPMetaSymbol}) = true
is_simple_t(::Type{Nothing}) = true
is_simple_t(::Type{PPPair{A1, A2}}) where {A1, A2} = is_simple_t(A1) && is_simple_t(A2)
is_simple(::T) where T = is_simple_t(T)

is_atom_t(_) = false
is_atom_t(::Type{<:Number}) = true
is_atom_t(::Type{<:AbstractString}) = true
is_atom_t(::Type{<:AbstractChar}) = true
is_atom_t(::Type{<:Symbol}) = true
is_atom_t(::Type{<:Complex}) = true
is_atom_t(::Type{<:PPMetaSymbol}) = true
is_atom_t(::Type{Nothing}) = true
is_atom(::T) where T = is_atom_t(T)

is_static_t(a) = isprimitivetype(a)
is_static_t(::Type{String}) = true
is_static_t(::Type{Symbol}) = true
is_static_t(::Type{<:Complex}) = true

const Indentation = Int
const MaxIndentExpected = Ref{Int}(30)

print_repeat(n::Int, io::IO, e) = for i = 1:n
    print(io, e)
end

function pprint_for_seq(io, left::String, right::String, seq :: Union{Vector, Tuple}, indent::Int)::Indentation
    print(io, left)
    n = length(seq)
    is_empty = n === 0
    if is_empty
        print(io, right)
        return indent + length(left) + length(right)
    end

    if all(is_simple(e) for e in seq)
        indent += length(left)
        indent′ = indent
        max_indent = MaxIndentExpected[]
        indent′ = pp_impl(io, seq[1], indent)
        if n == 1
            print(io, ",")
            indent′ += 2
        end
        for i in 2:n
            e = seq[i]
            print(io, ", ")

            if indent′ > max_indent
                print(io, '\n')
                indent′ = indent
                print_repeat(indent, io, ' ')
            else
                indent′ += 2
            end

            indent′ = pp_impl(io, e, indent′)
        end

        print(io, right)
        indent′ += length(right)

        return indent′
    end

    print(io, '\n')
    indent′ = indent + 2
    for e in seq
        print_repeat(indent′, io, ' ')
        pp_impl(io, e, indent′)
        print(io, ',')
        print(io, '\n')
    end
    print_repeat(indent, io, ' ')
    print(io, right)
    return indent + length(right)
end

"""
The implementation of `pprint` and `pformat`.

You can extend pretty print support via such boilerplate:
```
    function PrettyPrint.pp_impl(io, data::YourType, indent::Int)
        print(io, data)
    end
```
"""
function pp_impl(io, seq::Vector, indent::Indentation)::Indentation
    pprint_for_seq(io, "[", "]", seq, indent)
end

function pp_impl(io, seq::Tuple, indent::Indentation)::Indentation
    pprint_for_seq(io, "(", ")", seq, indent)
end

function pp_impl(io, seq::AbstractSet, indent::Indentation)::Indentation
    pprint_for_seq(io, "{", "}", collect(seq), indent)
end

function pp_impl(io, seq::AbstractDict, indent::Indentation)::Indentation
    pprint_for_seq(io, "{", "}", [PPPair(k, " : ", v) for (k, v) in seq], indent)
end

function pp_impl(io, p::PPPair{K, V}, indent)::Indentation where {K, V}
    indent = pp_impl(io, p.name, indent)
    print(io, p.sep)
    indent += length(p.sep)
    pp_impl(io, p.content, indent)
end

function pp_impl(io, p::Pair, indent::Indentation)::Indentation
    pp_impl(io, PPPair(p.first, " => ", p.second), indent)
end

pp_impl(io, p, indent)::Indentation =
    if is_static_t(typeof(p))
        pp_impl_static(io, p, indent)
    else
        pp_impl_dynamic(io, p, indent)
    end

@generated function pp_impl_static(io, data::T, indent)::Indentation where T
    if is_atom_t(T)
        return quote s = repr(data); print(io, s); length(s) + indent  end
    end
    fields = fieldnames(T)
    exp = Expr(:tuple)
    for field in fields
        push!(exp.args, Expr(:call, PPPair, PPMetaSymbol(field), "=", :(data.$field)))
    end
    left = string(T) * "("
    right = ")"
    :(pprint_for_seq(io, $left, $right, $exp, indent))
end

function pp_impl_dynamic(io, data, indent)::Indentation
    if is_atom(data)
        s = repr(data)
        print(io, s)
        return length(s) + indent
    end

    fields = propertynames(data)
    left = string(typeof(data)) * "("
    right = ")"
    seq = [PPPair(PPMetaSymbol(field), "=", getproperty(data, field)) for field in fields]
    pprint_for_seq(io, left, right, seq, indent)
end

"""
pretty print to IO.
"""
function pprint(io, data, indent)
    pp_impl(io, data, indent)
    nothing
end

function pprint(io, data)
    pprint(io, data, 0)
end

function pprint(data)
    pprint(stdout, data)
end

"""
return pretty formatted string
"""
function pformat(data)
    io = IOBuffer()
    pprint(io, data)
    String(take!(io))
end

pp_impl(io, s :: String, indent) = begin
    s = repr(s)
    print(io, s)
    length(s) + indent
end

pprintln(io, s) = begin
    pprint(io, s)
    print(io, '\n')
end

pprintln(s) = pprintln(stdout, s)

@deprecate pprint(io, data, indent, newline) pprint(io, data, indent)
@deprecate pp_impl(io, data, indent, newline) pp_impl(io, data, indent)

end # module


text:
mutable struct Callback
    callback::Function
    name::String

    function Callback(callback,name)
        cb = new(callback,name)
        callbacks[name] = cb
        finalizer(cb, cb -> delete!(callbacks, cb.name))
        cb
    end
end
const callbacks = Dict{String,Callback}() # name => args -> f(args...)

# Empty callback to notify the Server that a new page is loaded and its WebSocket is ready.
Callback("connected") do client, route, id, data
    println("New connection established (ID: $(id)).")
    notify(conditions["connected"])
end

# Callback used to cleanup when the browser navigates away from the page.
Callback("unloaded") do client, route, id, data
    close(client)
    delete!(pages[route].sessions,id)
end

# Callback used to message a specified connected browser.
Callback("message") do client, route, id, data
    message(route,data["args"])
end

# Callback used to broadcast to all connected browsers.
Callback("broadcast") do client, route, id, data
    broadcast(route,data["args"])
end

# Callback used for blocking Julia control flow until notified by the WebSocket.
Callback("notify") do client, route, id, data
    name = data["args"]
    if haskey(conditions,name)
        notify(conditions[name])
    else
        warn("""Condition "$name" was not found.""")
    end
end


text:
# ROV visualization

"""
    $(SIGNATURES)

Brings model data into context, but generating or loading meshes, Geometries, etc as required by requested model,
and must still be drawn with the Visualizer object using mdl(vis).
"""
function loadmodel(model::Symbol=:rov;
    color=RGBA(0., 1.0, 0.5, 0.3),
    offset = Translation(0.,0,0) ∘ LinearMap(Rotations.Quat(1.0,0,0,0))  )
  #
  if model==:rov
    @show file = joinpath(dirname(@__FILE__), ".." , "data", "models", "rov3.obj")
    rov = load(file)
    rovdata = GeometryData(rov)
    rovdata.color = color
    offset = Translation(-0.5,0,0.25) ∘ LinearMap(Rotations.Quat(0,0,0,1.0))
    return DrawROV(rovdata, 99, :rov, offset)
  elseif model == :scene01
    boxdata = GeometryData(HyperRectangle(Vec(4.0,0,0), Vec(5.0,5.0,0)), RGBA(0.5,0.1,0.0,0.5))
    # boxdata.color =
    offset = Translation(0.0,0,0) ∘ LinearMap(CoordinateTransformations.AngleAxis(0.0,0,0,1.0))
    return DrawScene(boxdata, :scene01, offset)
  elseif model == :dock
    @show file = joinpath(dirname(@__FILE__), ".." , "data", "models", "dock.obj")
    rov = load(file)
    rovdata = GeometryData(rov)
    rovdata.color = color
    offset = Translation(0.0,0,0) ∘ LinearMap(Rotations.Quat(1.0,0,0,0.0))
    return DrawScene(rovdata, :dock, offset)
  else
    error("Don't recognize requested $(string(model)) model.")
  end
end

# t∈[0,1], about
function parameterizeArcAffineMap(t, as::ArcPointsRangeSolve; initrot::Rotation=Rotations.Quat(1.0,0,0,0))
  the = t*as.angle
  rot = Rotations.AngleAxis(the, as.axis...)
  dp = as.x1-as.center
  arc = LinearMap(rot) ∘ Translation(dp...)

  cent = Translation(as.center)
  return cent ∘ arc ∘ LinearMap(initrot)
end

# as = ArcPointsRangeSolve(
#       [0.0;0;0.0],
#       [0.0;6.0;6.0],
#       [0.0;12.0;0.0], 6.0)
# findaxiscenter!(as)
# # @show as.center, as.axis, as.angle, as.r
# animatearc(vc, rovt, as, initrot=Rotations.AngleAxis(pi/2,0,0,1.0))
#
#
# as = ArcPointsRangeSolve(
#       [2.5;-1;0.0],
#       [-4.5;4.0;0.0],
#       [9.5;6.0;0.0], 7.0)
# findaxiscenter!(as)
# # @show as.center, as.axis, as.angle, as.r
# animatearc(vc, rovt, as, initrot=Rotations.AngleAxis(pi/2,0,0,1.0))
#

#


text:
#=
Copyright (c) 2015, Intel Corporation

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 Intel Corporation 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 OWNER 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.
=#

include("utils.jl")

function pcg(x, A, b, M, tol, maxiter)
    r = b - A * x
    normr0 = norm(r)
    rel_err = 1.0
    z::Vector{Float64} = M \ r   # Somehow, types for "M\r" is not inferred.
                                 # The explicit typing makes it compile and shows
                                 # the call replacement, but has a running error
                                 # in type conversion. So this example is for
                                 # manual checking of call replacement only.
    p = copy(z) #NOTE: do not write "p=z"! That would make p and z aliased (the same variable)
    rz = dot(r, z)
    k = 1
    while k <= maxiter
        old_rz = rz
        Ap = A*p # Ap = Sparso.SpMV(A, p) # This takes most time. Compiler can reorder A to make faster
        alpha = old_rz / dot(p, Ap)
        x += alpha * p
        r -= alpha * Ap
        rel_err = norm(r)/normr0
        if rel_err < tol
            break
        end
        z = M \ r
        rz = dot(r, z)
        beta = rz/old_rz
        p = z + beta * p
        k += 1
    end
    return x, k, rel_err
end

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
"""
$(DocStringExtensions.README)
"""
module StochasticDiffEq

using DocStringExtensions
  import RandomNumbers: Xorshifts

  using Reexport
  @reexport using DiffEqBase

  import OrdinaryDiffEq

  using UnPack, RecursiveArrayTools, DataStructures
  using DiffEqNoiseProcess, Random, ArrayInterface
  using NLsolve, ForwardDiff, StaticArrays, MuladdMacro, FiniteDiff, Base.Threads

  import DiffEqBase: ODE_DEFAULT_NORM, ODE_DEFAULT_ISOUTOFDOMAIN,
         ODE_DEFAULT_PROG_MESSAGE, ODE_DEFAULT_UNSTABLE_CHECK

  using DiffEqBase: DiffEqArrayOperator

  using DiffEqBase: TimeGradientWrapper, UJacobianWrapper, TimeDerivativeWrapper, UDerivativeWrapper

  import RecursiveArrayTools: chain

  using Logging, SparseArrays

  using LinearAlgebra, Random, FillArrays

  import ForwardDiff.Dual

  import DiffEqBase: step!, initialize!, DEAlgorithm,
                     AbstractSDEAlgorithm, AbstractRODEAlgorithm, DEIntegrator, AbstractDiffEqInterpolation,
                     DECache, AbstractSDEIntegrator, AbstractRODEIntegrator, AbstractContinuousCallback,
                     Tableau

  # Integrator Interface
  import DiffEqBase: resize!,deleteat!,addat!,full_cache,user_cache,u_cache,du_cache,
                     rand_cache,ratenoise_cache,
                     resize_non_user_cache!,deleteat_non_user_cache!,addat_non_user_cache!,
                     terminate!,get_du, get_dt,get_proposed_dt,set_proposed_dt!,
                     u_modified!,savevalues!,add_tstop!,add_saveat!,set_reltol!,
                     set_abstol!, postamble!, last_step_failed, has_Wfact, has_jac

  using DiffEqBase: check_error!, is_diagonal_noise, @..

  using DiffEqBase: nlsolvefail, isnewton, set_new_W!, get_W, iipnlsolve, oopnlsolve, _vec, _reshape

  using DiffEqBase: NLSolver

  using DiffEqBase: FastConvergence, Convergence, SlowConvergence, VerySlowConvergence, Divergence

  import DiffEqBase: calculate_residuals, calculate_residuals!, nlsolve_f, unwrap_cache, islinear

  import DiffEqBase: iip_get_uf, oop_get_uf, build_jac_config

  using SparseDiffTools: forwarddiff_color_jacobian!, ForwardColorJacCache

  const CompiledFloats = Union{Float32,Float64}

  import DiffEqJump
  import DiffEqJump: JumpProblem

  import Base.Threads
  @static if VERSION < v"1.3"
    seed_multiplier() = Threads.threadid()
  else
    seed_multiplier() = 1
  end

  include("misc_utils.jl")
  include("algorithms.jl")
  include("options_type.jl")
  include("interp_func.jl")
  include("caches/cache_types.jl")
  include("caches/basic_method_caches.jl")
  include("caches/explicit_3s_mil_methods.jl")
  include("caches/lamba_caches.jl")
  include("caches/iif_caches.jl")
  include("caches/sdirk_caches.jl")
  include("caches/implicit_split_step_caches.jl")
  include("caches/sra_caches.jl")
  include("caches/rossler_caches.jl")
  include("caches/srk_weak_caches.jl")
  include("caches/kencarp_caches.jl")
  include("caches/predcorr_caches.jl")
  include("caches/SROCK_caches.jl")
  include("caches/tau_caches.jl")
  include("caches/dynamical_caches.jl")
  include("integrators/type.jl")
  include("dense.jl")
  include("alg_utils.jl")
  include("integrators/stepsize_controllers.jl")
  include("integrators/integrator_utils.jl")
  include("cache_utils.jl")
  include("integrators/integrator_interface.jl")
  include("iterator_interface.jl")
  include("composite_solution.jl")
  include("solve.jl")
  include("initdt.jl")
  include("perform_step/low_order.jl")
  include("perform_step/explicit_3s_mil_methods.jl")
  include("perform_step/lamba.jl")
  include("perform_step/iif.jl")
  include("perform_step/sri.jl")
  include("perform_step/sra.jl")
  include("perform_step/srk_weak.jl")
  include("perform_step/sdirk.jl")
  include("perform_step/implicit_split_step.jl")
  include("perform_step/kencarp.jl")
  include("perform_step/predcorr.jl")
  include("perform_step/split.jl")
  include("perform_step/composite.jl")
  include("perform_step/SROCK_perform_step.jl")
  include("perform_step/tau_leaping.jl")
  include("perform_step/dynamical.jl")
  include("tableaus.jl")
  include("SROCK_tableaus.jl")
  include("iterated_integrals.jl")
  include("SROCK_utils.jl")
  include("composite_algs.jl")
  include("weak_utils.jl")

  export StochasticDiffEqAlgorithm, StochasticDiffEqAdaptiveAlgorithm,
          StochasticCompositeAlgorithm

  export EM, LambaEM, PCEuler, RKMil, SRA, SRI, SRIW1,
         SRA1, SRA2, SRA3,
         SOSRA, SOSRA2, RKMilCommute, RKMilGeneral,
         SRIW2, SOSRI, SOSRI2, SKenCarp,
         SROCK1, SROCK2, SROCKEM, SKSROCK, TangXiaoSROCK2, KomBurSROCK2, SROCKC2,
         WangLi3SMil_A, WangLi3SMil_B, WangLi3SMil_C, WangLi3SMil_D, WangLi3SMil_E, WangLi3SMil_F,
         AutoSOSRI2, AutoSOSRA2,
         DRI1, DRI1NM, RI1, RI3, RI5, RI6, RDI1WM, RDI2WM, RDI3WM, RDI4WM,
         RS1, RS2,
         PL1WM, PL1WMA,
         NON, COM, NON2

  export SIEA, SMEA, SIEB, SMEB

  export EulerHeun, LambaEulerHeun

  export SimplifiedEM

  export SplitEM, IIF1M, IIF2M, IIF1Mil

  export ImplicitEM, ImplicitEulerHeun, ISSEM, ISSEulerHeun,
         ImplicitRKMil, STrapezoid, SImplicitMidpoint

  export TauLeaping, CaoTauLeaping

  export BAOAB

  export StochasticDiffEqRODEAlgorithm, StochasticDiffEqRODEAdaptiveAlgorithm,
         StochasticDiffEqRODECompositeAlgorithm

  export RandomEM

  export IteratedIntegralApprox, IICommutative, IIWiktorsson

  #General Functions
  export solve, init, solve!, step!

  #Misc Tools
  export checkSRIOrder, checkSRAOrder,  checkRIOrder, checkRSOrder,
         checkNONOrder,
         constructSRIW1, constructSRA1,
         constructDRI1, constructRI1, constructRI3, constructRI5, constructRI6,
         constructRDI1WM, constructRDI2WM, constructRDI3WM, constructRDI4WM,
         constructRS1, constructRS2,
         constructNON, constructNON2

end # module


text:
using Base.Broadcast: broadcastable, Broadcasted, materialize
using Swizzles.NamedArrays
using Swizzles.ArrayifiedArrays
using Swizzles

foo(x) = x

@testset "NamedArrays" begin
    A = rand(3,3)
    B = rand(3,3)
    C = rand(3,3)

    bc = Delay().(A .+ B .+ C)
    @test typeof(lift_names(bc)) <: Broadcasted{<:Any, <:Any, typeof(+), <:Tuple{<:Broadcasted{<:Any, <:Any, typeof(+), <:Tuple{<:NamedArray{<:Any, <:Any, <:Any, :obj1}, <:NamedArray{<:Any, <:Any, <:Any, :obj2}}}, NamedArray{<:Any, <:Any, <:Any, :obj3}}}
    @test name(lift_names(bc).args[1].args[1]) == :obj1
    @test name(lift_names(bc).args[1].args[2]) == :obj2
    @test name(lift_names(bc).args[2]) == :obj3
    @test materialize(lift_names(bc)) == materialize(bc)

    bc = Delay().(Swizzle(+).(A .+ B))
    @test typeof(lift_names(bc)) <: Swizzles.SwizzledArray{<:Any, <:Any, typeof(+), (), <:Any, <:ArrayifiedArray{<:Any, <:Any, <:Broadcasted{<:Any, <:Any, typeof(+), <:Tuple{<:NamedArray{<:Any, <:Any, <:Any, :obj1}, <:NamedArray{<:Any, <:Any, <:Any, :obj2}}}}}
end


text:
# ## Непосредственная реализация  проверки сопряженной задачи
#
using Test;
using LaTeXStrings;
using Plots;

include("test/adjoint_check.jl")

# Посмотрим на результат решения в сравнении с аналитическим решением

d = [missing, missing];
dd = [missing missing; missing missing];
make_gif(ψ, Xₙ, Tₘ, dd, dd, dd, d, d, ψ_model; name = "direct_check.gif")

# Найдем абсолютную погрешность численного решения

err = ψ .- ψ_model
heatmap(Xₙ, Tₘ, err', xlabel=L"X_n", ylabel=L"T_m", title="Absolute Error", size=(1200, 800))


text:
# testfile.jl

function passthrough(inp, outp)
    for oidx in 1:size(outp,2)
        # modulo output index to a usable input index
        iidx = mod(oidx, size(inp,2)) + 1

        # copy in-channel (iidx) to out-channel (oidx)
        outp[oidx,:] = inp[iidx,:]
    end
end


text:
module mq4gen1
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: Q4blockx
using Test
function test()
    connectivity = Q4blockx([0.0, 1.0, 3.0], [0.0, 1.0, 3.0])
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (9, 4)
    vtkwrite("mq4gen1a", connectivity)
    try rm("mq4gen1a.vtu"); catch end
    true
end
end
using .mq4gen1
mq4gen1.test()

module mq4gen1b
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: Q4blockx
using Test
function test()
    connectivity = Q4blockx([0.0, 1.0, 3.0], [0.0, 1.0, 3.0])
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (9, 4)
    vtkwrite("mq4gen1b", connectivity)
    try rm("mq4gen1b.vtu"); catch end
    true
end
end
using .mq4gen1b
mq4gen1b.test()

module mq4gen2
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: Q4block
using Test
function test()
    connectivity = Q4block(1.0, 2.0, 7, 9)
    # @show (nshapes(connectivity.right), nshapes(connectivity.left))
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (80, 63)
    vtkwrite("mq4gen2", connectivity)
    try rm("mq4gen2.vtu"); catch end
    true
end
end
using .mq4gen2
mq4gen2.test()

module mq4gen3
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: export_MESH
using MeshSteward: Q4blockx, gradedspace
using Test
function test()
    connectivity = Q4blockx([1.0, 2.0, 3.0], gradedspace(0.0, 5.0, 7, 2))
    # @show (nshapes(connectivity.right), nshapes(connectivity.left))
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (21, 12)
    vtkwrite("mq4gen3", connectivity)
    try rm("mq4gen3.vtu"); catch end
    export_MESH("mq4gen3", connectivity)
    try rm("mq4gen3.mesh"); catch end
    try rm("mq4gen3-xyz.dat"); catch end
    try rm("mq4gen3-conn.dat"); catch end
    true
end
end
using .mq4gen3
mq4gen3.test()

module mq4gen4
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: export_MESH
using MeshSteward: Q4quadrilateral
using Test
function test()
    connectivity = Q4quadrilateral([1.0 0.0; 1.5 1.7; -0.5 0.9; -0.1 -0.1], 3, 4)
    # @show (nshapes(connectivity.right), nshapes(connectivity.left))
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (20, 12)
    vtkwrite("mq4gen4", connectivity)
    try rm("mq4gen4.vtu"); catch end
    export_MESH("mq4gen4", connectivity)
    try rm("mq4gen4.mesh"); catch end
    try rm("mq4gen4-xyz.dat"); catch end
    try rm("mq4gen4-conn.dat"); catch end
    true
end
end
using .mq4gen4
mq4gen4.test()

module mq4gen5
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: export_MESH
using MeshSteward: Q4quadrilateral
using Test
function test()
    connectivity = Q4quadrilateral([1.0 0.0; 1.5 1.7], 6, 4)
    # @show (nshapes(connectivity.right), nshapes(connectivity.left))
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (35, 24)
    vtkwrite("mq4gen5", connectivity)
    try rm("mq4gen5.vtu"); catch end
    export_MESH("mq4gen5", connectivity)
    try rm("mq4gen5.mesh"); catch end
    try rm("mq4gen5-xyz.dat"); catch end
    try rm("mq4gen5-conn.dat"); catch end
    true
end
end
using .mq4gen5
mq4gen5.test()

module mq4gen6
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: export_MESH
using MeshSteward: Q4quadrilateral
using Test
function test()
    connectivity = Q4quadrilateral([1.0 0.0 0.0; 1.5 1.7 -0.1; 2.1 1.7 0.5; -1.0 1.0 0.3], 6, 4)
    # @show (nshapes(connectivity.right), nshapes(connectivity.left))
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (35, 24)
    vtkwrite("mq4gen6", connectivity)
    try rm("mq4gen6.vtu"); catch end
    export_MESH("mq4gen6", connectivity)
    try rm("mq4gen6.mesh"); catch end
    try rm("mq4gen6-xyz.dat"); catch end
    try rm("mq4gen6-conn.dat"); catch end
    true
end
end
using .mq4gen6
mq4gen6.test()

module mq4gen7
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: export_MESH
using MeshSteward: Q4quadrilateral, mergeirs
using Test
function test()
    N = 2
    c1 = Q4quadrilateral([-1 -1; -0.2 -1; -0.1 -0.2; -1 0.8], N, N)
    c2 = Q4quadrilateral([-0.2 -1; 1 -1; 1 -0.5; -0.1 -0.2], N, N)
    c3 = Q4quadrilateral([1 -0.5; 1 1; 0.3 1; -0.1 -0.2], N, N)
    c4 = Q4quadrilateral([0.3 1; -1 1; -1 0.8; -0.1 -0.2], N, N)
    c  = mergeirs(c1, c2, 0.001)
    c  = mergeirs(c, c3, 0.001)
    c  = mergeirs(c, c4, 0.001)
    @test nshapes(c.left) == 4 * N ^ 2

    vtkwrite("mq4gen7", c)
    try rm("mq4gen7.vtu"); catch end
    # export_MESH("mq4gen7", connectivity)
    # try rm("mq4gen7.mesh"); catch end
    # try rm("mq4gen7-xyz.dat"); catch end
    # try rm("mq4gen7-conn.dat"); catch end
    true
end
end
using .mq4gen7
mq4gen7.test()

module mq4gen8
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: export_MESH
using MeshSteward: Q4blockwdistortion
using Test
function test()
    c1 = Q4blockwdistortion(3.0, 1.0, 15, 6)

    vtkwrite("mq4gen8", c1)
    try rm("mq4gen8.vtu"); catch end
    # export_MESH("mq4gen8", connectivity)
    # try rm("mq4gen8.mesh"); catch end
    # try rm("mq4gen8-xyz.dat"); catch end
    # try rm("mq4gen8-conn.dat"); catch end
    true
end
end
using .mq4gen8
mq4gen8.test()


text:
#An N-dimensional PSF is modeled as a product of N gaussians oriented orthogonally
#The magnitude of the PSF at a given location in space can be found by
#indexing the Psf.
#The gaussians are centered at the origin; indexing at [0,0...] returns the maximum value
#See also: PsfFit
struct Psf{T,N} <: AbstractArray{T,N}
    sigmas::NTuple{N,T} #beam halfwidth, i.e. sigma in y = e^-(2x^2/sigma^2)
end
Psf(args...) = Psf((args...,))

function Psf(;NA=0.5, M=10.0, lambda_nm=500, um_per_pix=6.5, um_per_slice=5.0, n=1.333)
    lateral_optical_um = approx_sigma_lateral(lambda_nm/1000, NA)
    axial_optical_um = approx_sigma_axial(lambda_nm/1000, NA, n)
    lateral_pixel = lateral_optical_um/(um_per_pix/M)
    axial_pixel = axial_optical_um/um_per_slice
    return Psf((lateral_pixel, lateral_pixel, axial_pixel))
end

size_optics(psf::Psf{T,N}) where {T,N} = psf.sigmas
size(psf::Psf{T,N}) where {T,N} = size_optics(psf)
size_stats(psf::Psf{T,N}) where {T,N} = size_optics(psf) ./ 2

axes(psf::Psf{T,N}) where {T,N} = ntuple(i->nothing, N)
Base.show(io::IO, psf::Psf{T,N}) where {T,N} = print(io, "$N-dimensional Psf with size $(size(psf)) (halfwidths of Gaussian approximation)\n")
Base.show(io::IO, ::MIME"text/plain", psf::Psf) = show(io, psf)

getindex(psf::Psf{T,N}, args...) where {T,N} = getindex(psf, T.(args)...)
getindex(psf::Psf{T,N}, I::Vararg{T,N}) where {T,N} = prod(eval_gaussian.((I...,), size(psf) ./ 2))
setindex!(psf::Psf, args...) = error("setindex! is not appropriate for a Psf")

fwhm(psf::Psf, dim::Int) = fwhm_optics(size(psf)[dim])
fwhm_stats(psf::Psf, dim::Int) = fwhm_stats(size(psf)[dim])

lightsheet_psf(det_psf::Psf{T,3}, NA_illum, lambda_illum_nm, um_per_slice =1) where {T} =
    lightsheet_psf(det_psf, approx_sigma_lateral(lambda_illum_nm/1000, NA_illum)/um_per_slice)

function lightsheet_psf(det_psf::Psf{T,3}, sheet_sigma) where {T}
    dax = size(det_psf)[3]
    combinedax = lightsheet_axial_combined_psf(dax, sheet_sigma)
    return Psf((size(det_psf)[1:2]..., combinedax))
end

#PsfFit contains the result of fitting a point spread function to an image img
#See also: Psf
struct PsfFit{Tfit,N,Timg,TI}
    parent_inds::TI #image array to which we fit the psf
    bias::Timg #black level of image (100 for PCO cameras)
    psf::Psf{Tfit,N}
    shift::NTuple{N,Tfit} #amount by which fitted psf is shifted relative to the center of the img array
    scale::NTuple{N,Tfit} #amount by which fitted psf is scaled (scale of 1 means peak height is 1)
    quality::NTuple{N,Tfit} #fractions of variance explained by the gaussian fits in each dimension
end

psf(fit::PsfFit) = fit.psf
quality(fit::PsfFit) = fit.quality
size(fit::PsfFit) = size(psf(fit))
shift(fit::PsfFit) = fit.shift
scale(fit::PsfFit) = fit.scale
Base.parentindices(fit::PsfFit) = fit.parent_inds

function sum_otherdims(img, dim::Int, precision)
    sz = ntuple(d->d!=dim ? 1 : size(img,dim), Val(ndims(img)))
    reduced = Array{precision}(undef, sz...)
    sum_otherdims!(reduced, img)
end

function sum_otherdims!(reduced, img)
    fill!(reduced, 0)
    rmax = last(CartesianIndices(reduced))
    for I in CartesianIndices(img)
        reduced[min(rmax, I)] += img[I]
    end
    return reduced
end

#Fits a gaussian along each dimension independently
function PsfFit(img::OffsetArray{Timg,N}, bias;
                    guess_sigmas=ones(N),
                    guess_scales=ones(N),
                    precision=Float64) where {Timg,N}
    stds = zeros(precision,N)
    shifts, scales, quality = copy(stds), copy(stds), copy(stds)
    pimg = parent(img)
    for i = 1:N
        reduced = sum_otherdims(pimg, i, precision)
        newbias = precision(bias)*(div(length(img), size(img,i))) #sum the pixelwise bias in this slice
        curdata = reshape(reduced, length(reduced))
        curdata = OffsetArray(curdata, axes(img)[i])
        stds[i], shifts[i], scales[i], fval = fit_gaussian_1d(curdata, newbias;
                                                             guess_std=guess_sigmas[i] / 2,
                                                             guess_scale=guess_scales[i])
        quality[i] = 1.0-fval/sum(curdata.^2)
    end
    sigmas = stds .* 2 #convert from std to optics halfwidth convention
    psf = Psf((sigmas...,))
    fit = PsfFit(parentindices(parent(img)), eltype(img)(bias), psf, (shifts...,), (scales...,), (quality...,))
    return fit
end

max_stars_repo_name: HolyLab/ExtractPSF

max_forks_repo_name: HolyLab/ExtractPSF


text:
function readindent(io)
    indent = 0
    while !eof(io) && (p = Base.peek(io); p == 0x20 || p == 0x09)
        p == 0x20 && (indent += 1)
        p == 0x09 && (indent += 4)
        read(io, UInt8)
    end
    return indent
end

function get_lines(s::String)
    io = IOBuffer(s)
    lines = Tuple{Int,Int}[(0, readindent(io))] # startbyte, indent
    while !eof(io)
        c = read(io, Char)
        if c == '\r' && Base.peek(io) == 0x0a
            c = read(io, Char)
            push!(lines, (position(io), readindent(io)))
        elseif c == '\n'
            push!(lines, (position(io), readindent(io)))
        end
    end
    first(last(lines)) != sizeof(s) && push!(lines, (sizeof(s), 0))
    lines
end

function line_of(offset, lines)
    offset > first(last(lines)) && error()
    offset == 0 && return 1
    for i = 1:length(lines) - 1
        if lines[i][1] < offset <= lines[i + 1][1]
            return i
        end
    end
    return length(lines)
end

function issameline(offset1, offset2, lines)
    line_of(offset1, lines) == line_of(offset2, lines)
end

function get_expr(x, offset, pos = 0)
    if pos > offset
        return nothing
    end
    if x.args !== nothing
        for a in x.args
            if pos < offset <= (pos + a.fullspan)
                return get_expr(a, offset, pos)
            end
            pos += a.fullspan
        end
    elseif pos == 0
        return x
    elseif (pos < offset <= (pos + x.fullspan))
        return x
    end
end


text:
################################################################################
#
#   bgjt.jl : functions for the BGJT algorithm
#
################################################################################

export dlogBGJT, smsrField

# Polynomial manipulation

"""
    homogene{T <: PolyElem}(P::T, h0::T, h1::T)

Homogenize ``π`` with `h0` and `h1` seen as indeterminates. The polynomial
``π`` is the polynomial `P` with coefficients powered by ``q``, the
characteristic in the context.
"""
function homogene{T <: PolyElem}(P::T, h0::T, h1::T)

    # We compute the characteristic of the context and the degree of `P`
    R = parent(P)
    q::Int = sqrt(length(base_ring(R)))
    H = R()
    d = degree(P)

    # Then we compute the homogenized polymonial and return it
    for i in 0:d
        H += (coeff(P, i)^q)*(h0^i)*(h1^(d-i))
    end
    return H
end

"""
    makeEquation{T <: RingElem, Y <: PolyElem}(m::Nemo.GenMat{T},
                                               P::Y, h0::Y, h1::Y)

Compute the denomimator of the left-side of the equation generated by `m`.

See reference **[1]** in the documentation of the module for more informations.
"""
function makeEquation{T <: RingElem, Y <: PolyElem}(m::Nemo.GenMat{T},
                                                    P::Y, h0::Y, h1::Y)

    # We compute some constants depending on the context
    a, b, c, d = m[1, 1], m[1, 2], m[2, 1], m[2, 2]
    D = degree(P)
    q::Int = sqrt(length(base_ring(P)))

    # We compute the homogenized version of `P`
    H = homogene(P, h0, h1)

    # And we compute the denominator of the left side of the equation defined
    # in the paper of Barbulescu et al.
    return ((a^q)*H + (b^q)*h1^D)*(c*P+d) - (a*P+b)*((c^q)*H + (d^q)*h1^D)
end

"""
    fillMatrixBGJT!(M::MatElem, j::Integer, m::MatElem, F::Nemo.Field)

Fill the `j`-th column of the matrix `M`, using the action of `m` on ``P_1(F)``
and returns the unit generated in the process.

See reference **[1]** for more informations, the references are listed in the
documentation of the module.
"""
function fillMatrixBGJT!(M::MatElem, j::Integer, m::MatElem, F::Nemo.Field)

    # We compute `q`
    q::Int = sqrt(length(F))

    # We set some constants to the coefficients of `m`
    a, b, c, d = m[1, 1], m[1, 2], m[2, 1], m[2, 2]

    # We compute the one of our ring
    Z = parent(M[1, 1])
    u = Z(1)

    # We first set our unit to 1 and our index to 0
    i = 0
    unit = F(1)

    # We iterate over the projective line P_1(F), with elements of the form
    # (y, 1) where y is in F
    for y in F
        # In fact the columns of `M` can be indexed by elements of P_1(F)
        i += 1

        # We compute m⋅(y, 1) = (α, β)
        α = a*y + b
        β = c*y + d
        # If (α, β) is in P_1(GF(q)), meaning if δ = α/β is in GF(q) or (α, β) is
        # the infinite element, we set M[i, j] to 1
        if β != 0
            δ = divexact(α, β)
            if δ^q == δ
                M[i, j] = u

                # We compute the constant (λ in [2]) needed for the two
                # sides of the equation to match
                tmp = -c*δ+a
                if tmp != 0
                    unit *= tmp
                else
                    unit *= -d*δ + b
                end
            end

        # Else we are in the case where m⋅(y, 1) = (., 0) = infinite
        # element = (-1, 0), so we are in P_1(GF(q))
        else
            M[i, j] = u
            if c != 0
                unit *= c
            else
                unit *= d
            end
        end
    end

    # We do the same for the last element of P_1(F)
    # the infinite element (-1, 0)
    i += 1
    if c != 0
        δ = divexact(a, c)
        if δ^q == δ
            M[i, j] = u
        # In this case tmp = -c*δ + a = -c*a/c + a = 0, so no need to test
        unit *= -d*δ + b
        end
    else
        M[i, j] = u
        unit *= d
    end
    return unit
end

"""
    descentBGJT{T <: PolyElem}(L::FactorsList, i0::Integer, h0::T, h1::T,
                                    card::Nemo.fmpz)

The descent phase of the BGJT algorithm.
"""
function descentBGJT{T <: PolyElem}(L::FactorsList, i0::Integer, h0::T, h1::T,
                                    card::Nemo.fmpz)

    # We set some constants, arrays, matrices
    F = base_ring(h0)
    elem, coef = L[i0]
    deg = degree(elem)
    smoothBound = ceil(Integer, deg/2)
    numerators = Array{fq_nmod_poly, 1}()
    charac::Int = characteristic(F)
    units = Array{fq_nmod, 1}()
    x = gen(F)
    j = 1

    S = MatrixSpace(ZZ, charac^2+1,charac^3+charac+1)
    M = zero(S)
    Pq = pglCosets(x)

    # We iterate over Pq = PGL(P_1(F_q²))/PGL(P_1(F_q)) to create new equations
    # involving P and its translations P - μ with μ in F_q², we keep only the
    # one with a smooth left side and we fill a matrix to remember which
    # translations were used
    for m in Pq
        N = makeEquation(m, elem, h0, h1)

        if isSmooth(N, smoothBound)
            unit = fillMatrixBGJT!(M, j, m, F)
            push!(units, unit)
            j += 1
            push!(numerators, N)
        end
    end

    # We set the last column to the vector (1, 0, ..., 0), which
    # represent the polynomial P
    M[1,j] = 1
    M = subMatrix(M, charac^2 + 1, j)

    # We compute the row echelon form of M, such that M/det is reduced
    M, det = rref(M)

    # We compute the inverse of `det` mod `card`
    det %= card
    if det <= 0
        det += card
    end
    g, s = gcdinv(det, card)

    # If it is not invertible, we return the problematic factor
    if g != 1
        println("The following number was a factor")
        return g
    end

    # We compute a solution
    sol = Nemo.fmpz[(s*M[i,j])%card for i in 1:(charac^2+1)]
#    sol = Nemo.fmpz[M[i,j] for i in 1:(charac^2+1)]

    # We compute the coordinates of the pivots (because we have redundant
    # equations)
    piv = pivots(M, charac^2+1)

    # We add the new polynomials and their coefficients in our list
    for j in 1:(charac^2+1)
        fact = factor(numerators[piv[j]])
        for f in fact
            push!(L, f[1], f[2]*sol[j]*coef)
        end
        push!(L, h1, -deg*sol[j]*coef)
        leadcoef = coeff(numerators[piv[j]], degree(numerators[piv[j]]))
        L.unit *= (inv(units[piv[j]])*leadcoef)^(sol[j]*coef)
    end
    deleteat!(L, i0)
end

function descentBGJT{T <: PolyElem}(L::FactorsList, i0::Integer, F::Nemo.Field,
                                    h0::T, h1::T, card::Integer)
    return descentBGJT(L, i0, F, h0, h1, Nemo.fmpz(card))
end

"""
    dlogBGJT(P::Nemo.RingElem, K::bgjtContext, dlogs::Dict{Any, Any})

Compute the discrete logarithm of `P` in the field `K`.

The field `K` must have a sparse medium subfield representation. All the
informations needed in the algorithm, like the basis, are in the object `K`.
This algorithm need the discrete logarithms of the linear factors, they can be
computed using **linearDlog**.
"""
function dlogBGJT(P::Nemo.fq_nmod_poly, K::BgjtContext, dlogs::Dict{Any, Any})

    defPol = K.definingPolynomial
    R = parent(K.gen)
    # We first compute the discrete logarithm of `P` modulo
    # the small factors with Pohlig Hellman algorithm
    cardinality = length(base_ring(R))^degree(defPol)
    s, n = pohligHellman(cardinality, K.gen, P, defPol)

    # Then we work with the large factors, using the descent algorithm to
    # express log P in function of logarithms of polynomials with smaller degree
    # and we stop once we have only linear polynomials, which logarithms are
    # known
    N = div(Nemo.fmpz(cardinality-1), n)
    i = 1
    L = factorsList(P)

    while i < length(L.factors) + 1
        while degree(L.factors[i]) > 1
            descentBGJT(L, i, K.h0, K.h1, N)
        end
        i += 1
    end

    # Once we have only linear polynomials, we compute the linear sum
    # log P = Σ e_i × log P_i (mod N) where the e_is where calculated during the
    # descent, and where the log P_is are given
    S = sum([L.coefs[j]*dlogs[L.factors[j]] for j in 1:length(L.factors)])
    #S += dlogSmallField(K.characteristic, K.extensionDegree, K.gen,
    #                    R(L.unit), defPol)

    # Then we reconstruct the result using Chinese remindering theorem
    return crt(s, n, S, N)
end


text:
MLJBase.is_wrapper(::Type{<:EitherIteratedModel}) = true
MLJBase.caches_data_by_default(::Type{<:EitherIteratedModel}) = false
MLJBase.supports_weights(::Type{<:EitherIteratedModel{M}}) where M =
    MLJBase.supports_weights(M)
MLJBase.supports_class_weights(::Type{<:EitherIteratedModel{M}}) where M =
    MLJBase.supports_class_weights(M)
MLJBase.load_path(::Type{<:DeterministicIteratedModel}) =
    "MLJIteration.DeterministicIteratedModel"
MLJBase.load_path(::Type{<:ProbabilisticIteratedModel}) =
    "MLJIteration.ProbabilisticIteratedModel"
MLJBase.package_name(::Type{<:EitherIteratedModel}) = "MLJIteration"
MLJBase.package_uuid(::Type{<:EitherIteratedModel}) =
    "614be32b-d00c-4edb-bd02-1eb411ab5e55"
MLJBase.package_url(::Type{<:EitherIteratedModel}) =
    "https://github.com/JuliaAI/MLJIteration.jl"
MLJBase.package_license(::Type{<:EitherIteratedModel}) = "MIT"
MLJBase.is_pure_julia(::Type{<:EitherIteratedModel{M}}) where {T,M} =
    MLJBase.is_pure_julia(M)
MLJBase.input_scitype(::Type{<:EitherIteratedModel{M}}) where {T,M} =
    MLJBase.input_scitype(M)
MLJBase.output_scitype(::Type{<:EitherIteratedModel{M}}) where {T,M} =
    MLJBase.output_scitype(M)
MLJBase.target_scitype(::Type{<:EitherIteratedModel{M}}) where {T,M} =
    MLJBase.target_scitype(M)

max_stars_repo_name: ablaom/MLJIteration

max_issues_repo_name: ablaom/MLJIteration

max_forks_repo_name: ablaom/MLJIteration


text:
# to_inplace.jl - transform ExGraph to buffered/in-place operations

const Num = Number
const Vec = AbstractVector
const Mat = AbstractMatrix

const INPLACE_PHS = Set([:u, :v, :w, :x, :y, :z, :X, :Y, :Z])

const INPLACE_RULES = [
    Type[Mat, Mat] => :(Z = X' * Y) => :(mul!(Z, transpose(X), Y)),
    Type[Mat, Mat] => :(Z = X * Y') => :(mul!(Z, X, transpose(Y))),
    Type[Mat, Mat] => :(Z = X * Y) => :(mul!(Z, X, Y)),
]

function to_inplace(g::ExGraph)
    evaluate!(g)
    g = expand_fixed_sequences(g)
    g = fuse_broadcasting(g)
    res = :(begin end)
    for nd in g.tape
        if getcategory(nd) != :input
            vex = to_inplace(g, nd)
            # push!(res.args, simplify(subs_size(vex, sizes)))
            push!(res.args, vex)
        end
    end
    res = subs_bcast_with_dot(sanitize(res))
    return res
end

function to_inplace(g::ExGraph, nd::Union{ExNode{:call}, ExNode{:bcast}, ExNode{:opaque}})
    ex = expand_const(g, to_expr_kw(nd)) |> simplify
    if isa(nd, ExNode{:call}) && !iscall(ex.args[2])
        return to_inplace(g, convert_call(g, nd))
    end
    # try patterns
    dep_vals = [getvalue(g[dep]) for dep in dependencies(nd)]
    for (types, (pat, rpat)) in INPLACE_RULES
        if all(isa(val, T) for (val, T) in zip(dep_vals, types))
            pat = sanitize(pat)
            rex = tryrewrite(ex, pat, rpat; phs=INPLACE_PHS)
            if rex != nothing
                return rex
            end
        end
    end
    # try broadcasting
    if is_bcast_vec(nd)
        lhs_is_scalar = isa(getvalue(nd), Number)
        return make_elementwise(to_expr_kw(nd); lhs_is_scalar=lhs_is_scalar)
    end
    # if LHS is an array, use .= instead of =
    if isa(getvalue(nd), AbstractArray)
        rex = to_expr_kw(nd)
        rex.head = :.=
        return rex
    end
    return to_expr_kw(nd)
end

function to_inplace(g::ExGraph, nd::ExNode{:ctor})
    ex = to_expr_kw(nd)
    insert!(ex.args[2].args, 3, :mem)
    insert!(ex.args[2].args, 4, QuoteNode(varname(nd)))
    return ex
end

function to_inplace(g::ExGraph, nd::ExNode)
    return to_expr_kw(nd)
end

to_buffered = to_inplace

## @inplacerule

# TODO: copied from XGrad, fix
isparameters(a) = isa(a, Expr) && a.head == :parameters

function without_types(pat)
    rpat = copy(pat)
    for i=2:length(rpat.args)
        a = rpat.args[i]
        if !isparameters(a)  # parameters aren't modified
            rpat.args[i] = isa(a, Expr) ? a.args[1] : a
        end
    end
    return rpat
end

function get_arg_names(pat)
    return [isa(a, Expr) ? a.args[1] : a for a in pat.args[2:end] if !isparameters(a)]
end

function get_arg_types(pat)
    return [isa(a, Expr) ? eval(Base, a.args[2]) : Any for a in pat.args[2:end] if !isparameters(a)]
end

function add_inplace_rule(mod::Module, tpat, rpat)
    @assert tpat.head == :(=)
    @assert rpat.head == :call
    op = canonical(mod, tpat.args[2].args[1])
    rop = canonical(mod, rpat.args[1])
    types = get_arg_types(tpat.args[2])
    pat = :($(tpat.args[1]) = $(without_types(tpat.args[2])))
    # replace function names with canonical names
    pat = subs(pat, Dict(tpat.args[2].args[1] => op))
    rpat = subs(rpat, Dict(rpat.args[1] => rop))
    rule = types => pat => rpat
    push!(INPLACE_RULES, rule)
end

macro inplacerule(tpat, rpat)
    mod = @__MODULE__
    add_inplace_rule(mod, tpat, rpat)
    nothing
end


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

using Markdown
using InteractiveUtils

# ╔═╡ ce069b56-6212-11eb-359b-d567223d2ca7
using UnitfulAstro, Unitful

# ╔═╡ dabacae4-6211-11eb-24e6-0f01be7534c2
md"Assignment3: Stellar Collapse"

# ╔═╡ 13fff7ca-6212-11eb-077c-7947cfdc39c8
md"## Problem 1."

# ╔═╡ 2d4e997a-6212-11eb-3671-c55e57a35aa6
md"The energy released in a core collapse supernovae could be estimated by the difference in the gravitational energy pre and post the collapse."

# ╔═╡ 79329336-6213-11eb-08b6-21052f99e3d8
rini = 1e4 * u"km"

# ╔═╡ ec734190-6213-11eb-245d-49b11cbb0cdc
rfinal = 10 * u"km"

# ╔═╡ f4789400-6213-11eb-23f9-690d151ca69b
mass = 1.4 * u"Msun"

# ╔═╡ fe84e7fa-6213-11eb-3024-b50806285bd5
G = u"GMsun"/u"Msun"

# ╔═╡ 4cac0352-6212-11eb-2dd2-3370a67a7b01
function EngCoreCollapse(rini, rfinal, mass)
	G * mass^2 * ((1/rfinal) - (1/rini))
end

# ╔═╡ 2b5b3ef0-6214-11eb-1f2a-397bc6558586
Eng = uconvert(u"erg", EngCoreCollapse(rini, rfinal, mass))

# ╔═╡ 2ff74464-6216-11eb-3109-8b1189247127
md"## Problem 2"

# ╔═╡ 39137b50-6216-11eb-038d-932576307c40
md"We calculate the energy lost in the ejecta and photons and neutrons"

# ╔═╡ 50eb550e-6216-11eb-3b07-9f78ee7ea244
Mprogenitor = 10 * u"Msun"

# ╔═╡ cfb61ed2-6216-11eb-1422-b7ad3f259e40
Mejecta = Mprogenitor - mass

# ╔═╡ e21a9b20-6216-11eb-1115-c9e34718d03a
Vejecta = 10^4 * u"km/s"

# ╔═╡ f3ee18c2-6216-11eb-1c6b-6fdd191655e3
KE = uconvert(u"erg", (1/2)* Mejecta * Vejecta^2)

# ╔═╡ 1f000c50-6217-11eb-3b29-35e08d006c79
fractionKE = KE/Eng

# ╔═╡ 6261a06c-6217-11eb-0ed4-591915acb820
Luminosity = 2 * 10^8 * u"Lsun"

# ╔═╡ ab933340-6217-11eb-2388-9db82ff944a7
Tvisible = 60 * 24 * 3600 * u"s"

# ╔═╡ 90837fa6-6217-11eb-1801-839ebe2789a6
EMEnergy = uconvert(u"erg", Luminosity * u"s")

# ╔═╡ 3b24e44a-6218-11eb-396c-2932821c89e7
md"## Problem 3"

# ╔═╡ 42f1582a-6218-11eb-2999-bd0511099c52
Etotnutrino = Eng - KE - EMEnergy

# ╔═╡ 6644c924-6218-11eb-00ea-bda8f6c78709
Enutrino = uconvert(u"erg", 5 * 1e6 * u"eV")

# ╔═╡ 8899849c-6218-11eb-2773-d90d69839c9b
Nnutrino = Etotnutrino / Enutrino

# ╔═╡ 2985585e-6219-11eb-0a1d-13516cc0943a
dL = uconvert(u"m", 1e4 * u"pc")

# ╔═╡ 184cf092-6219-11eb-2d3a-75ef0c0f5580
Flux = Nnutrino / (4*pi*dL^2)

# ╔═╡ Cell order:
# ╟─dabacae4-6211-11eb-24e6-0f01be7534c2
# ╟─13fff7ca-6212-11eb-077c-7947cfdc39c8
# ╟─2d4e997a-6212-11eb-3671-c55e57a35aa6
# ╠═4cac0352-6212-11eb-2dd2-3370a67a7b01
# ╠═ce069b56-6212-11eb-359b-d567223d2ca7
# ╠═79329336-6213-11eb-08b6-21052f99e3d8
# ╠═ec734190-6213-11eb-245d-49b11cbb0cdc
# ╠═f4789400-6213-11eb-23f9-690d151ca69b
# ╠═fe84e7fa-6213-11eb-3024-b50806285bd5
# ╠═2b5b3ef0-6214-11eb-1f2a-397bc6558586
# ╟─2ff74464-6216-11eb-3109-8b1189247127
# ╟─39137b50-6216-11eb-038d-932576307c40
# ╠═50eb550e-6216-11eb-3b07-9f78ee7ea244
# ╠═cfb61ed2-6216-11eb-1422-b7ad3f259e40
# ╠═e21a9b20-6216-11eb-1115-c9e34718d03a
# ╠═f3ee18c2-6216-11eb-1c6b-6fdd191655e3
# ╠═1f000c50-6217-11eb-3b29-35e08d006c79
# ╠═6261a06c-6217-11eb-0ed4-591915acb820
# ╠═ab933340-6217-11eb-2388-9db82ff944a7
# ╠═90837fa6-6217-11eb-1801-839ebe2789a6
# ╟─3b24e44a-6218-11eb-396c-2932821c89e7
# ╠═42f1582a-6218-11eb-2999-bd0511099c52
# ╠═6644c924-6218-11eb-00ea-bda8f6c78709
# ╠═8899849c-6218-11eb-2773-d90d69839c9b
# ╠═2985585e-6219-11eb-0a1d-13516cc0943a
# ╠═184cf092-6219-11eb-2d3a-75ef0c0f5580


text:
using FermiCG
using Printf
using Test

@testset "TDMs" begin
    atoms = []
    push!(atoms,Atom(1,"H",[0,0,0]))
    push!(atoms,Atom(2,"H",[0,0,1]))
    push!(atoms,Atom(3,"H",[0,0,2]))
    push!(atoms,Atom(4,"H",[0,0,3]))
    push!(atoms,Atom(5,"H",[0,0,4]))
    push!(atoms,Atom(6,"H",[0,0,5]))
    basis = "6-31g"
    basis = "sto-3g"
    mol     = Molecule(0,1,atoms,basis)

    mf = FermiCG.pyscf_do_scf(mol)
    nbas = size(mf.mo_coeff)[1]
    ints = FermiCG.pyscf_build_ints(mol,mf.mo_coeff, zeros(nbas,nbas));
    e_fci, d1_fci, d2_fci = FermiCG.pyscf_fci(ints,1,1)
    @printf(" FCI Energy: %12.8f\n", e_fci)

    C = mf.mo_coeff
    Cl = FermiCG.localize(mf.mo_coeff,"lowdin",mf)
    FermiCG.pyscf_write_molden(mol,Cl,filename="lowdin.molden")
    S = FermiCG.get_ovlp(mf)
    U =  C' * S * Cl
    println(" Build Integrals")
    flush(stdout)
    ints = FermiCG.orbital_rotation(ints,U)
    println(" done.")
    flush(stdout)

    clusters    = [(1:4),(5:8),(9:12)]
    init_fspace = [(1,1),(1,1),(1,1)]

    clusters    = [(1:2),(3:4),(5:6)]
    init_fspace = [(1,1),(1,1),(1,1)]

    clusters    = [(1:4),(5:6)]
    init_fspace = [(2,1),(1,1)]

    max_roots = 20

    clusters = [Cluster(i,collect(clusters[i])) for i = 1:length(clusters)]
    display(clusters)

    #cluster_bases = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=0,
    #                                                   max_roots=2, init_fspace=init_fspace, delta_elec=2)
    cluster_bases = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=0, max_roots=100)

    #for cb in cluster_bases
    #    display(cb)
    #end

    cluster_ops = Vector{FermiCG.ClusterOps}()
    for ci in clusters
        push!(cluster_ops, FermiCG.ClusterOps(ci))
    end

    for ci in clusters
        cb = cluster_bases[ci.idx]

        cluster_ops[ci.idx]["A"], cluster_ops[ci.idx]["a"] = FermiCG.tdm_A(cb,"alpha")
        cluster_ops[ci.idx]["B"], cluster_ops[ci.idx]["b"] = FermiCG.tdm_A(cb,"beta")
        cluster_ops[ci.idx]["AA"], cluster_ops[ci.idx]["aa"] = FermiCG.tdm_AA(cb,"alpha")
        cluster_ops[ci.idx]["BB"], cluster_ops[ci.idx]["bb"] = FermiCG.tdm_AA(cb,"beta")
        cluster_ops[ci.idx]["Aa"] = FermiCG.tdm_Aa(cb,"alpha")
        cluster_ops[ci.idx]["Bb"] = FermiCG.tdm_Aa(cb,"beta")
        cluster_ops[ci.idx]["Ab"], cluster_ops[ci.idx]["Ba"] = FermiCG.tdm_Ab(cb)
        # remove BA and ba account for these terms
        @time cluster_ops[ci.idx]["AB"], cluster_ops[ci.idx]["ba"], cluster_ops[ci.idx]["BA"], cluster_ops[ci.idx]["ab"] = FermiCG.tdm_AB(cb)
        @time cluster_ops[ci.idx]["AAa"], cluster_ops[ci.idx]["Aaa"] = FermiCG.tdm_AAa(cb,"alpha")
        @time cluster_ops[ci.idx]["BBb"], cluster_ops[ci.idx]["Bbb"] = FermiCG.tdm_AAa(cb,"beta")
        #@time cluster_ops[ci.idx]["ABa"], cluster_ops[ci.idx]["Aba"] = FermiCG.tdm_ABa(cb,"alpha")
        #@time cluster_ops[ci.idx]["ABb"], cluster_ops[ci.idx]["Bba"] = FermiCG.tdm_ABa(cb,"beta")

    end

    ref1 = [     0.760602    0.194328
                -0.356922   -0.610415
                 0.540039   -0.637967
                -0.0494895   0.42735]
    test =cluster_ops[2]["A"][((1,1),(0,1))][1,:,:]
    @test isapprox(abs.(ref1), abs.(test), atol=1e-5)

    #for (ftrans,op) in cluster_ops[1]["B"]
    #    println(ftrans, "  :  ", size(op))
    #    display(op[1,:,:])
    #end

    ref2 = [ 0.692476    0.168256   -0.0394252    0.00396003
             0.633938    0.0369774   0.0770702   -0.0153966
            -0.0742221   0.646856   -0.636296     0.103912
            -0.308786   -0.0175901   0.00515862  -0.0633258
            -0.0776733   0.444367   -0.0201446   -0.621285
             0.0130287   0.0582555   0.348725    -0.653307]

    #display(ref2)
    #display(cluster_ops[1]["A"][((2, 4), (1, 4))][1,:,:])
    test =cluster_ops[1]["B"][((4, 2), (4, 1))][1,:,:]
    @test isapprox(abs.(ref2), abs.(test), atol=1e-5)

    ref = [  0.997916  0.064526
           -0.064526  0.997916]
    test =cluster_ops[2]["AA"][((2, 1), (0, 1))][1,2,:,:]
    @test isapprox(abs.(ref), abs.(test), atol=1e-5)

    ref = [ -0.210308     0.134555     -0.0105751   -0.0130626
           -0.146011     0.459975      0.0957622   -0.0409076
           -0.121296     0.404645      0.201137    -0.0233472
           -0.398129     0.129795     -0.0889375   -0.0373247
           0.160396    -0.276295      0.317716     0.172581
           -0.116331    -0.373659     -0.164785     0.0129931
           0.0401054   -0.17791      -0.248367    -0.0897069
           -0.12651      0.301717     -0.0496348   -0.248452
           -0.149806     0.0853372     0.433886     0.1877
           0.481946     0.0435557     0.128998     0.111449
           0.0926449    0.195847      0.134642     0.084457
           0.507474     0.109186      0.144895    -0.1903
           -0.0439681   -0.0443377    -0.338921    -0.364227
           0.307398     0.331592     -0.268884     0.0109789
           -0.00297574   0.0595918     0.112011     0.0919267
           -0.0601778    0.0736758    -0.022481     0.0812163
           0.247002     0.109472     -0.00769421  -0.383945
           0.126592     0.178113     -0.442779     0.222488
           0.0362439    0.0233627    -0.0142876   -0.190231
           -0.0222505    0.000279442   0.0915151    0.0481816
           0.0659491    0.159962     -0.310095     0.627363
           0.0115114    7.89468e-5   -0.0468136   -0.0231196
           -0.0121033   -0.0291859     0.0530118   -0.164273
           -0.00148844  -0.00618908    0.00372897  -0.0336697]
    test =cluster_ops[1]["aa"][((2, 1), (4, 1))][1,2,:,:]
    @test isapprox(abs.(ref), abs.(test), atol=1e-5)

    #for (ftrans,op) in cluster_ops[1]["aa"]
    #    println(ftrans, "  :  ", size(op))
    #    display(op[1,2,:,:])
    #end
    #for (ftrans,op) in cluster_ops[2]["AA"]
    #    println(ftrans, "  :  ", size(op))
    #    display(op[1,2,:,:])
    #end

    #for (ftrans,op) in cluster_ops[1]["Bb"]
    #    println(ftrans, "  :  ", size(op))
    #    display(op[1,2,:,:])
    #end
    ref = [ 0.46354     0.394923     -0.416441   -0.268368
           -0.199213    1.17757e-16   0.33466     0.590183
            0.356483    0.557353     -0.087595    0.503239
           -0.0456762  -0.18274      -0.0908969  -0.375945]
    test =cluster_ops[2]["Aa"][((1, 1), (1, 1))][1,2,:,:]
    @test isapprox(abs.(ref), abs.(test), atol=1e-5)

    ref = [ 0.370772    -0.0457363    0.0432994   0.264432    0.0339364   0.00585283
            0.510164    -1.20899e-5   0.17079     0.420218   -0.0816309   0.0788461
           -0.319318     0.237434     0.312764   -0.0502967  -0.433148    0.217789
           -0.356189     0.0186202   -0.0863461  -0.276715    0.0190341  -0.0341083
           -0.181094     0.0902997    0.0989763  -0.0682649  -0.155185    0.0736075
            0.00357835  -0.224128    -0.39496    -0.19786     0.456546   -0.251624]
    test =cluster_ops[1]["Bb"][((0, 2), (0, 2))][1,2,:,:]
    @test isapprox(abs.(ref), abs.(test), atol=1e-5)

    #for (ftrans,op) in cluster_ops[1]["Ab"]
    #    println(ftrans, "  :  ", size(op))
    #    display(op[1,2,:,:])
    #end
    ref = [-0.369922      0.0414125   0.0253631    0.00504684
           -0.41114      -0.0696597   0.12167     -0.00578434
           -0.0297536    -0.352606    0.298661    -0.0198366
           -8.22928e-5    0.280202   -0.218175     0.0411312
            0.0182344     0.0406758  -0.0132163    0.0410473
           -0.100527     -0.224248    0.072862    -0.226296
           -0.488761      0.025401   -0.00741915  -0.102493
            0.414459     -0.141834   -0.0148663   -0.119597
           -0.0286128     0.0514988  -0.0944393   -0.0892154
            0.065136      0.289568   -0.0509547    0.346402
            0.115824      0.25837    -0.0839488    0.260729
            0.16495       0.0135669  -0.0310352    0.0110329
            0.173482      0.136309   -0.339572    -0.318982
           -0.359717      0.324396   -0.232706    -0.0128002
           -0.167805     -0.0148749   0.0542493    0.0245395
            0.0704252     0.157099   -0.0510441    0.158533
           -0.0074902     0.0978906  -0.280457    -0.323457
            0.00575834   -0.0406361   0.0847375    0.0828477
           -0.0123439    -0.179338    0.0875409   -0.115216
           -0.145663      0.0755987   0.0108951    0.0881159
           -0.0112101    -0.113188   -0.131103    -0.380087
           -0.0261508     0.0420964  -0.261946    -0.377486
            0.000736687  -0.0341842   0.0390686    0.0156487
            0.00543429    0.0115098   0.0303084    0.0676828]
    test =cluster_ops[1]["Ab"][((2, 3), (1, 4))][1,2,:,:]
    @test isapprox(abs.(ref), abs.(test), atol=1e-5)

    display(cluster_ops[1]["AA"][((3,1),(1,1))][1,2,1,1])
    display(cluster_ops[1]["Aa"][((3,1),(3,1))][1,2,1,1])
    display(cluster_ops[1]["Ab"][((3,1),(2,2))][1,2,1,1])
    display(cluster_ops[1]["Ba"][((1,3),(2,2))][1,3,1,1])
    display(cluster_ops[1]["AB"][((3,3),(2,2))][1,3,1,2])
    display(cluster_ops[1]["AAa"][((3,2),(2,2))][1,2,3,1,2])
    display(cluster_ops[1]["Bbb"][((2,1),(2,2))][1,2,3,1,2])
    println(length(cluster_ops[1].data))
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG