text:
"""
We want to analyze only those data blocks that are active - so neither full nor empty

mainWorkQueue - basically we need to add to the one dimensional queue all of the  active blocks we find in first pass- crerating basic queue that will be processed by normals passes
of course during those normal passes some blocks will be added to the queue and some will be removed in order to mage it we would need additional helper structures
so we will have:
mainQuesCounter - telling us how many entries we have in work queue - we will divide this by number of thread blocks + some constant
    - to get some amount of the data blocks to be processed + tail queue that will be accessed in atomic way - but it can be accessed by all blocks
    so if some block will finish work before others it will start processing this tail queue
mainActiveCounterNow,mainActiveCounterNext - at first main queue will have only active blocks but progressively it will have more and more empty spots
    so we need to get the second counter that will keep track on the ramaining active blocks in curent iteration
        -mainActiveCounterNow - will be reduced every time thread block finish processing block - we will sync grid and start next iteration when it will reach 0
        -mainActiveCounterNext - will be increased every time we activate some block - will become the  mainActiveCounterNow in next iteration if it will reach 0 we will
            call it the end and finish kernel
"""
module ScheduleBlocksHD
using CUDA, ..CUDAGpuUtils, Logging,StaticArrays

"""
run periodically in normal pass - each thread block is reponsible to update workplan for itself
    - and add some additional to the general work queue so all thread blocks will get work for next iteration if possible
    we will get data about what blocks are active from analyzing metadata chunk that this block is responsible for
 """
function updateWorkScheduleNormalPass(metadata
                                     ,mataDataDims)

end#updateWorkScheduleNormalPass

"""
we have 3 dimensional data and in order to distribute work evenly and using consecutive memory chunks
    based on https://cs.calvin.edu/courses/cs/374/CUDA/CUDA-Thread-Indexing-Cheatsheet.pdf
"""
function getLinearIndexFrom3d()::UInt32
    blockId = blockIdx().x + blockIdx().y * gridDim.x
    + gridDim().x * gridDim().y * blockIdx().z

    return blockId * (blockDimX() * blockDimY() * blockDimZ())
    + threadIdxZ() * blockDimX() * blockDimY()
    + threadIdxY() * blockDimX() + threadIdxX()
end

end#ScheduleBlocksHD


text:
"""
    MinTimeZ <: AbstractSequencingAlgorithm

Algorithm that decides the sequence based on minimizing the time from the current
stop to the next stop, but subject to the additional constraint that in the final
sequence all zones are performed in order.
"""
struct MinTimeZ{T} <: AbstractSequencingAlgorithm
    algz::T
end

MinTimeZ() = MinTimeZ(MinTime())

function solve(alg::MinTimeZ, route::AbstractRoute, times::TravelTimes;
               initial_zone=zone(route.stops[1]),
               zones_list=zones(route, remove_station=false))

    route_id = name(route)
    T = times.travel_times

    algz = alg.algz

    # zonas que quedan por recorrer
    waiting_list = zones_list

    idx = findfirst(x -> x == initial_zone, waiting_list)

    sequence = Sequence(route_id, Vector{BareStop}())

    while length(waiting_list) > 0
        # filtrar solo las paradas de la zona actual
        actual_zone = waiting_list[idx]
        idx_stops = findall(s -> zone(s) == actual_zone, route.stops)

        # obtener las paradas y la matriz de tiempos asociada a esa zona
        subroute = BareRoute(route_id, route.stops[idx_stops])

        # indexamos la matriz de tiempos por nombres de parada
        name_stops = [name(s) for s in route.stops[idx_stops]]
        subtimes = view(T, name_stops, name_stops)

        # calcular la secuencia de acuerdo a los tiempos minimos
        initial_stop = _next_stop(alg, subroute)
        subsequence = solve(algz, subroute, subtimes, initial_stop=initial_stop)

        # append stops to the result
        for s in subsequence.stops
            push!(sequence, s)
        end

        deleteat!(waiting_list, idx)
    end

    return sequence
end

# simply return the first element of the list
function _next_stop(::MinTimeZ, subroute)
    return subroute.stops[1]
end

ext: jl


text:
function applyPauliExp!(sv::CuStateVec, theta::Float64, paulis::Vector{<:Pauli}, targets::Vector{Int32}, controls::Vector{Int32}, controlValues::Vector{Int32}=fill(one(Int32), length(controls)))
    cupaulis = CuStateVecPauli.(paulis)
    custatevecApplyExp(handle(), pointer(sv.data), sv.nbits, theta, cupaulis, targets, length(targets), controls, controlVals, length(controls))
    sv
end

function applyMatrix!(sv::CuStateVec, matrix::Union{Matrix, CuMatrix}, adjoint::Bool, targets::Vector{Int32}, controls::Vector{Int32}, controlValues::Vector{Int32}=fill(one(Int32), length(controls)))
    function bufferSize()
        out = Ref{Csize_t}()
        custatevecApplyMatrix_bufferSize(handle(), eltype(sv), sv.nbits, pointer(matrix), eltype(matrix), CUSTATEVEC_MATRIX_LAYOUT_COL, Int32(adjoint), length(targets), length(controls), compute_type(eltype(sv), eltype(matrix)), out)
        out[]
    end
    with_workspace(bufferSize) do buffer
        custatevecApplyMatrix(handle(), pointer(sv.data), eltype(sv), sv.nbits, pointer(matrix), eltype(matrix), CUSTATEVEC_MATRIX_LAYOUT_COL, Int32(adjoint), targets, length(targets), controls, length(controls), controlValues, compute_type(eltype(sv), eltype(matrix)), pointer(buffer), length(buffer))
    end
    sv
end

function expectation(sv::CuStateVec, matrix::Union{Matrix, CuMatrix}, basis_bits::Vector{Int32})
    function bufferSize()
        out = Ref{Csize_t}()
        custatevecExpectation_bufferSize(handle(), eltype(sv), sv.nbits, pointer(matrix), eltype(matrix), CUSTATEVEC_MATRIX_LAYOUT_COL, length(basis_bits), compute_type(eltype(sv), eltype(matrix)), out)
        out[]
    end
    expVal = Ref{Float64}()
    residualNorm = Ref{Float64}()
    with_workspace(bufferSize) do buffer
        custatevecExpectation(handle(), pointer(sv.data), eltype(sv), sv.nbits, expVal, Float64, residualNorm, pointer(matrix), eltype(matrix), CUSTATEVEC_MATRIX_LAYOUT_COL, basis_bits, length(basis_bits), compute_type(eltype(sv), eltype(matrix)), buffer, length(buffer))
    end
    return expVal[], residualNorm[]
end

function sample(sv::CuStateVec, sampled_bits::Vector{Int32}, shot_count)
    sampler = CuStateVecSampler(sv, UInt32(shot_count))
    bitstrings = Vector{custatevecIndex_t}(undef, shot_count)
    with_workspace(sampler.ws_size) do buffer
        custatevecSampler_preprocess(handle(), Ref(sampler.handle), buffer, length(buffer))
        custatevecSampler_sample(handle(), Ref(sampler.handle), bitstrings, sampled_bits, length(sampled_bits), rand(shot_count), shot_count, CUSTATEVEC_SAMPLER_OUTPUT_RANDNUM_ORDER)
    end
    return bitstrings
end


text:
SRC_DIR = joinpath(Base.source_dir(), "..", "src")
AUX_DIR = joinpath(SRC_DIR, "auxiliary")
OUT_DIR=joinpath(Base.source_dir(), "..", "output")
mkpath(OUT_DIR)

#include(joinpath(SRC_DIR, "BridgeSDEInference.jl"))
#using Main.BridgeSDEInference
#include(joinpath(SRC_DIR, "fitzHughNagumo.jl"))
#include(joinpath(SRC_DIR, "fitzHughNagumo_conjugateUpdt.jl"))

include(joinpath(SRC_DIR, "types.jl"))
include(joinpath(SRC_DIR, "vern7.jl"))
include(joinpath(SRC_DIR, "tsit5.jl"))
include(joinpath(SRC_DIR, "rk4.jl"))
include(joinpath(SRC_DIR, "ralston3.jl"))
include(joinpath(SRC_DIR, "priors.jl"))
include(joinpath(SRC_DIR, "guid_prop_bridge.jl"))

include(joinpath(SRC_DIR, "bounded_diffusion_domain.jl"))
include(joinpath(SRC_DIR, "radial_ornstein_uhlenbeck.jl"))
include(joinpath(SRC_DIR, "euler_maruyama_dom_restr.jl"))
include(joinpath(SRC_DIR, "prokaryotic_autoregulatory_gene_network.jl"))

include(joinpath(SRC_DIR, "random_walk.jl"))
include(joinpath(SRC_DIR, "blocking_schedule.jl"))
include(joinpath(SRC_DIR, "starting_pt.jl"))
include(joinpath(SRC_DIR, "mcmc.jl"))
include(joinpath(SRC_DIR, "path_to_wiener.jl"))

using StaticArrays
using Distributions # to define priors
using Random        # to seed the random number generator

# Let's generate the data
# -----------------------
using Bridge
include(joinpath(AUX_DIR, "data_simulation_fns.jl"))
Random.seed!(4)
# values taken from table 1 of Golithly and Wilkinson
θ₀ = [0.1, 0.7, 0.35, 0.2, 0.1, 0.9, 0.3, 0.1]
K = 10.0
Pˣ = Prokaryote(θ₀..., K)

x0, dt, T = ℝ{4}(7.0, 10.0, 4.0, 6.0), 1/5000, 50.0
tt = 0.0:dt:T
XX, _ = simulate_segment(ℝ{4}(0.0, 0.0, 0.0, 0.0), x0, Pˣ, tt)

using Plots

x₁ = [x[1] for x in XX.yy]
x₂ = [x[2] for x in XX.yy]
x₃ = [x[3] for x in XX.yy]
x₄ = [x[4] for x in XX.yy]

skip = 10
plot(XX.tt[1:skip:end], x₁[1:skip:end])
plot(XX.tt[1:skip:end], x₂[1:skip:end])
plot(XX.tt[1:skip:end], x₃[1:skip:end])
plot(XX.tt[1:skip:end], x₄[1:skip:end])

# Let's try proposal now
idx₁, idx₂ = 100001, 105001
t0, T = XX.tt[idx₁], XX.tt[idx₂]
x0, xT = XX.yy[idx₁], XX.yy[idx₂]

P̃ = ProkaryoteAux(θ₀..., K, t0, x0, T, xT, Val{(true, true, true, true)}())
tt = t0:dt:T

XX2, _ = simulate_segment(ℝ{4}(0.0, 0.0, 0.0, 0.0), x0, P̃, tt)

[x[1] for x in XX2.yy]

plot(XX.tt[idx₁:idx₂], x₁[idx₁:idx₂])
plot!(XX2.tt, [x[1] for x in XX2.yy])
plot(XX.tt[idx₁:idx₂], x₂[idx₁:idx₂])
plot!(XX2.tt, [x[2] for x in XX2.yy])
plot(XX.tt[idx₁:idx₂], x₃[idx₁:idx₂])
plot!(XX2.tt, [x[3] for x in XX2.yy])
plot(XX.tt[idx₁:idx₂], x₄[idx₁:idx₂])
plot!(XX2.tt, [x[4] for x in XX2.yy])

τ(t₀,T) = (x) ->  t₀ + (x-t₀) * (2-(x-t₀)/(T-t₀))
tt = τ(t0,T).(tt)
Σdiagel = 10^-4
Σ = SMatrix{4,4}(1I)*Σdiagel
L = SMatrix{4,4}(1I)

Pᵒ = GuidPropBridge(eltype(x0), tt, Pˣ, P̃, L, xT, Σ; solver=Vern7())

XX3, _ = simulate_segment(ℝ{4}(0.0, 0.0, 0.0, 0.0), x0, Pᵒ, tt)

for i in 1:1000
    XX3, _ = simulate_segment(ℝ{4}(0.0, 0.0, 0.0, 0.0), x0, Pᵒ, tt)
end

plot(XX.tt[idx₁:idx₂], x₁[idx₁:idx₂])
plot!(XX3.tt, [x[1] for x in XX3.yy])
plot(XX.tt[idx₁:idx₂], x₂[idx₁:idx₂])
plot!(XX3.tt, [x[2] for x in XX3.yy])
plot(XX.tt[idx₁:idx₂], x₃[idx₁:idx₂])
plot!(XX3.tt, [x[3] for x in XX3.yy])
plot(XX.tt[idx₁:idx₂], x₄[idx₁:idx₂])
plot!(XX3.tt, [x[4] for x in XX3.yy])


text:
module StateSpace
using Distributions
using ForwardDiff
import Distributions: mean, var, cov, rand

import Base: mean, filter, show
import StatsBase: loglikelihood

export
	# model types
	AbstractStateSpaceModel,
	AbstractSSM,
	AbstractGaussianSSM,
	LinearGaussianSSM,
	NonlinearGaussianSSM,
	NonlinearSSM,
	# filter types
	AbstractStateSpaceFilter,
	AbstractKalmanFilter,
	LinearKalmanFilter,
	NonlinearKalmanFilter,
	NonlinearFilter,
	KalmanFilter,
	KF,
	ExtendedKalmanFilter,
	EKF,
	UnscentedKalmanFilter,
	UKF,
	EnsembleKalmanFilter,
	EnKF,
	ParticleFilter,

	FilteredState,
	show,
	process_matrix,
	observation_matrix,
	predict,
	observe,
	update,
	update!,
	filter,
	smooth,
	loglikelihood,
	simulate

include("matrix_utils.jl")
include("model_types.jl")
include("filtered_states.jl")
include("filter_types.jl")
include("common.jl")
include("particle_filter.jl")

end # module


text:
module VulkanCore

include("LibVulkan.jl")
using .LibVulkan

const api = LibVulkan
include(joinpath(@__DIR__, "..", "docs", "vulkan_doc.jl"))

const vk = api
export vk

#### External definitions

# X11/X.h
const Window = UInt32 # unsigned long
const VisualID = UInt32 # unsigned long
# X11/Xlib.h
const Display = Cvoid # Opaque struct
# xcb.h
const xcb_connection_t = Cvoid # opaque struct
const xcb_window_t = UInt32
const xcb_visualid_t = UInt32
# Wayland
const wl_display = Cvoid # TODO: make opaque for now
const wl_surface = Cvoid # TODO: make opaque for now
# Mir
const MirConnection = Cvoid # TODO: make opaque for now
const MirSurface = Cvoid # TODO: make opaque for now
# Android
const ANativeWindow = Cvoid # TODO: make opaque for now
# Windows
const HINSTANCE = Ptr{Cvoid}
const HWND = Ptr{Cvoid}

end # module


text:
using GaussianDistributions

@inline tuplejoin(x) = x
@inline tuplejoin(x, y) = (x..., y...)
@inline tuplejoin(x, y, z...) = tuplejoin(tuplejoin(x, y), z...)

#NOTE it seems to be IMPERATIVE to define the function `num_non_hypo` PRIOR to
# running the @generated functions below due to the `world age problem`
# to this end in BridgeSDEInference all `examples` are imported before this file

"""
    φ_old(::Val{T}, args...)

Old way to compute the φ function appearing in the Girsanov formula and needed for
sampling from the full conditional distribution of the parameters (whose
indices are specified by the `Val`) conditional on the path,
observations and other parameters.
"""
@generated function φ_old(::Val{T}, args...) where T
    z = Expr(:tuple, (:(phi(Val($i), args...)) for i in 1:length(T) if T[i])...)
    return z
end

"""
    φ(::Val{T}, args...)

New way t0 compute the φ function appearing in the Girsanov formula and needed for
sampling from the full conditional distribution of the parameters (whose
indices are specified by the `Val`) conditional on the path,
observations and other parameters.
"""
@generated function φ(::Val{T}, t, x, P::S) where {T,S}
    data = Expr(:call, :tuplejoin, (:(phi(Val($i), t, x, P)) for i in 1:length(T) if T[i])...)
    mat = Expr(:call, SMatrix{num_non_hypo(S),sum(T)}, data)
    return mat
end

"""
    φᶜ(::Val{T}, args...)

Compute the φᶜ function appearing in the Girsanov formula. This function
complements φ.
"""
# by default define as a linear map
φᶜ(p, θ, args...) = φᶜlinear(p, args...) * θ

@generated function φᶜlinear(::Val{T}, t, x, P::S) where {T,S}
    data = Expr(:call, :tuplejoin, (:(phi(Val($i), t, x, P)) for i in 0:length(T) if i==0 || !T[i])...)
    mat = Expr(:call, SMatrix{num_non_hypo(S),length(T)-sum(T)+1}, data)
    return mat
end

"""
    conjugate_draw(θ, XX, PT, prior, updtIdx)

Draw from the full conditional distribution of the parameters whose indices are
specified by the object `updtIdx`, conditionally on the path given in container
`XX`, and conditionally on all other parameter values given in vector `θ`.
"""
function conjugate_draw(θ, XX, PT, prior, updtIdx)
    μ = mustart(updtIdx)
    𝓦 = μ*μ'
    ϑ = SVector(thetaex(updtIdx, θ))
    μ, 𝓦 = _conjugate_draw(ϑ, μ, 𝓦, XX, PT, updtIdx)

    Σ = inv(𝓦 + inv(Matrix(prior.Σ)))
    Σ = (Σ + Σ')/2 # eliminates numerical inconsistencies
    μ_post = Σ * (μ + Vector(prior.Σ\prior.μ))
    ϑ = rand(Gaussian(μ_post, Σ))
    move_to_proper_place(ϑ, θ, updtIdx)     # align so that dimensions agree
end
"""
    mustart(::Val{T}) where {T} = @SVector zeros(sum(T))

mustart initialize the vector μ used in the functions ϕ, ϕᶜ
"""
mustart(::Val{T}) where {T} = @SVector zeros(sum(T))
#NOTE already defined in coordinate_access.jl
#@generated function thetaex(::Val{T}, θ) where T
#    z = Expr(:tuple, 1.0, (:(θ[$i]) for i in 1:length(T) if  !T[i])...)
#    return z
#end

function _conjugate_draw_old(ϑ, μ, 𝓦, XX, PT, updtIdx)
    for X in XX
        for i in 1:length(X)-1
            φₜ = SVector(φ(updtIdx, X.tt[i], X.yy[i], PT))
            φᶜₜ = SVector(φᶜ(updtIdx, X.tt[i], X.yy[i], PT))
            dt = X.tt[i+1] - X.tt[i]
            dy = X.yy[i+1][2]-X.yy[i][2]
            μ = μ + (φₜ*dy - φₜ*dot(ϑ, φᶜₜ)*dt)/PT.σ^2 #safe to use a(X.tt[i], X.yy[i], PT)
            𝓦 = 𝓦 + (φₜ*φₜ'*dt)/PT.σ^2
        end
    end
    μ, 𝓦
end

"""
    hypo_a_inv(P, t, x)

Base definition, assumes no hypoellipticity and no closed form expression for
the inverse of `a`
"""
hypo_a_inv(P, t, x) = inv(a(t, x, P))
nonhypo(::Any, x) = x

function _conjugate_draw(ϑ, μ, 𝓦, XX, PT, updtIdx)
    for X in XX
        for i in 1:length(X)-1
            φₜ = φ(updtIdx, X.tt[i], X.yy[i], PT)
            φᶜₜ = φᶜ(updtIdx, ϑ, X.tt[i], X.yy[i], PT)
            Γ⁻¹ = hypo_a_inv(PT, X.tt[i], X.yy[i])
            dt = X.tt[i+1] - X.tt[i]
            dy = nonhypo(PT, X.yy[i+1])-nonhypo(PT, X.yy[i])
            μ = μ + φₜ'*Γ⁻¹*dy - φₜ'*Γ⁻¹*φᶜₜ*dt
            𝓦 = 𝓦 + φₜ'*Γ⁻¹*φₜ*dt
        end
    end
    μ, 𝓦
end


text:
struct GroupReserve <: AbstractReservesFormulation end
############################### Reserve Variables` #########################################
"""
This function checks if the variables for reserves were created
"""
function check_activeservice_variables(
    optimization_container::OptimizationContainer,
    contributing_services::Vector{<:PSY.Service},
)
    for service in contributing_services
        get_variable(optimization_container, PSY.get_name(service), typeof(service))
    end
    return
end

################################## Reserve Requirement Constraint ##########################
"""
This function creates the requirement constraint that will be attained by the apropriate services
"""
function service_requirement_constraint!(
    optimization_container::OptimizationContainer,
    service::SR,
    ::ServiceModel{SR, GroupReserve},
    contributing_services::Vector{<:PSY.Service},
) where {SR <: PSY.StaticReserveGroup}
    parameters = model_has_parameters(optimization_container)
    initial_time = model_initial_time(optimization_container)
    @debug initial_time
    time_steps = model_time_steps(optimization_container)
    name = PSY.get_name(service)
    constraint =
        get_constraint(optimization_container, make_constraint_name(REQUIREMENT, SR))
    use_slacks = get_services_slack_variables(optimization_container.settings)
    reserve_variables = [
        get_variable(optimization_container, PSY.get_name(r), typeof(r)) for
        r in contributing_services
    ]

    requirement = PSY.get_requirement(service)
    for t in time_steps
        resource_expression = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}()
        for reserve_variable in reserve_variables
            JuMP.add_to_expression!(resource_expression, sum(reserve_variable[:, t]))
        end
        if use_slacks
            resource_expression += slack_vars[t]
        end
        constraint[name, t] = JuMP.@constraint(
            optimization_container.JuMPmodel,
            resource_expression >= requirement
        )
    end

    return
end


text:
# Modelo de tanque de armazenamento com quatro alimentacoes e perdas
type Tanque4perdas
	Tanque4perdas()=begin
		PP2=outers.PP2
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"Steam tables"
			)),
			Fracao(Dict{Symbol,Any}(
				:Brief=>"Fracao massica de perdas"
			)),
			VazaoMassica(),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in1}",
				:PosX=>0.1,
				:PosY=>0.16
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in2}",
				:PosX=>0.1,
				:PosY=>0.35
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in3}",
				:PosX=>0.1,
				:PosY=>0.55
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in4}",
				:PosX=>0.1,
				:PosY=>0.75
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out}",
				:PosX=>1,
				:PosY=>0.45
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{perdas}",
				:PosX=>0.5,
				:PosY=>1
			)),
			[
				:(Fout.F = Fin1.F + Fin2.F + Fin3.F + Fin4.F - Fperda.F),
				:(Fsalida= Fin1.F + Fin2.F + Fin3.F + Fin4.F),
				:(Fperda.F=y*Fsalida),
				:(Fout.F * Fout.H = Fin1.F * Fin1.H + Fin2.F * Fin2.H + Fin3.F * Fin3.H + Fin4.F * Fin4.H - Fperda.F*Fperda.H),
				:([Fout.S,Fout.T] = PP2.propPH(Fout.P,Fout.H)),
				:(Fperda.S = Fout.S),
				:(Fperda.T = Fout.T),
				:(Fperda.P = Fout.P),
				:(Fperda.H = Fout.H),
			],
			[
				"","","","","","","","","",
			],
			[:PP2,],
			[:y,:Fsalida,:Fin1,:Fin2,:Fin3,:Fin4,:Fout,:Fperda,]
		)
	end
	PP2::DanaPlugin
	y::Fracao
	Fsalida::VazaoMassica
	Fin1::Corrente
	Fin2::Corrente
	Fin3::Corrente
	Fin4::Corrente
	Fout::Corrente
	Fperda::Corrente
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Tanque4perdas
function setEquationFlow(in::Tanque4perdas)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
end
function atributes(in::Tanque4perdas,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/tanque4perdas"
	drive!(fields,_)
	return fields
end
Tanque4perdas(_::Dict{Symbol,Any})=begin
	newModel=Tanque4perdas()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# # The Theta'-function of Erdos-Renyi graphs
# Let $q$ be an odd prime, and let $V = \mathrm{GF}(q)^3$ be a three dimensional
# vector space over the finite field of order $q$. The set of one dimensional
# subspaces, i.e. the projective plane, of $V$ is denoted by $\mathrm{PG}(2,q)$.
# There are $q^2+q+1$ such subspaces, which are the vertices of the Erdos-Renyi
# graph $\mathrm{ER}(q)$. Two vertices are adjacent if they are distinct and
# orthogonal, i.e. for two representing vectors $x$ and $y$ we have $x^Ty=0$. We
# are interested in the size of a maximum stable set of these graphs,
# specifically upper bounds for this value. Note that these are not the equally named
# random graphs.

using LinearAlgebra #hide
q = 7
PG2q = vcat([[0, 0, 1]],
[[0, 1, b] for b = 0:q-1],
[[1, a, b] for a = 0:q-1 for b = 0:q-1])
Adj = [x' * y % q == 0 && x != y for x in PG2q, y in PG2q]
size(Adj)

#
using Plots #hide
spy(Adj)

# ## The Theta'-function
# The Theta'-function $\vartheta(G)$ of a graph $G=(V,E)$ is such an upper bound, based on
# semidefinite programming:
#
# ```math
# \vartheta(G)\coloneqq \sup\{\langle X,J\rangle : \langle X,A\rangle = 0, X\succcurlyeq 0, X\geq 0\}.
# ```
#
# In vectorized standard form this is simply

N = length(PG2q)
C = ones(N^2)
A = vcat(vec(Adj)', vec(Matrix{Float64}(I, N, N))')
b = [0.0, 1.0];

# ## Determining the symmetry reduction
# We can now apply the Jordan reduction method to the problem.
# First, we need to determine an (optimal) admissible subspace.

using SDPSymmetryReduction
P = admPartSubspace(C, A, b, true)
P.n

using Test #src
@test P.n == 18 #src

# Running `admPartSubspace` returns a `Partition` object. `P.n` are the number of orbits (and thus
# variables), and `P.P` is a matrix with integer values from `1` trough `P.n`. Here, `P.P` looks like this
# (different color shades = different orbits):

heatmap(reverse(P.P, dims=1)) #hide

# Now we can block-diagonalize the algebra (numerically)
blkD = blockDiagonalize(P, true);
@test sort(blkD.blkSizes) == [2,2,2,2,3] #src

# ## Building the reduced SDP
# Since `blkD.blks[i]` is the block-diagonalized image of `P.P .== i`,
# we obtain the new, symmetry reduced SDP by
using SparseArrays
PMat = hcat([sparse(vec(P.P .== i)) for i = 1:P.n]...)
newA = A * PMat
newB = b
newC = C' * PMat;

# ## Solving the SDP with JuMP and CSDP

using JuMP, CSDP
m = Model(CSDP.Optimizer)

## Initialize variables corresponding parts of the partition P
## >= 0 because the original SDP-matrices are entry-wise nonnegative
x = @variable(m, x[1:P.n] >= 0)

@constraint(m, newA * x .== newB)
@objective(m, Max, newC * x)

psdBlocks = sum(blkD.blks[i] .* x[i] for i = 1:P.n)
for blk in psdBlocks
    if size(blk, 1) > 1
        @constraint(m, blk in PSDCone())
    else
        @constraint(m, blk .>= 0)
    end
end

optimize!(m)
#
termination_status(m)
#
objective_value(m)

@test objective_value(m) ≈ 15.743402681126568 atol = 5 #src


text:
import .AST: ASTMDP, ASTMetrics

"""
    print_metrics(planner)
    print_metrics(mdp::ASTMDP)
    print_metrics(metrics::ASTMetrics)

Print failure metrics including:
- iteration of first failure
- number of failures
- failure rate
"""
print_metrics(planner) = print_metrics(planner.mdp.metrics)
print_metrics(mdp::ASTMDP) = print_metrics(mdp.metrics)
function print_metrics(metrics::ASTMetrics)
    E = metrics.event

    if findfirst(E) === nothing
        @info "No failures recorded."
        return 0
    else
        first_failure = findfirst(E)
        num_evals = length(E)
        println("First failure: ", first_failure, " of ", num_evals)
        num_failures = sum(E)
        println("Number of failures: ", num_failures)
        failure_rate = sum(E)/length(E) * 100
        println("Failure rate: ", round(failure_rate, digits=5), "%")
        return failure_rate
    end
end


text:
using Kinetic
using KitBase.OffsetArrays, KitBase.ProgressMeter

cd(@__DIR__)
D = read_dict("mixture_shock.txt")
for key in keys(D)
    s = Symbol(key)
    @eval $s = $(D[key])
end

begin
    γ = heat_capacity_ratio(inK, 1)
    set = Setup(matter, case, space, flux, collision, nSpecies, interpOrder, limiter, boundary, cfl, maxTime)
    pSpace = PSpace1D(x0, x1, nx, nxg)
    ue0 = umin * sqrt(mi / me)
    ue1 = umax * sqrt(mi / me)
    vSpace = MVSpace1D(umin, umax, ue0, ue1, nu, vMeshType, nug)
    kne = knudsen * (me / mi)
    gas = Mixture(
        [knudsen, kne],
        mach,
        prandtl,
        inK,
        γ,
        mi,
        ni,
        me,
        ne,
    )
    wL, primL, hL, bL, bcL, wR, primR, hR, bR, bcR = ib_rh(gas.Ma, gas.γ, gas.K, mi, me, ni, ne, vSpace.u)
    ib = IB2F(wL, primL, hL, bL, bcL, wR, primR, hR, bR, bcR)

    ks = SolverSet(set, pSpace, vSpace, gas, ib, pwd())
end

begin
    ctr = OffsetArray{ControlVolume1D2F}(undef, axes(ks.pSpace.x, 1))
    face = Array{Interface1D2F}(undef, ks.pSpace.nx + 1)

    idx0 = (eachindex(pSpace.x) |> collect)[1]
    idx1 = (eachindex(pSpace.x) |> collect)[end]

    for i in eachindex(ctr)
        if i <= ks.pSpace.nx ÷ 2
            ctr[i] = ControlVolume1D2F(ks.pSpace.x[i], ks.pSpace.dx[i], ks.ib.wL, ks.ib.primL,
            ks.ib.hL, ks.ib.bL)
        else
            ctr[i] = ControlVolume1D2F(ks.pSpace.x[i], ks.pSpace.dx[i], ks.ib.wR, ks.ib.primR,
            ks.ib.hR, ks.ib.bR)
        end
    end

    face = Array{Interface1D2F}(undef, ks.pSpace.nx+1)
    for i=1:ks.pSpace.nx+1
        face[i] = Interface1D2F(ks.ib.wL, ks.ib.hL)
    end
end

begin
    iter = 0
    res = zeros(3)
    simTime = 0.0
    dt = Kinetic.timestep(ks, ctr, simTime)
    nt = Int(floor(ks.set.maxTime / dt))
end

res = zeros(3, 2)
begin
    Kinetic.reconstruct!(ks, ctr)
    Kinetic.evolve!(ks, ctr, face, dt)
    Kinetic.update!(ks, ctr, face, dt, res; bc=:fix)
end


text:
const SymStr = Union{Symbol,AbstractString}

import DataFrames: ByRow

iscolumn(tbl, var) = var_key(var) in propertynames(tbl)

var_key(var) = var # needed to handle, eg, `linewidth = 2`
var_lab(var::SymStr) = string(var)
var_trans(var::SymStr) = identity

var_key(var::Pair{<: SymStr, <: SymStr}) = first(var)
var_lab(var::Pair{<: SymStr, <: SymStr}) = string(last(var))
var_trans(var::Pair{<: SymStr, <: SymStr}) = identity

var_key(var::Pair{<: SymStr, <: Function}) = first(var)
var_lab(var::Pair{<: SymStr, <: Function}) = funname(var_trans(var)) * "(" * string(first(var)) * ")"
var_trans(var::Pair{<: SymStr, <: Function}) = last(var)

var_key(var::Pair{<: SymStr, <: Pair{<: Function, <: SymStr}}) = first(var)
var_lab(var::Pair{<: SymStr, <: Pair{<: Function, <: SymStr}}) = string(last(last(var)))
var_trans(var::Pair{<: SymStr, <: Pair{<: Function, <: SymStr}}) = first(last(var))

funname(f) = string(Symbol(f))
funname(f::ByRow) = funname(f.fun)

function get(tbl, var)
    f = var_trans(var)
    col = getproperty(tbl, var_key(var))
    if f isa ByRow
        (f.fun).(col)
    else
        f(col)
    end
end

max_stars_repo_name: greimel/TabularMakie

max_issues_repo_name: greimel/TabularMakie

max_forks_repo_name: greimel/TabularMakie


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using BenchmarkTools
using CUDA
using Oceananigans
using Benchmarks
using Plots
pyplot()
# Benchmark function

function benchmark_nonhydrostatic_model(Arch, FT, N)
    grid = RegularRectilinearGrid(FT, size=(N, N, N), extent=(1, 1, 1))
    model = NonhydrostaticModel(architecture=Arch(), grid=grid)

    time_step!(model, 1) # warmup

    trial = @benchmark begin
        @sync_gpu time_step!($model, 1)
    end samples=10

    return trial
end

# Benchmark parameters

Architectures = has_cuda() ? [CPU, GPU] : [CPU]
Float_types = [Float32, Float64]
Ns = [32, 64, 128, 256]

# Run and summarize benchmarks

print_system_info()
suite = run_benchmarks(benchmark_nonhydrostatic_model; Architectures, Float_types, Ns)

plot_num = length(Ns)
CPU_Float32 = zeros(Float64, plot_num)
CPU_Float64 = zeros(Float64, plot_num)
GPU_Float32 = zeros(Float64, plot_num)
GPU_Float64 = zeros(Float64, plot_num)
plot_keys = collect(keys(suite))
sort!(plot_keys, by = v -> [Symbol(v[1]), Symbol(v[2]), v[3]])

for i in 1:plot_num
    CPU_Float32[i] = mean(suite[plot_keys[i]].times) / 1.0e6
    CPU_Float64[i] = mean(suite[plot_keys[i+plot_num]].times) / 1.0e6
    GPU_Float32[i] = mean(suite[plot_keys[i+2plot_num]].times) / 1.0e6
    GPU_Float64[i] = mean(suite[plot_keys[i+3plot_num]].times) / 1.0e6
end

plt = plot(Ns, CPU_Float32, lw=4, label="CPU Float32", xaxis=:log2, yaxis=:log, legend=:topleft,
          xlabel="Nx", ylabel="Times (ms)", title="Nonhydrostatic Model Benchmarks: CPU vs GPU")
plot!(plt, Ns, CPU_Float64, lw=4, label="CPU Float64")
plot!(plt, Ns, GPU_Float32, lw=4, label="GPU Float32")
plot!(plt, Ns, GPU_Float64, lw=4, label="GPU Float64")
display(plt)
savefig(plt, "nonhydrostatic_times.png")

plt2 = plot(Ns, CPU_Float32./GPU_Float32, lw=4, xaxis=:log2, legend=:topleft, label="Float32",
            xlabel="Nx", ylabel="Speedup Ratio", title="Nonhydrostatic Model Benchmarks: CPU/GPU")
plot!(plt2, Ns, CPU_Float64./GPU_Float64, lw=4, label="Float64")
display(plt2)
savefig(plt2, "nonhydrostatic_speedup.png")

df = benchmarks_dataframe(suite)
sort!(df, [:Architectures, :Float_types, :Ns], by=(string, string, identity))
benchmarks_pretty_table(df, title="Nonhydrostatic model benchmarks")

if GPU in Architectures
    df_Δ = gpu_speedups_suite(suite) |> speedups_dataframe
    sort!(df_Δ, [:Float_types, :Ns], by=(string, identity))
    benchmarks_pretty_table(df_Δ, title="Nonhydrostatic model CPU to GPU speedup")
end


text:
using QuantumAlgebra
using QuantumAlgebra: δ, prodtuples, prodtuple, sumtuple, distribute_indices!
using Test

@testset "QuantumAlgebra.jl" begin
    for with_σpm in (false,true)
        QuantumAlgebra.use_σpm(with_σpm)

        @test QuantumAlgebra.using_σpm == with_σpm

        @test QuantumAlgebra.SpatialIndex(QuantumAlgebra.x) == QuantumAlgebra.x

        # all equal numbers should be equal scalars (ignore type)
        @test scal(0) == scal(0.0)
        @test scal(0) == scal(0im)
        @test scal(1.5) == scal(3//2)

        # test params and their complex conjugation
        @test param(:ω) == param(:ω,()) == param(:ω,'n') == param(:ω,'n',())
        @test param(:ω,'n',(:i)) == param(:ω,:i)
        @test param(:ω,'n',(:i,:j,2)) == param(:ω,:i,:j,2)
        @test_throws ErrorException param(:ω,'g')
        @test_throws MethodError param(:ω,2,:i,"a")
        @test_throws MethodError a("a")

        @test adjoint(param(:g,'r')) == param(:g,'r')
        @test adjoint(param(:g,'n')) == param(:g,'c')
        @test adjoint(param(:g,:i)) == param(:g,'c',:i)
        @test adjoint(param(:g,'c',(:i,1,:m))) == param(:g,:i,1,:m)

        tmp1 = param(:g,2)*param(:g,1)*param(:g,'c',3)*param(:b)*param(:a)*param(:d,'c',3)*param(:f,'r',1,:i)
        tmpc = param(:g,'c',2)*param(:g,'c',1)*param(:g,'n',3)*param(:b,'c')*param(:a,'c')*param(:d,'n',3)*param(:f,'r',(1,:i))
        tmp2 = param(:a)*param(:b)*param(:g,1)*param(:g,2)*param(:g,'c',3)*param(:d,'c',3)*param(:f,'r',(1,:i))
        @test tmp1 == tmp2
        @test tmp1' == tmpc

        @test a()*5 == 5*a() == scal(5)*a()
        @test a()+5 == 5+a() == scal(5)+a()
        @test 5-a() == scal(5)+scal(-1)*a()
        @test a()-5 == scal(-5)+a()

        @test_throws MethodError σx(:i,"a")
        @test σx(:i,:b) == σx((:i,:b))
        @test σx(:i,:j) != σx(:i,:b)

        @test σx() == σp() + σm()
        @test σx(:i) == σp(:i) + σm(:i)
        @test σx(:i,:j) == σp((:i,:j)) + σm((:i,:j))

        @test σy(:i) == scal(-1im)*(σp(:i) - σm(:i))
        @test σz(:i) == scal(2)*σp(:i)*σm(:i) - scal(1)

        for s=(σx,σy,σz)
            @test s(:i)*s(:i) == scal(1)
            @test s(:i)' == s(:i)
        end
        @test σx(:i)*σy(:i) == scal(1im)*σz(:i)
        @test σy(:i)*σx(:i) == scal(-1im)*σz(:i)
        @test σy()*σz() == scal(1im)*σx()
        @test σx(:i,:j)*σz(:i,:j) == scal(-1im)*σy(:i,:j)
        @test σx(:i,:j)*σz(:i,:j) != scal(-1im)*σy(:i,:k)

        @test a()' == adag()
        @test a(:m)*adag(:m) == adag(:m)*a(:m) + scal(1)
        @test (a(:m)*a(1))' == adag(:m)*adag(1)
        @test (a(1,2,:k)*a(:m))' == adag(1,2,:k)*adag(:m)

        @test fdag(:a)*f(:b) + f(:b)*fdag(:a) == δ(:a,:b)
        @test f(:a)*f(:b) + f(:b)*f(:a) == scal(0)
        @test fdag(:a)*fdag(:b) + fdag(:b)*fdag(:a) == scal(0)

        @test f()' == fdag()
        @test fdag()' == f()
        @test (fdag()*f())' == fdag()*f()
        @test comm(fdag(),f()) == -comm(f(),fdag())

        @test scal(1+1im) < scal(2-1im)
        @test scal(1+1im) > scal(1-1im)
        @test a(1) < σx(1)
        @test adag(1) < σx(1)
        @test !(a(1) < adag(1))
        @test adag(1) < a(1)
        @test a(5) < a(:i)
        @test !(adag(:m) < adag(2))
        @test fdag() < f()
        @test adag() < fdag()
        @test fdag(:i,:j) < fdag(:k)

        @test a(1) * (σy(1) * a(1))' == a(1) * (adag(1) * σy(1)) == adag(1)*a(1)*σy(1) + σy(1)

        @test ∑(:j,a(:j))*∑(:i,a(:i)) == ∑(:i,∑(:j,a(:i)*a(:j)))

        tmp = OpSumAnalytic(:i,adag(:i)*a(:i))
        @test tmp' == tmp
        @test a(:i)*OpSumAnalytic(:i,a(:i)) == OpSumAnalytic(:i_1,a(:i_1)*a(:i))
        @test adag(:n)*tmp == OpSumAnalytic(:i,adag(:n)*adag(:i)*a(:i))
        @test a(:n)   *tmp == OpSumAnalytic(:i,adag(:i)*a(:n)*a(:i)) + a(:n)
        @test param(:g,:i)*OpSumAnalytic(:i,a(:i)) == OpSumAnalytic(:i_1,param(:g,:i)*a(:i_1))
        @test param(:g,:i_1)*a(:i)*OpSumAnalytic(:i,a(:i)) == OpSumAnalytic(:i_2,param(:g,:i_1)*a(:i)*a(:i_2))
        @test param(:g,:n)*OpSumAnalytic(:i,a(:i)) == ∑(:i,param(:g,:n)*a(:i))

        @test Pr"gz_i,mu" == param(:gz,'r',(:i,:mu))
        @test Pc"gz_i,mu" == param(:gz,'n',(:i,:mu))
        @test Pc"gz_i,mu"' == param(:gz,'c',(:i,:mu))

        @test prodtuple(a(5)) == (a(5),)
        # tuples come out ordered!
        @test prodtuple(a(5)*a(4)) == (a(4),a(5))
        @test prodtuple(a(5)*a(4)) != (a(5),a(4))

        @test sumtuple(a(5)) == (a(5),)
        # tuples come out ordered!
        @test sumtuple(a(5)+a(4)) == (a(4),a(5))
        @test sumtuple(a(5)+a(4)) != (a(5),a(4))

        @test_throws ErrorException prodtuples(a(5)+a(4))
        # tuples come out ordered!
        if QuantumAlgebra.using_σpm
            tmp1 = scal(3)*param(:ω)*param(:g)*ExpVal(σp(:k))*σp(:k)*adag(5)*a(5)
            tmp2 = ( (scal(3),param(:g),param(:ω)), (ExpVal(σp(:k)),), (adag(5),a(5),σp(:k)) )
            @test prodtuples(tmp1) == tmp2
        else
            tmp1 = scal(3)*param(:ω)*param(:g)*ExpVal(σz(:k))*σz(:k)*adag(5)*a(5)
            tmp2 = ( (scal(3),param(:g),param(:ω)), (ExpVal(σz(:k)),), (adag(5),a(5),σz(:k)) )
            @test prodtuples(tmp1) == tmp2
        end

        @test δ(:i,:k)*a(:k) == a(:i)*δ(:k,:i)
        @test δ(:i,:k)*a(:k,:i) == a(:i,:i)*δ(:k,:i)
        @test δ(:i,:k)*δ(:i,:j) == δ(:k,:i)*δ(:j,:k)
        # k cannot be equal to 1 and 3 at the same time
        @test δ(1,:k)*δ(:k,3)*σx(:k) == scal(0)
        @test δ(1,:k)*δ(:k,1)*σx(:k) == δ(1,:k)*σx(1)

        @test comm(σx(5),σy(3)) == scal(0)
        @test comm(σx(5),σx(5)) == scal(0)
        @test comm(σx(1),σz(1)) == scal(-2im)*σy(1)
        @test comm(σx(:mu),σy(:muuu)) == scal(2im)*δ(:mu,:muuu)*σz(:mu)
        @test scal(1//2im)*comm(σx(:m),σy(:m)) == σz(:m)
        @test σx(:a)*σy(:a)*σz(:a) == scal(1im)

        @test comm(param(:g),a(5)+a(3)) == scal(0)
        @test comm(param(:g),a(5)*a(3)) == scal(0)
        @test comm(a(5)+a(3),param(:g)) == scal(0)
        @test comm(a(5)*a(3),param(:g)) == scal(0)

        @test comm(a(5)+a(3),adag(5)) == comm(a(5),adag(5))+ comm(a(3),adag(5))
        @test comm(a(5)+a(3),adag(5)*a(3)) == comm(a(5),adag(5)*a(3))+ comm(a(3),adag(5)*a(3))

        @test adag(2)*σy(:i) - scal(14)*param(:ω) == scal(-14)*param(:ω) + σy(:i)*adag(2)
        @test adag(2)*σy(:i) == σy(:i)*adag(2)

        @test σp(1) * σm(1) == scal(1//2)*σz(1) + scal(1//2)
        @test σm(1)*σp(1) == scal(1) - σp(1)*σm(1)
        @test σm(1)*σp(1) == σp(1)*σm(1) + comm(σm(1),σp(1))
        @test σz(1) == σp(1)*σm(1) - σm(1)*σp(1)

        @test comm(σp(1),σp(1)) == scal(0)
        @test comm(σp(:n),σm(:n)) == σz(:n)

        @test comm(a(1),   adag(1)*a(1)) == a(1)
        @test comm(adag(1),adag(1)*a(1)) == -adag(1)

        @test ExpVal(scal(3)) == scal(3)
        @test Corr(scal(3)) == scal(3)
        @test Corr(scal(3)+adag(2)) == scal(3) + Corr(adag(2))

        @test ascorr(scal(3)) == scal(3)
        @test ascorr(a(2)) == ExpVal(a(2))
        @test ascorr(scal(3)*param(:g)) == scal(3)*param(:g)
        @test ascorr(scal(3)*a(2)) == scal(3)*ExpVal(a(2))
        @test ascorr(scal(3)+adag(2)) == scal(3) + ExpVal(adag(2))

        @test ascorr(a(2)*a(2)) == Corr(a(2)*a(2)) + ExpVal(a(2))*ExpVal(a(2))
        @test ascorr(a(2)*a(:m))' == Corr(adag(2)*adag(:m)) + ExpVal(adag(:m))*ExpVal(adag(2))
        tmpas = a.(1:3)
        @test *(tmpas...) == a(1)*a(2)*a(3)
        tmpEVs = ExpVal.(tmpas)
        @test ascorr(*(tmpas...)) == Corr(*(tmpas...)) + *(tmpEVs...) + tmpEVs[1]*Corr(tmpas[2]*tmpas[3]) + tmpEVs[2]*Corr(tmpas[1]*tmpas[3]) + tmpEVs[3]*Corr(tmpas[1]*tmpas[2])

        @test a(1) < ascorr(a(1)*a(2)*a(3)*a(4))
        @test_throws ErrorException ascorr(a(1)*a(2)*a(3)*a(4)*a(5))

        if QuantumAlgebra.using_σpm
            @test ascorr(scal(-1)*param(:g,'r',1)*σp(1)) == -param(:g,'r',1)*ExpVal(σp(1))
            @test ascorr(OpSumAnalytic(:i,σp(:i)*σm(:n))) == OpSumAnalytic(:i,Corr(σp(:i)*σm(:n))) + OpSumAnalytic(:i,ExpVal(σp(:i))*ExpVal(σm(:n)))
        else
            @test ascorr(scal(-1)*param(:g,'r',1)*σz(1)) == -param(:g,'r',1)*ExpVal(σz(1))
            @test ascorr(OpSumAnalytic(:i,σy(:i)*σy(:n))) == OpSumAnalytic(:i,Corr(σy(:i)*σy(:n))) + OpSumAnalytic(:i,ExpVal(σy(:i))*ExpVal(σy(:n))) - ExpVal(σy(:n))*ExpVal(σy(:n))
        end

        @test CorrOrExp(a(5)) == ExpVal(a(5))
        @test CorrOrExp(a(5)*a(:i)) == Corr(a(5)*a(:i))

        H = ∑(:i,param(:ω,'r',:i)*adag(:i)*a(:i))
        @test comm(a(:i),H) == param(:ω,'r',:i)*a(:i)
        @test comm(a(:n),H) == param(:ω,'r',:n)*a(:n)
        @test comm(H,a(:n)) == -param(:ω,'r',:n)*a(:n)
        @test comm(adag(:n),H) == -param(:ω,'r',:n)*adag(:n)
        @test comm(adag(:n)*a(:m),H) == (param(:ω,'r',:m)-param(:ω,'r',:n))*adag(:n)*a(:m)

        @test a()*H == ∑(:i,param(:ω,'r',:i)*adag(:i)*a(:i)*a())
        @test a(:k)*H == param(:ω,'r',:k)*a(:k) + ∑(:i,param(:ω,'r',:i)*adag(:i)*a(:i)*a(:k))
        HH = ∑(:i,param(:ω,'r',:i,:i)*adag(:i,:i)*a(:i,:i))
        @test a(:k,:k)*HH == param(:ω,'r',:k,:k)*a(:k,:k) + ∑(:i,param(:ω,'r',:i,:i)*adag(:i,:i)*a(:i,:i)*a(:k,:k))

        @test Avac(H) == scal(0)
        @test vacA(H) == scal(0)
        @test vacA(adag(3)*σp(1)*σm(1)) == scal(0)
        @test vacA(fdag(:n)) == scal(0)
        @test vacA(f(:n)) == f(:n)
        @test vacA(a(:n)) == a(:n)
        @test Avac(fdag(:n)) == fdag(:n)
        @test Avac(f(:n)) == scal(0)
        @test Avac(a(3)*σp(1)*σm(1)) == scal(0)
        @test Avac(σm(1)*σp(1)) == scal(1)
        @test Avac(σp(1)*σm(1)) == scal(0)
        @test Avac(σp(1)) == σp(1)
        @test vacA(σm(1)) == σm(1)
        if QuantumAlgebra.using_σpm
            @test Avac(σm(1)) == scal(0)
            @test vacA(σp(1)) == scal(0)
        else
            @test Avac(σx(1)) == vacA(σx(1)) == σx(1)
        end
        @test vacExpVal(σx(1)) == scal(0)
        @test vacExpVal(σp(1)) == scal(0)
        @test vacExpVal(OpSumAnalytic(:i,σp(:i))) == scal(0)
        @test vacExpVal(σp(:i)*σm(:k)) == scal(0)

        @test a(:n)*adag(:n)*a(:n)*adag(:n) == scal(1) + scal(3)*adag(:n)*a(:n) + adag(:n)*adag(:n)*a(:n)*a(:n)

        S = scal(1/√(2*6))*adag(:n)*adag(:n)*adag(:n) + scal(1/√2)*adag(:m)
        for (A,val) in [(scal(1),scal(1)),
                        (adag(:n)*a(:n),scal(1.5) + 0.5 * δ(:n,:m)),
                        (adag(:n)*adag(:n)*a(:n)*a(:n),scal(3))]
            @test vacExpVal(A,S) ≈ val
        end

        tmp = scal(1+2im)*OpSumAnalytic(:i,a(:i)*adag(:i)*ascorr(adag(:n)*a(:m)))
        @test latex(scal(1+2im)) == "(1+2i)"
        @test latex(scal(1+2im//5)) == "\\left(1+\\frac{2}{5}i\\right)"

        tmp = OpSumAnalytic(:i,ascorr(adag(:n)*a(:i)))
        tmplatex = "\\sum_{i}\\langle a_{n}^\\dagger a_{i} \\rangle_{c} + \\sum_{i}\\langle a_{n}^\\dagger \\rangle \\langle a_{i} \\rangle"
        @test latex(tmp) == tmplatex
        @test ascorr(tmp) == tmp
        @test sprint(show,"text/latex",tmp) == "\$$(tmplatex)\$"
        if QuantumAlgebra.using_σpm
            @test latex(σp()) == "\\sigma^+"
        else
            @test latex(σz()) == "\\sigma_{z}"
        end

        inds = [:a,:b,:c,:d,:e,:f,:g,:h,:i,:j,:k,:l,:m,:n]
        tmp1 = param(:ω,:y)*a(1)*adag(1)*a(3)*adag(4)*ExpVal(a(5))*Corr(adag(5)*a(9))
        tmp2 = param(:ω,:a)*ExpVal(a(:b))*Corr(adag(:c)*a(:d))*adag(:e)*a(:f) + param(:ω,:g)*ExpVal(a(:h))*Corr(adag(:i)*a(:j))*adag(:k)*adag(:l)*a(:m)*a(:n)
        @test distribute_indices!(copy(inds),tmp1) == tmp2
        if QuantumAlgebra.using_σpm
            tmp1 = a(1,:n)*adag()*σm(1,:n)*σp()
            @test distribute_indices!(copy(inds),tmp1) == adag()*a(:a,:b)*σp()*σm(:c,:d)
        else
            tmp1 = a(1,:n)*adag()*σz(1,:n)*σy()
            @test distribute_indices!(copy(inds),tmp1) == adag()*a(:a,:b)*σy()*σz(:c,:d)
        end

        @test_throws MethodError distribute_indices!(copy(inds),OpSumAnalytic(:i,a(:i)))
        @test_throws ArgumentError distribute_indices!([:a,:b],tmp1)

        @test QuantumAlgebra.exchange_inds(adag(:j)*a(:k),:k,:j) == adag(:k)*a(:j)
        @test QuantumAlgebra.extindices(∑(:i,adag(:i)*a(:k))) == [:k]
        @test QuantumAlgebra.symmetric_index_nums(adag(:i)*adag(:j)*a(:k)*a(:l)) == [2,2]

        @test string(OpSumAnalytic(:i,a(:i)) * adag(:n)) == "1 + ∑_i a†(n) a(i)"
        if QuantumAlgebra.using_σpm
            @test string(a(5)*adag(5)*σp(3)*ascorr(adag(5,:i)*a(5))) == "⟨a†(5i)⟩ ⟨a(5)⟩ σ+(3) + ⟨a†(5i) a(5)⟩c σ+(3) + ⟨a†(5i)⟩ ⟨a(5)⟩ a†(5) a(5) σ+(3) + ⟨a†(5i) a(5)⟩c a†(5) a(5) σ+(3)"
        else
            @test string(a(5)*adag(5)*σz(3)*ascorr(adag(5,:i)*a(5))) == "⟨a†(5i)⟩ ⟨a(5)⟩ σz(3) + ⟨a†(5i) a(5)⟩c σz(3) + ⟨a†(5i)⟩ ⟨a(5)⟩ a†(5) a(5) σz(3) + ⟨a†(5i) a(5)⟩c a†(5) a(5) σz(3)"
        end
    end
end


text:
module Bijectors

#=
  NOTE: Codes below are adapted from
  https://github.com/brian-j-smith/Mamba.jl/blob/master/src/distributions/transformdistribution.jl
  The Mamba.jl package is licensed under the MIT License:
  > Copyright (c) 2014: Brian J Smith and other contributors:
  >
  > https://github.com/brian-j-smith/Mamba.jl/contributors
  >
  > Permission is hereby granted, free of charge, to any person obtaining
  > a copy of this software and associated documentation files (the
  > "Software"), to deal in the Software without restriction, including
  > without limitation the rights to use, copy, modify, merge, publish,
  > distribute, sublicense, and/or sell copies of the Software, and to
  > permit persons to whom the Software is furnished to do so, subject to
  > the following conditions:
  >
  > The above copyright notice and this permission notice shall be
  > included in all copies or substantial portions of the Software.
  >
  > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  > EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  > MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
  > IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
  > CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
  > TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
  > SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#

using Reexport, Requires
@reexport using Distributions
using StatsFuns
using LinearAlgebra
using MappedArrays
using Base.Iterators: drop
using LinearAlgebra: AbstractTriangular
import Functors
import NonlinearSolve
import ChainRulesCore

export  TransformDistribution,
        PositiveDistribution,
        UnitDistribution,
        SimplexDistribution,
        PDMatDistribution,
        link,
        invlink,
        logpdf_with_trans,
        isclosedform,
        transform,
        forward,
        logabsdetjac,
        logabsdetjacinv,
        Bijector,
        ADBijector,
        Inverse,
        Composed,
        compose,
        Stacked,
        stack,
        Identity,
        bijector,
        transformed,
        UnivariateTransformed,
        MultivariateTransformed,
        logpdf_with_jac,
        logpdf_forward,
        PlanarLayer,
        RadialLayer,
        CouplingLayer,
        InvertibleBatchNorm

if VERSION < v"1.1"
    using Compat: eachcol
end

const DEBUG = Bool(parse(Int, get(ENV, "DEBUG_BIJECTORS", "0")))
_debug(str) = @debug str

_eps(::Type{T}) where {T} = T(eps(T))
_eps(::Type{Real}) = eps(Float64)
_eps(::Type{<:Integer}) = eps(Float64)

function _clamp(x, a, b)
    T = promote_type(typeof(x), typeof(a), typeof(b))
    ϵ = _eps(T)
    clamped_x = ifelse(x < a, convert(T, a), ifelse(x > b, convert(T, b), x))
    DEBUG && _debug("x = $x, bounds = $((a, b)), clamped_x = $clamped_x")
    return clamped_x
end

function mapvcat(f, args...)
    out = map(f, args...)
    init = vcat(out[1])
    return reshape(reduce(vcat, drop(out, 1); init = init), size(out))
end

function maphcat(f, args...)
    out = map(f, args...)
    init = reshape(out[1], :, 1)
    return reduce(hcat, drop(out, 1); init = init)
end
function eachcolmaphcat(f, x1, x2)
    out = [f(x1[:,i], x2[i]) for i in 1:size(x1, 2)]
    init = reshape(out[1], :, 1)
    return reduce(hcat, drop(out, 1); init = init)
end
function eachcolmaphcat(f, x)
    out = map(f, eachcol(x))
    init = reshape(out[1], :, 1)
    return reduce(hcat, drop(out, 1); init = init)
end
function sumeachcol(f, x1, x2)
    # Using a view below for x1 breaks Tracker
    return sum(f(x1[:,i], x2[i]) for i in 1:size(x1, 2))
end

# Distributions

link(d::Distribution, x) = bijector(d)(x)
invlink(d::Distribution, y) = inv(bijector(d))(y)
function logpdf_with_trans(d::Distribution, x, transform::Bool)
    if ispd(d)
        return pd_logpdf_with_trans(d, x, transform)
    elseif isdirichlet(d)
        l = logpdf(d, x .+ eps(eltype(x)))
    else
        l = logpdf(d, x)
    end
    if transform
        return l - logabsdetjac(bijector(d), x)
    else
        return l
    end
end

## Univariate

const TransformDistribution = Union{
    T,
    Truncated{T},
} where T <: ContinuousUnivariateDistribution
const PositiveDistribution = Union{
    BetaPrime, Chi, Chisq, Erlang, Exponential, FDist, Frechet, Gamma, InverseGamma,
    InverseGaussian, Kolmogorov, LogNormal, NoncentralChisq, NoncentralF, Rayleigh, Weibull,
}
const UnitDistribution = Union{Beta, KSOneSided, NoncentralBeta}

function logpdf_with_trans(d::UnivariateDistribution, x, transform::Bool)
    if transform
        return map(x -> logpdf(d, x), x) - logabsdetjac(bijector(d), x)
    else
        return map(x -> logpdf(d, x), x)
    end
end

## Multivariate

const SimplexDistribution = Union{Dirichlet}
isdirichlet(::SimplexDistribution) = true
isdirichlet(::Distribution) = false

###########
# ∑xᵢ = 1 #
###########

function link(
    d::Dirichlet,
    x::AbstractVecOrMat{<:Real},
    proj::Bool = true,
)
    return SimplexBijector{proj}()(x)
end

function link_jacobian(
    d::Dirichlet,
    x::AbstractVector{T},
    proj::Bool = true,
) where {T<:Real}
    return jacobian(SimplexBijector{proj}(), x)
end

function invlink(
    d::Dirichlet,
    y::AbstractVecOrMat{<:Real},
    proj::Bool = true
)
    return inv(SimplexBijector{proj}())(y)
end
function invlink_jacobian(
    d::Dirichlet,
    y::AbstractVector{T},
    proj::Bool = true
) where {T<:Real}
    return jacobian(inv(SimplexBijector{proj}()), y)
end

## Matrix

#####################
# Positive definite #
#####################

const PDMatDistribution = Union{MatrixBeta, InverseWishart, Wishart}
ispd(::Distribution) = false
ispd(::PDMatDistribution) = true

function logpdf_with_trans(
    d::MatrixDistribution,
    X::AbstractArray{<:AbstractMatrix{<:Real}},
    transform::Bool,
)
    return map(X) do x
        logpdf_with_trans(d, x, transform)
    end
end
function pd_logpdf_with_trans(d, X::AbstractMatrix{<:Real}, transform::Bool)
    T = eltype(X)
    Xcf = cholesky(X, check = false)
    if !issuccess(Xcf)
        Xcf = cholesky(X + max(eps(T), eps(T) * norm(X)) * I)
    end
    lp = getlogp(d, Xcf, X)
    if transform && isfinite(lp)
        U = Xcf.U
        d = dim(d)
        lp += sum((d .- (1:d) .+ 2) .* log.(diag(U)))
        lp += d * log(T(2))
    end
    return lp
end
function getlogp(d::MatrixBeta, Xcf, X)
    n1, n2 = params(d)
    p = dim(d)
    return ((n1 - p - 1) / 2) * logdet(Xcf) + ((n2 - p - 1) / 2) * logdet(I - X) + d.logc0
end
function getlogp(d::Wishart, Xcf, X)
    return 0.5 * ((d.df - (dim(d) + 1)) * logdet(Xcf) - tr(d.S \ X)) + d.logc0
end
function getlogp(d::InverseWishart, Xcf, X)
    Ψ = Matrix(d.Ψ)
    return -0.5 * ((d.df + dim(d) + 1) * logdet(Xcf) + tr(Xcf \ Ψ)) + d.logc0
end

include("utils.jl")
include("interface.jl")
include("chainrules.jl")

# Broadcasting here breaks Tracker for some reason
maporbroadcast(f, x::AbstractArray{<:Any, N}...) where {N} = map(f, x...)
maporbroadcast(f, x::AbstractArray...) = f.(x...)

# optional dependencies
function __init__()
    @require LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" begin
        function maporbroadcast(f, x1::LazyArrays.BroadcastArray, x...)
            return copy(f.(x1, x...))
        end
        function maporbroadcast(f, x1, x2::LazyArrays.BroadcastArray, x...)
            return copy(f.(x1, x2, x...))
        end
        function maporbroadcast(f, x1, x2, x3::LazyArrays.BroadcastArray, x...)
            return copy(f.(x1, x2, x3, x...))
        end
    end
    @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("compat/forwarddiff.jl")
    @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" include("compat/tracker.jl")
    @require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" include("compat/zygote.jl")
    @require ReverseDiff="37e2e3b7-166d-5795-8a7a-e32c996b4267" include("compat/reversediff.jl")
    @require DistributionsAD="ced4e74d-a319-5a8a-b0ac-84af2272839c" include("compat/distributionsad.jl")
end

end # module


text:
# # In-place mutation with Referenceables.jl

using CUDA
using Folds
using FoldsCUDA
using Referenceables: referenceable

using Test                                                             #src

# ## Simple mutation

function increment_folds!(xs)
    Folds.foreach(referenceable(xs)) do x
        x[] += 1
    end
    return xs
end

if has_cuda_gpu()
    xs = CuArray(1:5)
else
    xs = Array(1:5)
end

@test begin
    collect(increment_folds!(xs))
end == 2:6

# This can also be written with FLoops.jl:

using FLoops

function increment_floops!(xs, ex = nothing)
    @floop ex for x in referenceable(xs)
        x[] += 1
    end
    return xs
end

@test begin
    collect(increment_floops!(xs))
end == 3:7

# ## Fusing reduction and mutationg

# Computing `sum(f, xs)` and `f.(xs)` in one go:

function mutation_with_folds(f, xs)
    ys = similar(xs)
    s = Folds.sum(zip(referenceable(ys), xs)) do (r, x)
        r[] = y = f(x)
        return y
    end
    return s, ys
end
nothing  # hide
#-

if has_cuda_gpu()
    xs = CuArray(1:5)
else
    xs = Array(1:5)
end

s, ys = mutation_with_folds(x -> x^2, xs)
@test begin
    s
end == 55
#-

@test begin
    collect(ys)
end == collect(xs) .^ 2
#-

# An equilvalent computaton with FLoops.jl:

using FLoops

function mutation_with_floop(f, xs)
    ys = similar(xs)
    z = zero(eltype(ys))
    @floop for (r, x) in zip(referenceable(ys), xs)
        r[] = y = f(x)
        @reduce(s = z + y)
    end
    return s, ys
end
nothing  # hide
#-

@assert mutation_with_folds(x -> x^2, xs) == (s, ys)
#-


text:
using Oceananigans

include("ConvergenceTests/ConvergenceTests.jl")

run_xy = ConvergenceTests.ForcedFlowFixedSlip.setup_and_run_xy
run_xz = ConvergenceTests.ForcedFlowFixedSlip.setup_and_run_xz

# Run 4 simulations:
Nx = [16, 32, 64, 128]
stop_time = 0.01
h = π / maximum(Nx)

Δt = 0.01 .* h^2 # Initial guess

# Modify so simulations stop on a dime
stop_iteration = round(Int, stop_time / Δt)
Δt = stop_time / stop_iteration

# Run (x, y) and (x, z) simulations
for N in Nx
    run_xy(Nx=N, Δt=Δt, stop_iteration=stop_iteration)
    run_xz(Nx=N, Δt=Δt, stop_iteration=stop_iteration)
end


text:
"""
The `ViscousFlow` module is here

"""
module ViscousFlow

using DocStringExtensions
using Reexport

@reexport using CartesianGrids
@reexport using ImmersedLayers
@reexport using RigidBodyTools
@reexport using ConstrainedSystems
@reexport using GridUtilities

using LinearAlgebra
using SparseArrays
using UnPack

export NavierStokes, PulseParams, PointForce, SpatialDGaussian,
       ExternalFlow, InternalFlow, ExternalInternalFlow,
       setstepsizes, timestep, timerange, newstate, flowside,
       update_immersion_operators!,
       vorticity, velocity, velocity!, streamfunction, streamfunction!,
       scalarpotential, scalarpotential!, convective_derivative, convective_derivative!,
       pressure, traction, force, moment, pressurejump

abstract type PointMotionType end
abstract type StaticPoints <: PointMotionType end
abstract type MovingPoints <: PointMotionType end

abstract type FreestreamType end
abstract type StaticFreestream <: FreestreamType end
abstract type VariableFreestream <: FreestreamType end

abstract type FlowSide end
abstract type ExternalFlow <: FlowSide end
abstract type InternalFlow <: FlowSide end
abstract type ExternalInternalFlow <: FlowSide end

include("utils/pulses.jl")
include("types.jl")
include("navier_stokes.jl")
include("directmotion.jl")
include("plot_recipes.jl")

end


text:
""" Struct for 2D spectra specialized for EXPRES and its additional data products.
    Drafted, but not yet folded into examples, as still using BasicSpectra2D instead.
    For now extra info is showing up in metadata.
    """
struct Spectra2DEXPRES{T1<:Real,T2<:Real,T3<:Real,
     AA1<:AbstractArray{T1,2},AA2<:AbstractArray{T2,2},AA3<:AbstractArray{T3,2},
     InstT<:AbstractInstrument
    } <: AbstractSpectra2D
    λ::AA1
    flux::AA2
    var::AA3
    blaze::AA2
    continuum::AA2
    tellurics::AA2
    inst::InstT
    metadata::Dict{Symbol,Any} # Dict{Symbol,Any}
end

function Spectra2DEXPRES(λ::A1, flux::A2, var::A3, inst::InstT; blaze::A2, continuum::A2, tellurics::A2,
     metadata::Dict{Symbol,Any} = Dict{Symbol,Any}() ) where {
     T1<:Real, T2<:Real, T3<:Real, A1<:AbstractArray{T1,2}, A2<:AbstractArray{T2,2}, A3<:AbstractArray{T3,2},
     InstT<:AbstractInstrument  }
    @assert size(λ) == size(flux)
    @assert size(λ) == size(var)
    @assert 1 <= size(λ,1) <= max_pixel_in_order(inst)-min_pixel_in_order(inst)+1
    @assert 1 <= size(λ,2) <= max_order(inst)-min_order(inst)+1
    Spectra2DEXPRES{eltype(λ),eltype(flux),eltype(var),typeof(λ),typeof(flux),typeof(var),typeof(inst)}(λ,flux,var,blaze,continuum,tellurics,inst,metadata)
end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
# here we define and document the prototypes that are implemented

"""
`AbstractAtoms`

theabstractsupertype for storing atomistic
configurations. A basic implementation might simply store a list of positions.
"""
abstract type AbstractAtoms end

"""
abstract supertype of neighbourlists.

The standard contructor should be of the form
```
   NeighbourList(at::AbstractAtoms, rcut)
```
where `rcut` should be either a scalar or a vector of cut-offs.
"""
abstract type AbstractNeighbourList end

abstract type AbstractConstraint end

"""
`Dofs`: dof vector
"""
const Dofs = Vector{Float64}

"""
`AbstractCalculator`: theabstractsupertype of calculators. These
store model information, and are linked to the implementation of energy,
forces, and so forth.
"""
abstract type AbstractCalculator end

"""
`Preconditioner`:abstractbase type for preconditioners

Preconditioners need to implement the following three functions:
* `update!(P, at::AbstractAtoms)`
* `A_mul_B!(out::Dof, P, x::Dof)`
* `A_ldiv_B!(out::Dof, P, f::Dof)`
"""
abstract type Preconditioner end # <: AbstractMatrix{Float64} end

"""
`macro protofun(fsig::Expr)`

generates "function prototypes" which throw an error message
with a little extra information in case a certain function has not been
implemented. This is in some way duplicating default Julia behaviour, but it
has the advantage that we can write documention for non-existing functions and
it also emphasizes that this specific function is part of theabstractatoms
interface.

## Usage

```julia
"Returns number of atoms"
@protofun length(::AbstractAtoms)
```

This will create a function `length(at::AbstractAtoms)`; if somebody
calls `length(at)` on a concrete sub-type of `AbstractAtoms` for which
`length` hasn't been implemented, it will throw an error.
"""
macro protofun(fsig::Expr)
    @assert fsig.head == :call
    fname = fsig.args[1]
    argnames = Any[]
    for idx in 2:length(fsig.args)
        arg = fsig.args[idx]
        if isa(arg, Expr) && arg.head == :kw
            arg = arg.args[1]
        end
        if isa(arg, Symbol)
            push!(argnames, arg)
        elseif isa(arg, Expr) && arg.head == :(::)
            if length(arg.args) != 2
                @gensym s
                insert!(arg.args, 1, s)
            end
            push!(argnames, arg.args[1])
        end
    end
    body = quote
        error(string("JuLIP: ", $fname,
                     ($([:(typeof($(esc(arg)))) for arg in argnames]...),),
                     " ) has no implementation."))
    end
    return Expr(:function, esc(fsig), body)
end

# function defined primarily on AbstractAtoms
export AbstractAtoms,
      positions, get_positions, set_positions!,
      momenta, get_momenta, set_momenta!,
      masses, get_masses, set_masses,
      cell, get_cell, set_cell!, is_cubic,
      pbc, get_pbc, set_pbc!,
      set_data!, get_data, has_data,
      set_calculator!, calculator, get_calculator,
      set_constraint!, constraint, get_constraint,
      neighbourlist, cutoff,
      defm, set_defm!

# length is used for several things
import Base: length, A_ldiv_B!, A_mul_B!, gradient

export AbstractCalculator,
      energy, potential_energy, forces, gradient, hessian,
      site_energies,
      stress, virial, site_virials

export AbstractNeighbourList,
       sites, bonds

export AbstractConstraint, NullConstraint,
         position_dofs, set_position_dofs!,
         momentum_dofs, set_momentum_dofs!,
         dofs, set_dofs!

export Preconditioner, preconditioner

# TODO: probably rename Preconditioner to AbstractPreconditioner and
#       AMGPrecon to Preconditioner
# TODO: iterators for angles, dihedrals

"Return number of atoms"
@protofun length(::AbstractAtoms)

"Return copy of positions of all atoms as a `3 x N` array."
@protofun positions(::AbstractAtoms)

"alias for `positions`"
get_positions = positions

"Set positions of all atoms"
@protofun set_positions!(::AbstractAtoms, ::Any)

set_positions!(at::AbstractAtoms, p::Matrix) = set_positions!(at, vecs(p))
set_positions!(at::AbstractAtoms, x, y, z) = set_positions!(at, [x'; y'; z'])

xyz(at::AbstractAtoms) = xyz(positions(at))

"Return copy of momenta of all atoms as a `3 x N` array."
@protofun momenta(::AbstractAtoms)

"alias for `momenta`"
get_momenta = momenta

"Set momenta of all atoms as a `3 x N` array."
@protofun set_momenta!(::AbstractAtoms, ::Any)

set_momenta!(at::AbstractAtoms, p::Matrix) = set_momenta!(at, vecs(p))

@protofun masses(::AbstractAtoms)
get_masses = masses
@protofun set_masses!(::AbstractAtoms, ::Any)

"get computational cell (the rows are the lattice vectors)"
@protofun cell(::AbstractAtoms)

"alias for `cell`"
get_cell = cell

"set computational cell"
@protofun set_cell!(::AbstractAtoms, ::AbstractMatrix)

# TODO: deprecate these!
"deformation matrix; `defm(at) = cell(at)'`"
defm(at::AbstractAtoms) = JMat(cell(at)')

"""
`set_defm!(at::AbstractAtoms, F::AbstractMatrix; updatepositions=false) -> at`

set the deformation matrix
"""
function set_defm!(at::AbstractAtoms, F::AbstractMatrix; updatepositions=false)
   if updatepositions
      A = JMatF(F * inv(defm(at)))
      X = [A * x for x in positions(at)]
      set_positions!(at, X)
   end
   set_cell!(at, Matrix(F'))
end

"set periodic boundary conditions"
@protofun set_pbc!(::AbstractAtoms, ::Union{AbstractVector, Tuple})
set_pbc!(at::AbstractAtoms, p::Bool) = set_pbc!(at, (p,p,p))

"get array (or tuple) determining which directions are periodic"
@protofun pbc(::AbstractAtoms)

"alias for `pbc`"
get_pbc = pbc

"determines whether a cubic cell is used (i.e. cell is a diagonal matrix)"
is_cubic(a::AbstractAtoms) = isdiag(cell(a))

"""
`set_data!(at, name, value)`:
associate some data with `at`; to be stored in a Dict within `at`

if `name::Union{Symbol, AbstractString}`, then `setindex!` is an alias
for `set_data!`.
"""
@protofun set_data!(a::AbstractAtoms, name::Any, value::Any)
Base.setindex!(at::AbstractAtoms, value,
               name::Union{Symbol, AbstractString}) = set_data!(at, name, value)

"""
`get_data(a, name)`:
obtain some data stored with `set_data!`

if `name::Union{Symbol, AbstractString}`, then `getindex` is an alias
for `get_data`.
"""
@protofun get_data(a::AbstractAtoms, name::Any)

Base.getindex(at::AbstractAtoms,
               name::Union{Symbol, AbstractString}) = get_data(at, name)

"check whether some data with id `name` is already stored"
@protofun has_data(a::AbstractAtoms, name::Any)

"delete an atom"
@protofun deleteat!(a::AbstractAtoms, n::Integer)

"alias for `deleteat!`"
delete_atom! = deleteat!

"return an attached calculator"
@protofun calculator(at::AbstractAtoms)

"attach a calculator"
@protofun set_calculator!(at::AbstractAtoms, calc::AbstractCalculator)

"return attached constraint"
@protofun constraint(at::AbstractAtoms)

"attach a constraint"
@protofun set_constraint!(at::AbstractAtoms, cons::AbstractConstraint)

#######################################################################
#  NEIGHBOURLIST
#######################################################################

"""
`neighbourlist(a::AbstractAtoms, rcut)`

construct a suitable neighbourlist for this atoms object. `rcut` should
be allowed to be either a scalar or a vector.
"""
@protofun neighbourlist(a::AbstractAtoms, rcut::AbstractFloat)

# """
# `sites(::AbstractNeighbourList)` or `sites(::AbstractAtoms, rcut)`
#
# Returns an iterator over atomic sites.
# ```julia
# for (idx, neigs, r, R, S) in sites(nlist)
#     # do something at this site
# end
# ```
# Here `idx` is the current center-atom, `neigs` a collection
# indexing the neighbouring atoms, `r` a vector if distances and
# `R` a vector of distance vectors. `S` stores information
# about which copy of the cell the neighbours belong to
#
# A quicker way, if `nlist` won't be reused is
# ```julia
# for (n, ...) in sites(at, rcut)
# ```
#
# It should be assumed that `neigs, r, R, S` are views and therefore
# should not be modified!
# """
# @protofun sites(::AbstractNeighbourList)
#
# sites(at::AbstractAtoms, rcut::Float64) = sites(neighbourlist(at, rcut))
#
#
# """
# `bonds` : iterator over (pair-) bonds
#
# E.g., the a pair potential can be implemented as follows:
# ```julia
# ϕ(r) = r^(-12) - 2.0 * r^(-6)
# dϕ(r) = -12 * (r^(-13) - r^(-7))
# function lj(at::AbstractAtoms)
#    E = 0.0; dE = zeros(3, length(at)) |> vecs
#    for (i, j, r, R, S) in bonds(at, 4.1)
#       E += ϕ(r)
#       dE[j] += (dϕ(r)/r) * R
#       dE[i] -= (dϕ(r)/r) * R
#    end
#    return E, dE
# end
# ```
# """
# @protofun bonds(::AbstractNeighbourList)
#
# bonds(at::AbstractAtoms, cutoff::Float64) = bonds(neighbourlist(at, cutoff))

#######################################################################
#     CALCULATOR
#######################################################################

type  NullCalculator <: AbstractCalculator end

"Returns the cut-off radius of the potential."
@protofun cutoff(::AbstractCalculator)

cutoff(at::AbstractAtoms) = cutoff(calculator(at))

"""
`energy`: can be called in various ways
*  `energy(calc, at)`: base definition; this is normally the only method that
     needs to be overloaded when a new calculator is implemented.
* `energy(at) = energy(calculator(at), at)`: if a calculator is attached to `at`
* `energy(calc, at, const, dof) = energy(calc, dof2at!(at,const,dof))`
* `energy(calc, at, dof) = energy(calc, at, constraint(dof), dof)`
* `energy(at, dof) = energy`

Return the total potential energy of a configuration of atoms `a`, using the calculator
`c`.
"""
@protofun energy(c::AbstractCalculator, a::AbstractAtoms)
energy(at::AbstractAtoms) = energy(calculator(at), at)

"`potential_energy` : alias for `energy`"
potential_energy = energy

@protofun site_energies(c::AbstractCalculator, a::AbstractAtoms)
site_energies(a::AbstractAtoms) = site_energies(calculator(a), a)

# """
# energy difference between two configurations; default is to just compute the
# two energies, but this allows implementation of numerically robust
# differences,w hich can be important for very large problems.
# """
# energy_difference(c::AbstractCalculator, a::AbstractAtoms, aref::AbstractAtoms) =
#    energy(c, a) - energy(c, aref)

"""
forces in `Vector{JVecF}`  (negative gradient w.r.t. atom positions only)
"""
@protofun forces(c::AbstractCalculator, a::AbstractAtoms)
forces(at::AbstractAtoms) = forces(calculator(at), at)

"""
hessian with respect to all atom positions
"""
@protofun hessian_pos(calc::AbstractCalculator, at::AbstractAtoms)
hessian_pos(at::AbstractAtoms) = hessian_pos(calculator(at), at)

"""
* `virial(c::AbstractCalculator, a::AbstractAtoms) -> JMatF`
* `virial(a::AbstractAtoms) -> JMatF`

returns virial, (- ∂E / ∂F) where `F = defm(a)`
"""
@protofun virial(c::AbstractCalculator, a::AbstractAtoms)

virial(a::AbstractAtoms) = virial(calculator(a), a)

"""
* `stress(c::AbstractCalculator, a::AbstractAtoms) -> JMatF`
* `stress(a::AbstractAtoms) -> JMatF`

stress = - virial / volume; this function should *not* be overloaded;
instead overload virial
"""
stress(c::AbstractCalculator, a::AbstractAtoms) = - virial(c, a) / det(defm(a))

stress(at::AbstractAtoms) = stress(calculator(at), at)

# remove these for now; not clear they are useful.
# @protofun site_virials(c::AbstractCalculator, a::AbstractAtoms)
# site_virials(a::AbstractAtoms) = site_virials(calculator(a), a)

#######################################################################
#  Constraints and DoF Handling
#######################################################################

type NullConstraint <: AbstractConstraint end

"""
`position_dofs(at::AbstractAtoms, cons::AbstractConstraint) -> Dofs`
`position_dofs(at::AbstractAtoms) -> Dofs`
`dofs(at::AbstractAtoms, cons::AbstractConstraint) -> Dofs`

Take an atoms object `at` and return a Dof-vector that fully describes the
state given the constraint `cons`
"""
@protofun position_dofs(at::AbstractAtoms, cons::AbstractConstraint)
position_dofs(at::AbstractAtoms) = position_dofs(at, constraint(at))

"""
`dofs` is an alias for `positions_dofs`
"""
dofs(args...) = position_dofs(args...)

"""
* `set_position_dofs!(at::AbstractAtoms, cons::AbstractConstraint, x::Dofs) -> at`
* `set_position_dofs!(at::AbstractAtoms, x::Dofs) -> at`
* `set_dofs!(at::AbstractAtoms, cons::AbstractConstraint, x::Dofs) -> at`
* `set_dofs!(at::AbstractAtoms, cons::AbstractConstraint) -> at`

change configuration stored in `at` according to `cons` and `x`.
"""
@protofun set_position_dofs!(at::AbstractAtoms, cons::AbstractConstraint, x::Dofs)
set_position_dofs!(at::AbstractAtoms, x::Dofs) = set_position_dofs!(at, constraint(at), x)
"""
`set_dofs!` is an alias for `set_position_dofs!`
"""
set_dofs!(args...) = set_position_dofs!(args...)

"""
`momentum_dofs(at::AbstractAtoms, cons::AbstractConstraint) -> Dofs`
`momentum_dofs(at::AbstractAtoms) -> Dofs`

Take an atoms object `at` and return a Dof-vector that fully describes the
momenta (given the constraint `cons`)
"""
@protofun momentum_dofs(at::AbstractAtoms, cons::AbstractConstraint)
momentum_dofs(at::AbstractAtoms) = momentum_dofs(at, constraint(at))

"""
* `set_momentum_dofs!(at::AbstractAtoms, cons::AbstractConstraint, p::Dofs) -> at`
* `set_momentum_dofs!(at::AbstractAtoms, p::Dofs) -> at`

change configuration stored in `at` according to `cons` and `p`.
"""
@protofun set_momentum_dofs!(at::AbstractAtoms, cons::AbstractConstraint, q::Dofs)
set_momentum_dofs!(at::AbstractAtoms, q::Dofs) = set_momentum_dofs!(at, constraint(at), q)

"""
* `project!(at::AbstractAtoms, cons::AbstractConstraint) -> at`
* `project!(at::AbstractAtoms) -> at`

project the `at` onto the constraint manifold.
"""
@protofun project!(at::AbstractAtoms, cons::AbstractConstraint)

project!(at::AbstractAtoms) = project!(at, constraint(at))

# converting calculator functionality

"""
* `gradient(at, cons::AbstractConstraint) -> Float64`
* `gradient(at, x::Dofs) -> Float64`

`potential_energy`; with potentially added cell terms, e.g., if there is
applied stress or applied pressure.
"""
@protofun energy(at::AbstractAtoms, cons::AbstractConstraint)

energy(at::AbstractAtoms, x::Dofs) = energy(set_dofs!(at, x), constraint(at))

"""
* `gradient(at) -> Dofs`
* `gradient(at, cons::AbstractConstraint) -> Dofs`
* `gradient(at, x::Dofs) -> Dofs`

gradient of `potential_energy` in dof-format;
depending on the constraint this could include e.g. a gradient w.r.t. cell
shape
"""
@protofun gradient(at::AbstractAtoms, cons::AbstractConstraint)

gradient(at::AbstractAtoms, x::Dofs) = gradient(set_dofs!(at, x), constraint(at))

gradient(at::AbstractAtoms) = gradient(at, constraint(at))

"""
`hessian`: compute hessian of total energy with respect to DOFs!

TODO: write docs
"""
@protofun hessian(at::AbstractAtoms, cons::AbstractConstraint)

hessian(at::AbstractAtoms, x::Dofs) = hessian(set_dofs!(at, x), constraint(at))
hessian(at::AbstractAtoms) = hessian(at, constraint(at))

#######################################################################
#                     PRECONDITIONER
#######################################################################

"""
`update!(precond::Preconditioner, at::AbstractAtoms)`

Update the preconditioner with the new geometry information.
"""
@protofun update!(precond::Preconditioner, at::AbstractAtoms)
update!(precond::Preconditioner, at::AbstractAtoms, x::Dofs) =
            update!(precond, set_dofs!(at, x))
# TODO: this is a bit of a problem with the niceabstractframework we
#       are constructing here; it can easily happen now that we
#       update positions multiple times

"""
Identity preconditioner, i.e., no preconditioner.
"""
type Identity <: Preconditioner
end

A_ldiv_B!(out::Dofs, P::Identity, x::Dofs) = copy!(out, x)
A_mul_B!(out::Dofs, P::Identity, f::Dofs) = copy!(out, f)
update!(P::Identity, at::AbstractAtoms) = P

"construct a preconditioner suitable for this atoms object"
preconditioner(at::AbstractAtoms) = preconditioner(at, calculator(at), constraint(at))
preconditioner(at::AbstractAtoms, calc::AbstractCalculator, con::AbstractConstraint) =
         Identity()


text:
# Monkey patch `Base.package_slug`
#
# See also:
# * Suggestion: Use different precompilation cache path for different
#   system image -- https://github.com/JuliaLang/julia/pull/29914
#
Base.eval(Base, quote
    function package_slug(uuid::UUID, p::Int=5)
        crc = _crc32c(uuid)
        crc = _crc32c(unsafe_string(JLOptions().image_file), crc)
        crc = _crc32c(get(ENV, "JLM_PRECOMPILE_KEY", ""), crc)
        return slug(crc, p)
    end
end)


text:
export ptiread

"""
    ptiread

Read signals from *.pti files measured using B&K system.
### Fields:
* `filename`    - file name of the *.pti file.
### Available Constructors:
* `Signal, Freq, Date, Time = ptiread('recording-1.1.pti')`
"""
function ptiread(filename::AbstractString)

    # extract header lines
    header = Array{Union{Nothing, String}}(nothing, 297)
    open(filename) do io
        for i=1:297
            header[i] = readline(io)
        end
    end

    # determine start header line
    st_ref = 2
    for i=1:5
        if cmp("[SETUP START]",header[i]) == 0 # position of this line
            st_ref = i
            #println(st_ref)
        end
    end

    # extract general info
    RECInfoSectionSize = parse(Int64, match(r"(?<=RECInfoSectionSize=)(\d+)", header[st_ref+2]).match)
    RECInfoSectionPos = parse(Int64, match(r"(?<=RECInfoSectionPos=)(\d+)", header[st_ref+3]).match)
    SampleFrequency = parse(Int64, match(r"(?<=SampleFrequency=)(\d+)", header[st_ref+4]).match)
    NoChannels = parse(Int8, match(r"(?<=NoChannels=)(\d+)", header[st_ref+5]).match)
    OffsetStopSample = parse(Int64, match(r"(?<=OffsetStopSample=)(\d+)", header[st_ref+11]).match)
    Date = header[st_ref+12][6:end]
    Time = header[st_ref+13][6:end]

    # extract correction for each channels
    CorrectionFactor = Array{Float64,1}(undef, NoChannels)

    for channel=1:NoChannels
        k = 15 + (channel-1)*10 + 4 + (st_ref-2)

        cfactor = match(r"(?<=CorrectionFactor=)(-?\ *[0-9]+\.?[0-9]*(?:[Ee]\ *-?\ *[0-9]+)?)",
            header[k])
        CorrectionFactor[channel] = parse(Float64, cfactor.match)
    end

    # Read binary data
    data = open(filename) do io
        seek(io, RECInfoSectionPos + RECInfoSectionSize + 20)
        dsize = read(io,Int16)
        cols = Int(OffsetStopSample/(dsize-4)*NoChannels)
        data = Array{Int32}(undef, (dsize, cols))

        seek(io, RECInfoSectionPos + RECInfoSectionSize + 20)
        read!(io, data)
    end

    #reshape channels
    Signal = Array{Float64}(undef, (OffsetStopSample, NoChannels))
    factor = CorrectionFactor/(2^16)

    for i=1:NoChannels
        Signal[:,i] = reshape(data[5:end,i:NoChannels:end],:,1) * factor[i]
    end

    return Signal, SampleFrequency, Date, Time

end

#filename = "R:\\CMPH-Windfarm Field Study\\Hallett\\set1\\Recording-3.8.pti"
#@time ptiread(filename)


text:
function get_mesh_1d(i)

    coordinates_min = (-1,)
    coordinates_max = (1,)

    mesh1_1d = TreeMesh(coordinates_min, coordinates_max,
                    initial_refinement_level=4,
                    n_cells_max=10_000
                    )

    coordinates_min = (-1,)
    coordinates_max = (1,)
    # refinement_patches = (
    #  (type="box", coordinates_min=(-1.0), coordinates_max=(1.0)),
    # )
    mesh2_1d = TreeMesh(coordinates_min, coordinates_max,
                    initial_refinement_level=5,
                    # refinement_patches=refinement_patches,
                    n_cells_max=10_000
                    )

    coordinates_min = (0,)
    coordinates_max = (2,)
    refinement_patches = (
    (type = "box", coordinates_min = (0.5), coordinates_max = (1.5)),
    )
    mesh3_1d = TreeMesh(coordinates_min, coordinates_max,
                    initial_refinement_level=6,
                    refinement_patches=refinement_patches,
                    n_cells_max=10_000
                    )

    coordinates_min = (-1,)
    coordinates_max = (1,)
    refinement_patches = (
    (type = "box", coordinates_min = (-0.5), coordinates_max = (0.5)),
    (type = "box", coordinates_min = (-0.2), coordinates_max = (0.2)),
    )
    mesh4_1d = TreeMesh(coordinates_min, coordinates_max,
                    initial_refinement_level=5,
                    refinement_patches=refinement_patches,
                    n_cells_max=10_000
                    )
    if i == 1
        mesh = [mesh1_1d]
    elseif i == 2
        mesh = [mesh2_1d]
    elseif i == 3
        mesh = [mesh3_1d]
    elseif i == 4
        mesh = [mesh4_1d]
    elseif i == 5
        return [mesh1_1d,mesh2_1d,mesh3_1d,mesh4_1d]
    end
    return mesh
end


text:
using DiffEqBiological, Plots
gr(); default(fmt = :png);

bistable_switch = @reaction_network begin
    d,    (X,Y) → ∅
    hillR(Y,v1,K1,n1), ∅ → X
    hillR(X,v2,K2,n2), ∅ → Y
end d v1 K1 n1 v2 K2 n2
d = 0.01;
v1 = 1.5; K1 = 30; n1 = 3;
v2 = 1.; K2 = 30; n2 = 3;
bistable_switch_p = [d, v1 ,K1, n1, v2, K2, n2];

ss = steady_states(bistable_switch, bistable_switch_p)

stability(ss,bistable_switch, bistable_switch_p)

rn1 = @reaction_network begin
    p, ∅ → X
    hill(X,v,K,n), X → ∅
end p v K n
p1 = [1.,2.5,1.5,1.5]
steady_states(rn1,p1)

rn2 = @reaction_network begin
    p, ∅ → X
    log(X), X → ∅
end p
p2 = [1.]
steady_states(rn2,p2)

bif = bifurcations(bistable_switch, bistable_switch_p, :v1, (.1,5.))
plot(bif,ylabel="[X]",label="")
plot!([[],[]],color=[:blue :red],label = ["Stable" "Unstable"])

plot(bif,2,ylabel="[Y]")
plot!([[],[]],color=[:blue :red],label = ["Stable" "Unstable"])

bif = bifurcations(bistable_switch, bistable_switch_p,:v1,(.1,10.))
plot(bif,linewidth=1.,title="A bifurcation diagram",ylabel="Steady State concentration")
plot!([[],[]],color=[:blue :red],label = ["Stable" "Unstable"])

bif = bifurcation_grid(bistable_switch, bistable_switch_p,:n1,1.:5.)
plot(bif)
scatter!([[],[]],color=[:blue :red],label = ["Stable" "Unstable"])

bif = bifurcation_grid_diagram(bistable_switch, bistable_switch_p,:n1,0.:4.,:v1,(.1,5.))
plot(bif)
plot!([[],[]],color=[:blue :red],label = ["Stable" "Unstable"])

bif = bifurcation_grid_2d(bistable_switch, bistable_switch_p,:n1,1.:3.,:n2,1.:10.)
plot(bif)
scatter!([[],[]],color=[:blue :red],label = ["Stable" "Unstable"])

brusselator = @reaction_network begin
    A, ∅ → X
    1, 2X + Y → 3X
    B, X → Y
    1, X → ∅
end A B;
A = 0.5; B = 4.;
brusselator_p = [A, B];

bif = bifurcations(brusselator,brusselator_p,:B,(0.1,2.5))
plot(bif,2)
plot!([[],[],[],[]],color=[:blue :cyan :orange :red],label = ["Stable Real" "Stable Complex" "Unstable Complex" "Unstable Real"])

bif = bifurcation_grid_diagram(brusselator,brusselator_p,:B,0.5:0.02:5.0,:A,(0.2,5.0))
plot(bif)
plot!([[],[],[],[]],color=[:blue :cyan :orange :red],label = ["Stable Real" "Stable Complex" "Unstable Complex" "Unstable Real"])

using DiffEqTutorials
DiffEqTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file], remove_homedir=true)


text:
function smolyak_derivative_finite_difference(weights::Array{T,1},node::Array{T,1},multi_index::Array{S,2}) where {T<:AbstractFloat,S<:Integer}

  m = length(node)
  e  = eps(T)^(1/3)*maximum(abs,[node;one(T)])
  dh = Matrix{T}(I,m,m)*e

  evaluated_derivative = zeros(1,m)

  for i = 1:m
    f1 = smolyak_evaluate(weights,node+2*dh[:,i],multi_index)
    f2 = smolyak_evaluate(weights,node+dh[:,i],multi_index)
    f3 = smolyak_evaluate(weights,node-dh[:,i],multi_index)
    f4 = smolyak_evaluate(weights,node-2*dh[:,i],multi_index)
    evaluated_derivative[i] = (-f1+8*f2-8*f3+f4)/(12*e)
  end

  return evaluated_derivative

end

function smolyak_derivative_finite_difference(weights::Array{T,1},node::Array{T,1},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,S<:Integer}

  if size(domain,2) != length(node)
    error("domain is inconsistent with the number of dimensions")
  end

  m = length(node)
  e  = eps(T)^(1/3)*maximum(abs,[node;one(T)])
  dh = eye(m)*e

  evaluated_derivative = zeros(1,m)

  for i = 1:m
    f1 = smolyak_evaluate(weights,node+2*dh[:,i],multi_index,domain)
    f2 = smolyak_evaluate(weights,node+dh[:,i],multi_index,domain)
    f3 = smolyak_evaluate(weights,node-dh[:,i],multi_index,domain)
    f4 = smolyak_evaluate(weights,node-2*dh[:,i],multi_index,domain)
    evaluated_derivative[i] = (2.0/(domain[1,i]-domain[2,i]))*(-f1+8*f2-8*f3+f4)/(12*e)
  end

  return evaluated_derivative

end


text:
function solve()
    n, a, m = [parse(Int, x) for x in split(readline())]
    d = Dict{Int, Int}(a=>1)
    v = [a]
    i = 2
    while i ≤ n
        a = a^2 % m
        if a ∈ keys(d)
            break
        end
        d[a] = i
        push!(v, a)
        i += 1
    end
    if i > n
        sum(v)
    else
        cyclestart = d[a]
        cyclelen = i - cyclestart
        ret = sum(v[1:cyclestart-1])
        n = n - cyclestart
        ret += sum(v[cyclestart:end])* (n ÷ cyclelen)
        ret += sum(v[cyclestart:cyclestart+(n%cyclelen)])
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text: println("I am being built...")

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
using StatsFuns: logistic
using Plots

βs = [0, 0.2, 0.5, 0.8, 1.0, 1.2, 1.5, 2, 3, 5, 7, 8, 10, 20]

plot([x -> logistic(β * x) for β in βs],
        -2, 2,
        alpha = 0.8,
        legend = :none,
        ylims = [0, 1])

max_stars_repo_name: kogpsy/vpompd

max_issues_repo_name: kogpsy/vpompd

max_forks_repo_name: kogpsy/vpompd


text:
function get_sim_data_direct(mechanism::Mechanism,env::Environment,Δt::Real;
                             relax_comp=false)
    num_contacts = length(env.contacts)

    vs = VariableSelector()
    add_var!(vs, :qnext, num_positions(mechanism))
    add_var!(vs, :vnext, num_velocities(mechanism))
    if relax_comp
        add_var!(vs, :slack, 1)
    end

    cs = ConstraintSelector()
    add_eq!(cs, :kin, num_positions(mechanism))
    add_eq!(cs, :dyn, num_velocities(mechanism))

    state_cache = [StateCache(mechanism) for n = 1:2]
    envj_cache = [EnvironmentJacobianCache(env) for n = 1:2]

    generate_solver_fn = :generate_solver_fn_sim_direct
    extract_sol = :extract_sol_sim_direct

    normal_vs = VariableSelector()
    for i = 1:length(env.contacts)
        add_var!(normal_vs, Symbol("c_n", i), 1)
    end

    normal_cs = ConstraintSelector()
    for i = 1:length(env.contacts)
        add_ineq!(normal_cs, Symbol("c_n_pos", i), 1)
        add_ineq!(normal_cs, Symbol("ϕ", i), 1)
    end

    normal_options = Dict{String, Any}()
    normal_options["num_fosteps"] = 1
    normal_options["num_sosteps"] = 10
    normal_options["c"] = 10
    normal_options["c_fos"] = 10
    normal_options["c_sos"] = 1

    fric_vs = VariableSelector()
    for i = 1:length(env.contacts)
        β_dim = size(env.contacts[i].obstacle.basis,2)
        add_var!(fric_vs, Symbol("β", i), β_dim)
    end

    fric_cs = ConstraintSelector()
    for i = 1:length(env.contacts)
        β_dim = size(env.contacts[i].obstacle.basis,2)
        add_ineq!(fric_cs, Symbol("β_pos", i), β_dim)
        add_ineq!(fric_cs, Symbol("fric_cone", i), 1)
    end

    fric_options = Dict{String, Any}()
    fric_options["num_fosteps"] = 1
    fric_options["num_sosteps"] = 10
    fric_options["c"] = 10
    fric_options["c_fos"] = 10
    fric_options["c_sos"] = 1

    fric_P = sparse(ones(4,4))
    fric_q = ones(4)
    fric_A = sparse(ones(5,4))
    fric_l = ones(5)
    fric_u = ones(5)
    fric_model = OSQP.Model()
    OSQP.setup!(fric_model, P=fric_P, q=fric_q, A=fric_A, l=fric_l, u=fric_u, verbose=0)

    normal_P = sparse(ones(1,1))
    normal_q = ones(1)
    normal_A = sparse(ones(2,1))
    normal_l = ones(2)
    normal_u = ones(2)
    normal_model = OSQP.Model()
    OSQP.setup!(normal_model, P=normal_P, q=normal_q, A=normal_A, l=normal_l, u=normal_u, verbose=0)

    SimData(mechanism,env,
            state_cache,envj_cache,
            Δt,vs,cs,generate_solver_fn,extract_sol,
            [],[],[],[normal_vs],[normal_cs],[normal_options],[fric_vs],[fric_cs],[fric_options],1,[],[],[fric_model],[normal_model])
end

function extract_sol_sim_direct(sim_data::SimData, results::AbstractArray{T,2}) where T
    vs = sim_data.vs
    relax_comp = haskey(vs.vars, :slack)
    env = sim_data.env
    N = size(results,2)

    qtraj = Array{Array{Float64,1},1}(undef, 0)
    vtraj = Array{Array{Float64,1},1}(undef, 0)
    utraj = Array{Array{Float64,1},1}(undef, 0)
    contact_traj = Array{Array{Float64,1},1}(undef, 0)
    slack_traj = Array{Array{Float64,1},1}(undef, 0)
    x = MechanismState(sim_data.mechanism)
    for n = 1:N
        set_configuration!(x, vs(results[:,n], Symbol("qnext")))
        normalize_configuration!(x)
        qnext = configuration(x)
        push!(qtraj, qnext)
        push!(vtraj, vs(results[:,n], Symbol("vnext")))
        if relax_comp
            push!(slack_traj, vs(results[:,n], Symbol("slack")))
        end
        # TODO compute contact force for comparison purposes
    end

    # some other usefull vectors
    ttraj = [(i-1)*sim_data.Δt for i = 1:N]
    qv_mat = vcat(hcat(qtraj...),hcat(vtraj...))

    qtraj, vtraj, utraj, contact_traj, slack_traj, ttraj, qv_mat, results
end

function generate_solver_fn_sim_direct(sim_data,q0,v0,u0)
    x0 = sim_data.state_cache[1][Float64]
    envj = sim_data.envj_cache[1][Float64]
    Δt = sim_data.Δt
    vs = sim_data.vs
    cs = sim_data.cs

    relax_comp = haskey(vs.vars, :slack)
    num_contacts = length(sim_data.env.contacts)
    num_vel = num_velocities(sim_data.mechanism)

    set_configuration!(x0, q0)
    set_velocity!(x0, v0)
    setdirty!(x0)
    H = mass_matrix(x0)
    Hi = inv(H)

    contact_jacobian!(envj, x0)
    dyn_bias0 = dynamics_bias(x0) # TODO preallocate

    function eval_obj(x::AbstractArray{T}) where T
        f = 0.

        if relax_comp
            slack = vs(x, :slack)
            f += .5 * slack' * slack
        end

        f
    end

    function eval_cons(x::AbstractArray{T}) where T
        xn = sim_data.state_cache[2][T]

        normal_bias = Vector{T}(undef, num_vel)
        contact_bias = Vector{T}(undef, num_vel)
        g = Vector{T}(undef, cs.num_eqs + cs.num_ineqs) # TODO preallocate

        qnext = vs(x, :qnext)
        vnext = vs(x, :vnext)
        if relax_comp
            slack = vs(x, :slack)
        end

        set_configuration!(xn, qnext)
        set_velocity!(xn, vnext)
        setdirty!(xn)

        normalize_configuration!(xn)

        if (num_contacts > 0)
            # compute normal forces
            # x_normal = contact_normal_τ_direct!(normal_bias, sim_data, sim_data.Δt, Hi, envj, dyn_bias0, u0, v0, x, 1)
            x_normal = contact_normal_τ_direct_osqp!(normal_bias, sim_data, sim_data.Δt, Hi, envj, dyn_bias0, u0, v0, x, 1)

            # compute friction forces
            # contact_friction_τ_direct!(contact_bias, sim_data, sim_data.Δt, Hi, envj, dyn_bias0, u0, v0, x, x_normal, 1)
            contact_friction_τ_direct_osqp!(contact_bias, sim_data, sim_data.Δt, Hi, envj, dyn_bias0, u0, v0, x, x_normal, 1)
        end
        config_derivative = configuration_derivative(xn) # TODO preallocate
        dyn_bias = dynamics_bias(xn) # TODO preallocate

        g[cs(:kin)] .= qnext .- q0 .- Δt .* config_derivative
        g[cs(:dyn)] .= H * (vnext - v0) .- Δt .* (u0 .- dyn_bias .- contact_bias)

        g
    end

    generate_autodiff_solver_fn(eval_obj,eval_cons,cs.eqs,cs.ineqs,vs.num_vars)
end


text:
struct MSARIMAParams{T <: Integer} <: AbstractParams
    c::Bool
    p::Vector{Vector{T}}
    d::Vector{T}
    q::Vector{Vector{T}}
    s::Vector{T}
    function MSARIMAParams(c::Bool, p::Vector{Vector{T}}, d::Vector{T}, q::Vector{Vector{T}}, s::Vector{T}) where T <: Integer
        length(p) == length(q) == length(d) == length(s) || throw(DimensionMismatch("c,p,q,d and s should have the same length"))
        return new{T}(c, p, d, q, s)
    end
end

MSARIMAParams(c::Bool, p::Vector{T},d::Vector{T},q::Vector{T},s::Vector{T}) where T <: Integer = MSARIMAParams(c, map(i -> collect(1:i), p), d, map(i -> collect(1:i), q), s)
MSARIMAParams(p::Vector{T},d::Vector{T},q::Vector{T},s::Vector{T}) where T <: Integer = MSARIMAParams(true, p, d, q, s)

function toarma(params::MSARIMAParams, armas::Vector{ARMAModel})
    M = length(armas)
    ϕp = Polynomial(1)
    θp = Polynomial(1)
    for i in 1:M
        ϕp *= Polynomial([1;.-armas[i].ϕ]) * Polynomial(1 - variable()^params.s[i])^params.d[i]
        θp *= Polynomial([1;.-armas[i].θ])
    end
    p0 = length(ϕp) - 1
    q0 = length(θp) - 1
    return ARMAModel{p0,q0}(armas[1].μ, .-coeffs(ϕp)[2:end], .-coeffs(θp)[2:end], armas[1].σ2)
end

function fit(params::MSARIMAParams, z::AbstractVector; n::Integer=10)
    M = length(params.p)
    N = length(z)
    N > maximum(params.d) * (1 + maximum(params.s)) || throw(DomainError("ARIMA requires at least d+D*s data points"))
    armas = ARMAModel[]
    ∇z = z .- mean(z)
    for i in 1:M
        ∇z = difference(∇z, d=params.d[i], s=params.s[i])
    end
    for i in 1:M
        arma = fit(ARMAParams(params.c, params.p[i] .* params.s[i], params.q[i] .* params.s[i]), ∇z, n=n)
        push!(armas, arma)
    end
    return toarma(params, armas)
end


text:
struct DomTreeNode
    level::Int
    children::Vector{Int}
end
DomTreeNode() = DomTreeNode(1, Vector{Int}())

struct DomTree
    idoms::Vector{Int}
    nodes::Vector{DomTreeNode}
end

"""
    Checks if bb1 dominates bb2
"""
function dominates(domtree, bb1, bb2)
    bb1 == bb2 && return true
    target_level = domtree.nodes[bb1].level
    source_level = domtree.nodes[bb2].level
    source_level < target_level && return false
    for _ in (source_level-1):-1:target_level
        bb2 = domtree.idoms[bb2]
    end
    return bb1 == bb2
end

function update_level!(domtree, node, level)
    domtree[node] = DomTreeNode(level, domtree[node].children)
    foreach(domtree[node].children) do child
        update_level!(domtree, child, level+1)
    end
end

struct DominatedBlocks
    domtree::DomTree
    worklist::Vector{Int}
end

function dominated(domtree::DomTree, root::Int)
    doms = DominatedBlocks(domtree, Vector{Int}())
    push!(doms.worklist, root)
    doms
end

function start(doms::DominatedBlocks)
    nothing
end

function next(doms::DominatedBlocks, state::Nothing)
    bb = pop!(doms.worklist)
    for dominated in doms.domtree.nodes[bb].children
        push!(doms.worklist, dominated)
    end
    (bb, nothing)
end

function done(doms::DominatedBlocks, state::Nothing)
    isempty(doms.worklist)
end

# Construct Dom Tree
# Simple algorithm - TODO: Switch to the fast version (e.g. https://tanujkhattar.wordpress.com/2016/01/11/dominator-tree-of-a-directed-graph/)
function construct_domtree(cfg)
    dominators = IdSet{Int}[n == 1 ? IdSet{Int}(n) : IdSet{Int}(1:length(cfg.blocks)) for n = 1:length(cfg.blocks)]
    changed = true
    while changed
        changed = false
        for n = 2:length(cfg.blocks)
            isempty(cfg.blocks[n].preds) && continue
            firstp, rest = Iterators.peel(cfg.blocks[n].preds)
            new_doms = copy(dominators[firstp])
            for p in rest
                intersect!(new_doms, dominators[p])
            end
            push!(new_doms, n)
            changed |= (new_doms != dominators[n])
            dominators[n] = new_doms
        end
    end
    # Compute idoms
    idoms = fill(0, length(cfg.blocks))
    for i = 2:length(cfg.blocks)
        for dom in dominators[i]
            i == dom && continue
            any(p->p !== i && p !== dom && dom in dominators[p], dominators[i]) && continue
            idoms[i] = dom
        end
    end
    # Compute children
    domtree = DomTreeNode[DomTreeNode() for _ = 1:length(cfg.blocks)]
    for (idx, idom) in Iterators.enumerate(idoms)
        (idx == 1 || idom == 0) && continue
        push!(domtree[idom].children, idx)
    end
    # Recursively set level
    update_level!(domtree, 1, 1)
    DomTree(idoms, domtree)
end

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
using Gurobi, Base.Test, MathOptInterface, MathOptInterface.Test, MathOptInterfaceGurobi

const MOIT = MathOptInterface.Test
const MOIGRB = MathOptInterfaceGurobi

@testset "MathOptInterfaceGurobi" begin
    @testset "Linear tests" begin
        linconfig = linconfig = MOIT.TestConfig()
        solver = GurobiOptimizer(OutputFlag=0)
        MOIT.contlineartest(solver, linconfig, ["linear10","linear12","linear8a","linear8b","linear8c"])

        solver_nopresolve = GurobiOptimizer(OutputFlag=0, InfUnbdInfo=1)
        MOIT.contlineartest(solver_nopresolve, linconfig, ["linear10","linear12","linear8a"])

        linconfig_nocertificate = MOIT.TestConfig(infeas_certificates=false)
        MOIT.linear12test(solver, linconfig_nocertificate)
        MOIT.linear8atest(solver, linconfig_nocertificate)

        # 10 is ranged
    end

    @testset "Quadratic tests" begin
        quadconfig = MOIT.TestConfig(atol=1e-4, rtol=1e-4, duals=false, query=false)
        solver = GurobiOptimizer(OutputFlag=0)
        MOIT.contquadratictest(solver, quadconfig)
    end

    @testset "Linear Conic tests" begin
        linconfig = MOIT.TestConfig()
        solver = GurobiOptimizer(OutputFlag=0)
        MOIT.lintest(solver, linconfig, ["lin3","lin4"])

        solver_nopresolve = GurobiOptimizer(OutputFlag=0, InfUnbdInfo=1)
        MOIT.lintest(solver_nopresolve, linconfig)
    end

    @testset "Integer Linear tests" begin
        intconfig = MOIT.TestConfig()
        solver = GurobiOptimizer(OutputFlag=0)
        MOIT.intlineartest(solver, intconfig, ["int3"])

        # 3 is ranged
    end
    @testset "ModelLike tests" begin
        intconfig = MOIT.TestConfig()
        solver = GurobiOptimizer()
        MOIT.validtest(solver)
        MOIT.emptytest(solver)
        solver2 = GurobiOptimizer()
        MOIT.copytest(solver,solver2)
    end
end
;