text:
module Models

import Base: +, -, *, isapprox, convert, copy, show
using LinearAlgebra
using FlightMechanics

export
# Attitude
Attitude, get_euler_angles, get_quaternions,
# Position
Position, LLHPosition, ECEFPosition, EarthPosition,
get_position, get_llh, get_xyz_earth, get_xyz_ecef, get_height,
# State
State,
get_attitude,
get_body_velocity, get_horizon_velocity, get_flight_path_angle,
get_body_ang_velocity, get_turn_rate,
get_body_accel, get_horizon_accel,
get_body_ang_accel,
# Mass
RigidSolid, PointMass, get_mass, get_cg, get_inertia,
# Pont Forces and moments
PointForcesMoments, translate_forces_moments, rotate,
get_point, get_forces, get_moments,
# Atmosphere
AtmosphereISA,AtmosphereF16,
get_temperature, get_pressure, get_density, get_sound_velocity,
calculate_atmosphere,
# Wind
Wind, ConstantWind,
get_wind_NED, get_wind_body, get_wind_direction, get_wind_intensity,
get_wind_vertical, calculate_wind,
# Gravity
EarthConstantGravity, get_gravity_horizon, get_gravity_accel, get_gravity_body,
calculate_gravity,
# Environment
Environment, calculate_environment, get_atmos, get_wind, get_gravity,
# Aerostate
AeroState, get_alpha, get_beta, get_aero_angles, get_alpha_dot, get_tas,
get_eas, get_cas, get_ias, get_aero_speeds, get_qinf, get_mach,
generate_state_aerostate, check_state_aerostate_env_coherence,
# Controls
Control, RangeControl, DiscreteControl,
get_value, get_value_range, get_value_choices, set_value!,
# Inputs
Input, ConstantInput, StepInput, DoubletInput, InverseDoubletInput, RampInput,
SinusoidalInput,
get_value,
# FCS
FCS,
set_stick_lon!, set_stick_lat!, set_pedals!, set_thtl!,
get_allow_out_of_range_inputs, get_throw_error_on_out_of_range_inputs,
set_allow_out_of_range_inputs!, set_throw_error_on_out_of_range_inputs!,
get_controls_ranges_trimmer,
Controls, StickPedalsLeverControls,
set_controls!,
get_controls_array, get_n_controls, get_stick_lon, get_stick_lat, get_pedals, get_thrust_lever,
ControlsStream, StickPedalsLeverStream,
get_controls,
# Engine
Engine, get_engine_position, get_engine_orientation, get_engine_gyro_effects,
# Propulsion
Propulsion, get_propulsion_position, get_fuel_mass_props, calculate_propulsion,
get_engines, get_gyro_effects,
# Aerodynamics
Aerodynamics, aerodynamics_from_wind_total, aerodynamics_from_wind_coeff,
aerodynamics_from_body_total, aerodynamics_from_body_coeff, get_pfm,
get_wind_pfm, get_wind_adim_pfm, get_body_pfm, get_body_adim_pfm,
calculate_aerodynamics,
# Aircraft
Aircraft, get_mass_props, get_pfm, get_aerodynamics, get_propulsion, get_fcs, get_name,
get_wing_area, get_wing_span, get_chord, get_arp, get_empty_mass_props,
get_payload_mass_props, calculate_aircraft,
# DynamicSystem
SixDOFEulerFixedMass, SixDOFQuaternionFixedMass, SixDOFECEFQuaternionFixedMass,
SixDOFAeroEulerFixedMass,
get_state_equation, get_state_equation_ode_wrapper, get_x, get_x_count, get_state

include("attitude.jl")
include("position.jl")
include("state.jl")
include("mass.jl")
include("point_forces_moments.jl")
include("atmosphere.jl")
include("wind.jl")
include("gravity.jl")
include("environment.jl")
include("aero_state.jl")
include("controls.jl")
include("inputs.jl")
include("fcs.jl")
include("propulsion.jl")
include("aerodynamics.jl")
include("aircraft.jl")
include("dynamic_systems.jl")
end


text:
@i @inline function ⊕(angle)(r!::T, x::Complex{T}) where T<:GVar
    r! += atan(x.im, x.re)
end

@i @inline function ⊕(abs2)(y!::T, a::Complex{T}) where T<:GVar
    y! += a.re^2
    y! += a.im^2
end

@i @inline function ⊕(abs)(y!::T, a::Complex{T}) where T<:GVar
    @routine @invcheckoff begin
        y2 ← zero(y!)
        y2 += abs2(a)
    end
    y! += sqrt(y2)
    ~@routine
end

Base.zero(x::Complex{T}) where T<:GVar = Complex(zero(T), zero(T))
Base.zero(::Type{Complex{T}}) where T<:GVar = Complex(zero(T), zero(T))
Base.one(x::Complex{T}) where T<:GVar = Complex(one(T), zero(T))
Base.one(::Type{Complex{T}}) where T<:GVar = Complex(one(T), zero(T))

ext: jl


text:
@testset "DeviceModel Tests" begin
    @test_throws ArgumentError DeviceModel(PSY.ThermalGen, PSI.ThermalStandardUnitCommitment)
    @test_throws ArgumentError DeviceModel(PSY.ThermalStandard, PSI.AbstractDeviceFormulation)
end

@testset "OperationsProblem Tests" begin
    for p in [true, false]
        t = OperationsProblem(TestOpProblem, CopperPlatePowerModel, c_sys5; use_parameters = p)
        @test PSI.model_has_parameters(t.canonical) == p
    end
end


text:
using MatterPower
using Roots
using Plots, LaTeXStrings

# %% Specify a redshift
redshift = 0

# %% Define a function to return a linear matter power spectrum (in units of Mpc^3/h^3)
# as a function of the comoving wavenumber, kovh, in units of h/Mpc.
# Here is an example using Einstein & Hu's analytical transfer function

# Cosmological parameters
As, ns, kpivot = 2.097e-9, 0.9652, 0.05
Ωm, ΩΛ, Ωb = 0.315, 0.685, 0.049
Ωk = 1 - Ωm - ΩΛ
h0 = 0.674
ωm, fb = Ωm * h0^2, Ωb / Ωm

# Tabulate linear growth factor as a function of scale factor, a
sol = setup_growth(Ωm, ΩΛ)
a = 1 / (1 + redshift)
D1 = sol(a)[1]

pk(kovh) =
   D1^2 *
   As *
   (kovh * h0 / kpivot)^(ns - 1) *
   (2 * kovh^2 * 2998^2 / 5 / Ωm)^2 *
   t_nowiggle(kovh * h0, ωm, fb)^2 *
   2 *
   π^2 / kovh^3

# %% Alternatively you may read in pre-computed data and define a spline function
# using Dierckx
# pk = Spline1D(tabulated_kovh, tabulated_power_spectrum)

# %% Compute the r.m.s. mass fluctuation with a top-hat radius Rh
Rh = 8 # Mpc/h
σ2 = sigma2(pk, Rh)
println("σ(R) = ", √σ2, " for R = ", Rh, " Mpc/h")

# %% Compute the non-linear mass, M*, defined by σ(R*,z) = 1.6865
# and M* = (4π/3)ρm R*^3 = (4π/3)ρc Ωm R*^3
# ρc is the present-day critical density of the Universe in units of h^2 M⊙/Mpc^3
ρc = 2.775e11
f(x) = √sigma2(pk, x) - 1.6865
Rstar = find_zero(f, (1e-2, 3), Bisection())
println("R* = ", Rstar, " h/Mpc")
Mstar = 4π / 3 * ρc * Ωm * Rstar^3
println("Non-linear mass M* = ", Mstar, " M⊙/h")

# %% Compute the non-linear power spectrum using halofit
# Get three parameters [kσ, neff, C] needed for the halofit transform
p = kσ, neff, C = setup_halofit(pk)

# Calculate Ωm(z) = Ωm(z=0)(1+z)^3 / E^2(z), where E^2(z) = H^2(z) / H0^2
z = redshift
E2 = Ωm * (1 + z)^3 + Ωk * (1 + z)^2 + ΩΛ
Ωmz = Ωm * (1 + z)^3 / E2

# Define a function to return a halofit non-linear power spectrum
pknl(kovh) = halofit(pk, p, Ωmz, kovh) # Mpc^3/h^3

# %% Plot results and save to "pk.pdf"
lnk = log(3e-4):0.3:log(30)
kovh = exp.(lnk)
p = plot(
   kovh,
   pknl.(kovh),
   xaxis = :log,
   yaxis = :log,
   m = 2,
   lab = "Non-linear P(k)",
   ylab = L"P(k)~[h^{-3}~Mpc^3]",
   xlab = L"k~[h~Mpc^{-1}]",
)
p = plot!(kovh, pk.(kovh), ls = :dash, lab = "Linear P(k)")
savefig("pk.pdf")
display(p)


text:
export snr, smooth, smooth!, nextpow2, abs_max!, standardize!, mad, savitsky_golay,
     movingaverage, pws, std_threshold

"""
    snr(A,fs)

Signal to noise ratio of cross-correlations in matrix `A` with `fs'.

Follows method of Clarke et. al, 2011. Measures SNR at each point.
"""
function snr(A::AbstractArray,fs::Real)
    t,N = size(A)
    A_mean = mean(A,dims=2)

    # calculate noise and envelope functions
    sigma = mean(A.^2,dims=2) .- A_mean.^2
    sigma = sqrt.(sigma./ (N-1))
    return sigma
end
snr(C::CorrData) = snr(C.corr,C.fs)

"""
    smooth(A, half_win)

Smooth array `A` with half-window `half_win` (defaults to 3).
"""
function smooth!(A::AbstractArray, half_win::Int=3)
    if ndims(A) == 1
        return movingaverage(A,half_win)
    end

    Nrows, Ncols = size(A)

    for ii = 1:Ncols
        A[:,ii] .= movingaverage(A[:,ii],half_win)
    end
    return nothing
end
smooth(A::AbstractArray,half_win::Int=3) =
      (U = deepcopy(A);smooth!(U,half_win);return U)

"""
  movingaverage(A,half_win)

Smooth array `A` with moving average of half-window `half-win` (defaults to 3.)
"""
function movingaverage(A::AbstractArray, half_win::Int=3)
    prepend!(A,A[1:half_win])
    append!(A,A[end-half_win:end])
    N = length(A)
    B = zeros(eltype(A),N)
    window_len = 2 * half_win + 1
    s = sum(A[1:window_len])
    B[half_win+1] = s
    for ii = half_win+2:N-half_win
        s = s - A[ii-half_win] + A[ii+half_win]
        B[ii] = s
    end
    B ./= window_len
    return B[half_win+1:end-half_win-1]
end

"""
    nextpow2(x)

Returns the next power of 2 of real, positive number `x`.
"""
function nextpow2(x::Real)
    ceil(Int,log2(abs(x)))
end

"""
    abs_max!(A)

Returns array `A` divided by its absolute maximum value.
"""
function abs_max!(A::AbstractArray)
    A .= A ./ maximum(abs.(A),dims=1)
    return nothing
end
abs_max!(C::CorrData) = abs_max!(C.corr)

"""
    standardize!(A)

Demean and standardize array `A` to unit std.
"""
function standardize!(A::AbstractArray)
    A .= (A .- mean(A,dims=1)) ./ std(A,dims=1)
    return nothing
end
standardize!(C::CorrData) = standardize!(C.corr)

"""
    mad(A)

Median Absolute Deviation of array `A`.
MAD = median(|Xi- median(A)|)
"""
function mad(A::AbstractArray)
    return median(abs.(A .- median(A,dims=1)),dims=1)
end
mad(C::CorrData) = mad(C.corr)

"""
    savitsky_golay(x, window, polyOrder, [deriv])

Polynomial smoothing of vector `x` with Savitsky Golay filter.
Polynomial order `polyOrder' must be less than window length `N`.
Modified from https://github.com/BBN-Q/Qlab.jl/blob/master/src/SavitskyGolay.jl
"""
function savitsky_golay(x::Vector, N::Int, polyOrder::Int; deriv::Int=0)
    #Some error checking
    @assert isodd(N) "Window size must be an odd integer."
    @assert polyOrder < N "Polynomial order must me less than window size."

    halfWindow = Int((N-1)/2)

    #Setup the S matrix of basis vectors.
    S = zeros(N, polyOrder+1)
    for ct = 0:polyOrder
        S[:,ct+1] = Array(-halfWindow:halfWindow).^ct
    end

    #Compute the filter coefficients for all orders
    G = S*pinv(S'*S)

    #Slice out the derivative order we want
    filterCoeffs = G[:,deriv+1] * factorial(deriv);

    #Pad the signal with the endpoints and convolve with filter
    paddedX = vcat(x[1]*ones(halfWindow), x, x[end]*ones(halfWindow))
    y = conv(filterCoeffs[end:-1:1], paddedX)

    #Return the valid midsection
    return y[2*halfWindow+1:end-2*halfWindow]
end

"""

  std_threshold(A,thresh)

Returns indices of cols of `A` where max(abs.A[:,col]))/std(A[:,col]) < `max_std`.
"""
function std_threshold(A::AbstractArray, max_std::Float64)
    stds = std(A,dims=1)[1,:]
    maxs = maximum(abs.(A),dims=1)[1,:]
    threshs = maxs ./ stds
    ind = []
    for ii = 1:length(threshs)
        if threshs[ii] < max_std
            append!(ind,ii)
        end
    end
    return ind
end


text:
@testset "Metropolis Sampling" begin

function test_local_energy_change()
    system = produce_example_system()

    for _ in 1:3
        rand!(system)
        for _ in 1:50
            # Pick a random site, try to set it to a random spin
            randsite = CartesianIndex(Tuple(mod1.(rand(Int, 4), size(system))))
            newspin = Sunny._random_spin()

            func_diff = Sunny.local_energy_change(system, randsite, newspin)

            orig_energy = energy(system)
            system[randsite] = newspin
            new_energy = energy(system)

            actual_diff = new_energy - orig_energy

            if !(func_diff ≈ actual_diff)
                return false
            end
        end
    end

    return true
end

@test test_local_energy_change()

"Tests that set_temp!/get_temp behave as expected"
function test_set_get_temp_metropolis()
    system = produce_example_system()

    Random.seed!(1111)
    rand_kTs = rand(100)

    for sampler_type in [MetropolisSampler, IsingSampler]
        sampler = sampler_type(system, 1.0, 1)
        for kT in rand_kTs
            set_temp!(sampler, kT)
            # The approx here is due to numerical issues: the temp returned by
            #  get_temp is actually 1 / (1 / kT) from set_temp!.
            if !(get_temp(sampler) ≈ kT)
                return false
            end
        end
    end

    return true
end

@test test_set_get_temp_metropolis()

end


text:
module Debug

export @debug_mode, enable_debug, disable_debug, enable_trace_logging, disable_trace_logging

# `static_debug_mode` is a flag that enables/disables all debug mode checks
# This defaults to true, but you can disable it either by editing this file, or by
# setting this environment variable when compiling this package. This is useful
# for performance benchmarking and in well-tested production environments.
# e.g. SALSA_STATIC_DEBUG=false
const static_debug_mode = parse(Bool, get(ENV, "SALSA_STATIC_DEBUG", "true"))

"""
@debug_mode expr...

Execute `expr` only when static and runtime debug modes are enabled.
"""
macro debug_mode(instr)
    if static_debug_mode
        quote
            if debug_enabled()
                $(esc(instr))
            end
        end
    else
        :()
    end
end

if static_debug_mode
    # Runtime debug mode controls
    function enable_debug()
        _DBG[] = true
    end
    function disable_debug()
        _DBG[] = false
    end
    debug_enabled() = _DBG[]
    const _DBG = Ref(true)

    # Runtime trace logging controls
    function enable_trace_logging()
        _tracing[] = true
    end
    function disable_trace_logging()
        _tracing[] = false
    end
    trace_logging_enabled() = _tracing[]
    const _tracing = Ref(false)
else
    # Runtime Debugging is disabled.
    function _emit_debug_warning()
        @warn """
            Cannot enable runtime debug statements because debug is disabled statically.
            To enable, reload Salsa after setting `static_debug_mode = true` in:
            $(@__FILE__)
            """
    end

    enable_debug() = _emit_debug_warning()
    disable_debug() = _emit_debug_warning()
    debug_enabled() = false

    enable_trace_logging() = _emit_debug_warning()
    disable_trace_logging() = _emit_debug_warning()
    trace_logging_enabled() = false
end

macro dbg_log_trace(expr)
    if static_debug_mode
        quote
            if trace_logging_enabled()
                $(esc(expr))
            end
        end
    else
        :()
    end
end

end # module Debug


text:
using Revise
using CSV
using Geodesy
using Dates
using Plots
using StatsPlots
using MinimalRides
using MinimalRides: Pos, load_animal_data

filename = "/media/win/Data/Arctic fox Bylot - GPS tracking.csv"
tracks = load_animal_data(filename);
@assert length(tracks) == 20

length(tracks[1].route)
x = tracks[1].route[3].ts - tracks[1].route[2].ts

arr = map(2:length(tracks[1].route)) do i
    ts2 = tracks[1].route[i].ts
    ts1 = tracks[1].route[i - 1].ts
    p2 = tracks[1].route[i].ecef
    p1 = tracks[1].route[i-1].ecef
    return (seconds = Second(ts2 - ts1).value, dist = distance(p2, p1))
end

sum(x -> x.dist == 0., arr)
findall(x -> x.dist == 0., arr)
findall(x -> x.seconds > 1000., arr)

arr[5590:5610]
x = getfield.(tracks[1].route[2000:4000], :ts)
y = getfield.(arr[2000:4000], :dist)
flt = findall(x -> x > 20.0, y)
plot(x[flt], y[flt])
density(getfield.(arr[2000:4000], :dist))
histogram(getfield.(arr[2000:4000], :dist), bins = 1000)

filter(x -> Date(x.ts) == Date("2019-06-11"), tracks[1].route) |> length
minimum(getfield.(arr[1:349], :dist))
maximum(getfield.(arr[1:349], :seconds))
minimum(getfield.(arr[350:350+357], :dist))
maximum(getfield.(arr[350:350+357], :seconds))

plot(map(x -> (x.ll.lat, x.ll.lon), tracks[1].route))
plot!(map(x -> (x.ll.lat, x.ll.lon), tracks[2].route))
plot!(map(x -> (x.ll.lat, x.ll.lon), tracks[3].route))
plot!(map(x -> (x.ll.lat, x.ll.lon), tracks[4].route))
plot!(map(x -> (x.ll.lat, x.ll.lon), tracks[5].route))
plot!(map(x -> (x.ll.lat, x.ll.lon), tracks[6].route))

names = unique(map(x -> x.info.name_id, tracks))
for name in names
    println(name, ": ", filter(x -> x.info.name_id == name, tracks) |> length)
end

filter(x -> x.info.name_id == "BVOB", tracks)

data = CSV.File(filename);
data[1][Symbol("location-long")]
data[1][Symbol("location-lat")]
DateTime(data[1].timestamp[1:19], "yyyy-mm-dd HH:MM:SS")
data[1]

sort(data, by = x -> (x[Symbol("individual-local-identifier")], x.timestamp))

data[1]
Pos(data[1])
Pos(data[1]).ecef


text:
using REALMSupport
using Test, Optim, OffsetArrays, Statistics

@testset "REALMSupport" begin
    include("grid_test.jl")
    include("gauss_fit_test.jl")
    include("register_test.jl")
    include("segment_test.jl")
end


text:
mutable struct OperatorsXMLElement <: MyXMLElement
    el::XMLOrNothing
    els::Vector{<:MyXMLElement}

    OperatorsXMLElement(el::MyXMLElement) = new(nothing, [el])
    OperatorsXMLElement(els::Vector{<:MyXMLElement}) = new(nothing, els)

end

function make_xml(os_el::OperatorsXMLElement)
    for el in os_el.els
        make_xml(el)
    end

    el = new_element(bn.OPERATORS)
    set_attribute(el, bn.ID, bn.OPERATORS)
    for o_el in os_el.els
        if typeof(o_el) <: CompoundOperatorXMLElement
            for xml_el in o_el.els
                add_child(el, xml_el)
            end
        else
            add_child(el, o_el.el)
        end
    end
    os_el.el = el
    return el
end

function merge_operators!(ops1::OperatorsXMLElement, ops2::OperatorsXMLElement)
    ops1.els = [ops1.els; ops2.els]
    return ops1
end

mutable struct FireParameterOperatorXMLElement <: OperatorXMLElement
    el::XMLOrNothing
    param_provider::MyXMLElement
    value::Array{T} where T <: Real
    weight::Float64

    FireParameterOperatorXMLElement(my_xml::MyXMLElement) =
            new(nothing, my_xml, zeros(0), 1.0)
end

function make_xml(fpo_el::FireParameterOperatorXMLElement)

    make_xml(fpo_el.param_provider)

    el = new_element(bn.FIRE_PARAMETER_CHANGED)
    set_attribute(el, bn.WEIGHT, fpo_el.weight)

    if length(fpo_el.value) > 0
        value_string = join(fpo_el.value, ' ')
        set_attribute(el, bn.VALUE, value_string)
    end

    param_el = get_matrix_parameter(fpo_el.param_provider)
    add_ref_el(el, param_el)

    fpo_el.el = el

    return el
end

################################################################################
## scale operator
################################################################################

mutable struct ScaleOperator <: OperatorXMLElement
    el::XMLOrNothing
    param::MyXMLElement
    scale_factor::Float64
    inds::Vector{Int}
    weight::Float64

    function ScaleOperator(param::MyXMLElement, scale_factor::Float64,
                           weight::Float64)
        return new(nothing, param, scale_factor, Int[], weight)
    end
end

function ScaleOperator(param::MyXMLElement)
    return ScaleOperator(param, 0.5, 1.0)
end

function make_xml(so::ScaleOperator)
    el = new_element(bn.SCALE_OPERATOR)
    set_attribute(el, bn.SCALE_FACTOR, so.scale_factor)
    set_attribute(el, bn.WEIGHT, so.weight)
    if length(so.inds) == 0
        add_ref_el(el, so.param)
    else
        mask = zeros(length(so.param))
        mask[so.inds] .= 1.0
        mask_param = MaskedParameter(so.param, mask)
        add_child(el, make_xml(mask_param))
    end

    so.el = el
    return el
end

################################################################################
## shrinkage operators
################################################################################

abstract type CompoundOperatorXMLElement <: OperatorXMLElement end

mutable struct ShrinkageScaleOperators <: CompoundOperatorXMLElement
    els::Vector{XMLOrNothing}
    msl::MatrixShrinkageLikelihoods
    ifxml::IntegratedFactorsXMLElement
    global_weights::Vector{Float64}
    bb_weights::Vector{Float64}
    fix_first::Bool
    fix_globals::Bool

    function ShrinkageScaleOperators(msl::MatrixShrinkageLikelihoods,
                                    ifxml::IntegratedFactorsXMLElement)
        k = get_fac_dim(msl)
        return new(xml_vec(2 * k),
                    msl,
                    ifxml,
                    ones(k),
                    ones(k),
                    false,
                    false)
    end
end

function make_scale_operator(param::XMLElement,
                            scale_factor::Float64,
                            weight::Float64)
    return make_xml(ScaleOperator(param, scale_factor, weight))
end

function make_bb_operator(bb::XMLElement, weight::Float64)
    el = new_element(bn.BAYESIAN_BRIDGE_GIBBS_OP)
    set_attribute(el, bn.WEIGHT, weight)
    add_ref_el(el, bb)
    return el
end

function make_xml(xml::ShrinkageScaleOperators)
    msl = xml.msl
    make_xml(xml.ifxml)
    make_xml(msl)
    k = get_fac_dim(msl)
    if xml.fix_globals
        xml.els[1:k] .= nothing
    else
        if xml.fix_first
            xml.els[1] = nothing
        else
            op = NormalGammaPrecisionOperatorXMLElement(xml.msl, 1)
            make_xml(op)
            xml.els[1] = op.el
            # xml.els[1] = make_scale_operator(msl.gp_els[1], xml.scale_factors[1],
            #                                 xml.scale_weights[1])
        end

        for i = 2:k
            gibbs_provider = MultipilcativeGammaGibbsProvider(xml.msl, i)
            prior_provider = gibbs_provider
            op = NormalGammaPrecisionOperatorXMLElement(gibbs_provider,
                                                        prior_provider)
            make_xml(op)
            xml.els[i] = op.el
            # xml.els[i] = make_scale_operator(msl.mult_els[i - 1],
            #                                 xml.scale_factors[i],
            #                                 xml.scale_weights[i])
        end
    end

    for i = 1:k
        xml.els[i + k] = make_bb_operator(msl.bb_els[i], xml.bb_weights[i])
    end

end

################################################################################
## joint operator
################################################################################

mutable struct JointOperator <: OperatorXMLElement
    el::XMLOrNothing
    operators::Vector{<:OperatorXMLElement}
    weight::Float64

    function JointOperator(operators::Vector{<:OperatorXMLElement},
                           weight::Float64)
        return new(nothing, operators, weight)
    end
end

function make_xml(jo::JointOperator)
    el = new_element(bn.JOINT_OPERATOR)
    for op in jo.operators
        op_el = make_xml(op)
        add_child(el, op_el)
    end
    set_attribute(el, bn.WEIGHT, jo.weight)
    jo.el = el
    return el
end

################################################################################
## rejection operator
################################################################################

mutable struct RejectionOperator <: OperatorXMLElement
    el::XMLOrNothing
    operator::OperatorXMLElement
    condition::Union{String, MyXMLElement}

    function RejectionOperator(op::OperatorXMLElement, condition::Union{String, MyXMLElement})
        return new(nothing, op, condition)
    end
end

function name(::RejectionOperator)
    return bn.REJECTION_OPERATOR
end

function make_xml(ro::RejectionOperator)
    el = new_element(name(ro))
    add_condition(el, ro.condition)
    add_child(el, make_xml(ro.operator))
    ro.el = el
    return el
end

function add_condition(el::XMLElement, condition::String)
    set_attribute(el, bn.CONDITION, condition)
end

function add_condition(el::XMLElement, condition::MyXMLElement)
    c_el = make_xml(condition)
    add_child(el, c_el)
end

function get_parameter(ro::RejectionOperator)
    return get_parameter(ro.operator)
end

mutable struct DescendingAndSpacedCondition <: MyXMLElement
    el::XMLOrNothing
    spacing::Float64
    function DescendingAndSpacedCondition(spacing::Float64)
        return new(nothing, spacing)
    end
end

function make_xml(dsc::DescendingAndSpacedCondition)
    el = new_element(bn.DESCENDING_AND_SPACED)
    set_attribute(el, bn.SPACING, dsc.spacing)
    dsc.el = el
    return el
end


text:
module DamagedShearBand

  using Tensors
  using StaticArrays
  using ForwardDiff
  using Optim

  include("types.jl")
  include("utils.jl")
  include("damaged_rheology.jl")
  include("analysis/analysis.jl")
  include("simulations_helpers.jl")
  include("integrators.jl")

  #export Rheology
  #export set_plane_strain_oop_stress, get_damage_onset

end


text:
import ..SamplingConstraints:
    SequentialSamplingConstraint

import ..InterpolationAndGrids:
    InterpolationGrid

"""
    SequentialInterpolatedResampling{SequentialSamplingConstraint, InterpolationGrid}

Indicates that resampling should be done by first resampling sequentially, then
interpolating the sample to an interpolation grid.

## Fields

- `sequential_constraint::SequentialSamplingConstraint`. The sequential sampling constraint,
    for example `StrictlyIncreasing()`.
- `grid::InterpolationGrid`. The grid onto which the resampled draw (generated according to
    the sequential constraint) is interpolated, for example `RegularGrid(0, 100, 2.5)`.

## Examples

For example, `SequentialInterpolatedResampling(StrictlyIncreasing(), RegularGrid(0:2:100))`
indicates a sequential draw that is then interpolated to the grid 0:2:100.

"""
struct SequentialInterpolatedResampling{S, G} <: AbstractUncertainDataResampling where {
    S <: SequentialSamplingConstraint, G <: InterpolationGrid}
    sequential_constraint::S
    grid::G
end

function Base.show(io::IO, resampling::SequentialInterpolatedResampling{S, G}) where {S, G}
    print(io, "SequentialInterpolatedResampling{$S, $(resampling.grid)}")
end

export SequentialInterpolatedResampling


text:
stime=DateTime(0,7,1)

# time
ΔT=1; #seconds
ndays=10; # duration of simulation in days

n=Int(round(ndays*24*60*60/ΔT)); # number of time steps

dt = ΔT; # time step size [s]

#print times
PT = 0:n:24*60/5;


text:
@inline function add_dbdb_db(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    res = DoubleFloat{T}(add_dddd_dd(HILO(x), HILO(y)))
    return !isnan(HI(res))  ? res : add_dbdb_db_nonfinite(x,y)
end

Base.:(+)(x::DoubleFloat{T}, y::Tuple{T,T}) where {T<:AbstractFloat} = x + DoubleFloat{T}(y[1], y[2])
Base.:(+)(x::Tuple{T,T}, y::DoubleFloat{T}) where {T<:AbstractFloat} = DoubleFloat{T}(x[1], x[2]) + y

@inline function sub_dbdb_db(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    res = DoubleFloat{T}(sub_dddd_dd(HILO(x), HILO(y)))
    return !isnan(HI(res)) ? res : sub_dbdb_db_nonfinite(x,y)
end

Base.:(-)(x::DoubleFloat{T}, y::Tuple{T,T}) where {T<:AbstractFloat} = x - DoubleFloat{T}(y[1], y[2])
Base.:(-)(x::Tuple{T,T}, y::DoubleFloat{T}) where {T<:AbstractFloat} = DoubleFloat{T}(x[1], x[2]) - y

@inline function mul_dbdb_db(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    res = DoubleFloat{T}(mul_dddd_dd(HILO(x), HILO(y)))
    !isnan(HI(res))  ? res : mul_dbdb_db_nonfinite(x,y)
end

Base.:(*)(x::DoubleFloat{T}, y::Tuple{T,T}) where {T<:AbstractFloat} = x * DoubleFloat{T}(y[1], y[2])
Base.:(*)(x::Tuple{T,T}, y::DoubleFloat{T}) where {T<:AbstractFloat} = DoubleFloat{T}(x[1], x[2]) * y

@inline function dvi_dbdb_db(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    res = DoubleFloat{T}(dvi_dddd_dd(HILO(x), HILO(y)))
    !isnan(HI(res)) ? res : dvi_dbdb_db_nonfinite(x,y)
end

Base.:(/)(x::DoubleFloat{T}, y::Tuple{T,T}) where {T<:AbstractFloat} = x / DoubleFloat{T}(y[1], y[2])
Base.:(/)(x::Tuple{T,T}, y::DoubleFloat{T}) where {T<:AbstractFloat} = DoubleFloat{T}(x[1], x[2]) / y

@inline function add_dbdb_db_nonfinite(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    z = HI(x) + HI(y)
    return DoubleFloat{T}(z, zero(T))
end

@inline function sub_dbdb_db_nonfinite(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    z = HI(x) - HI(y)
    return DoubleFloat{T}(z, zero(T))
end

@inline function mul_dbdb_db_nonfinite(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    z = HI(x) * HI(y)
    return DoubleFloat{T}(z, zero(T))
end

@inline function dvi_dbdb_db_nonfinite(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    z = div(HI(x), HI(y))
    return DoubleFloat{T}(z, zero(T))
end


text:
##############
using Revise
using FDMIC_Geodynamics
using LoopVectorization, Printf, Base.Threads, Plots, Revise, LinearAlgebra, Statistics, SparseArrays
##############
function SetMarkers!( p::Markers, params::ModelParameters, domain::ModelDomain )    # Use this function to set up the model geometry
    R  = params.user[1]
    L  = domain.xmax - domain.xmin
    H  = domain.ymax - domain.ymin
    x1 = -1.0
    x2 =  1.0
    for k=1:p.nmark #@tturbo <---------rand does not work in @tturbo!!!!
        p.phase[k] = 1.0
    end
    # for k=1:p.nmark #@tturbo <---------rand does not work in @tturbo!!!!
    #     notch      = ( p.x[k] < 0.0) && (abs(p.y[k] - H/2.0) < R)
    #     in         = notch
    #     p.phase[k] = (in==0)* 1.0 + (in==1)*2.0
    # end
    return
end
##############
function Rheology!( f, Eps, Tau, params, materials, BC, ncx, ncy )
    # Compute effective viscosity for each phase
    f.etac .= 0.0
    f.Kc   .= 0.0
    f.Gc   .= 0.0
    # f.GDam .= 0.0
    # f.lDam .= 0.0
    for m=1:materials.nphase
        @tturbo eta    = materials.G[m]*params.dt
        @tturbo f.etac .+= f.phase_perc[m,:,:] .* eta
        @tturbo f.Kc   .+= f.phase_perc[m,:,:] .* materials.K[m]
        @tturbo f.Gc   .+= f.phase_perc[m,:,:] .* materials.G[m]
        # @tturbo f.GDam .+= f.phase_perc[m,:,:] .* materials.Gc[m]
        # @tturbo f.lDam .+= f.phase_perc[m,:,:] .* materials.l0[m]
    end
    @time CentroidsToVertices!( f.etav, f.etac, ncx, ncy, BC )
    @time CentroidsToVertices!( f.Gv,   f.Gc,   ncx, ncy, BC )
    println("min Eii: ", minimum(Eps.II), " --- max Eii: ", maximum(Eps.II))
    println("min eta: ", minimum(f.etac), " --- max eta: ", maximum(f.etac))
end
export Rheology!
##############
@views function main( )
    # Main routine
    println("#####################")
    println("###### M2Di.jl ######")
    println("#####################")
    # Name of the experiment
    experiment       = "ElasticDamage"
    # Domain
    domain           = ModelDomain()
    domain.xmin      = -0.5e-3
    domain.xmax      =  0.5e-3
    domain.ymin      = -0.0e-3
    domain.ymax      =  1.0e-3
    # Parameters
    params           = ModelParameters()
    params.Ebg       = 2.0e-6      # reference strain rate (>0 horizontal compression)
    # Scales
    scale            = CharScales()
    scale.L          = 1e-3
    scale.T          = 1.0
    scale.t          = 1.0/abs(params.Ebg)
    scale.S          = 1e7
    # Spatial discretisation
    ncx              = 150;  domain.ncx = ncx
    ncy              = 150;  domain.ncy = ncy
    # Time discretisation
    params.advection = 0
    params.nt        = 1#800
    params.dt        = 5.0
    params.Courant   = 0.01  # Courant number
    params.dt_var    = 0     # variable dt
    params.t         = 0.0   # time
    # Boundary conditions
    params.BC_type   = 1     # 1: Pure shear / 2: Simple shear / 3: Simple shear periodic
    params.ALE       = 1     # Deform box
    params.comp      = 1     # Compressible
    # Solver
    params.solver    = 1     # 0: coupled --- 1: decoupled
    params.gamma     = 1e8   # Penalty factor
    params.Dir_scale = 1.0   # Dirichlet scaling factor
    # Non-linear iterations
    params.niter_nl  = 5    # max. number of non-linear iterations
    params.tol_nl    = 1e-4  # non-linear tolerance
    params.JFNK      = 0     # Jacobian-Free Newton_Krylov
    # Visualisation
    params.show_figs = 1     # activates visualisation...
    params.nout      = 10    # ... every nout
    # Material setup
    params.user[1]    = 15e-6/2.0 # layer half thickness
    params.user[2]    = 1e-3 # numerical parameter for phase-field algorithm
    # Material parameters
    materials         = MaterialParameters()
    materials.nphase  = 2
    materials.Tref    = zeros(Float64, materials.nphase)
    materials.eta0    = zeros(Float64, materials.nphase)
    materials.n       = zeros(Float64, materials.nphase)
    materials.G       = zeros(Float64, materials.nphase)
    materials.K       = zeros(Float64, materials.nphase)
    # Constant material parameters (all phases)
    materials.gDam    = 2.70e3#zeros(Float64, materials.nphase)
    materials.lDam    = 15.0e-6#zeros(Float64, materials.nphase)
    materials.eDam    = 1.0e3#zeros(Float64, materials.nphase)
    # Material 1
    materials.Tref[1] = 1e-1  # reference flow stress
    materials.n[1]    = 1.0
    materials.K[1]    = 174.9e9
    materials.G[1]    = 80.77e9
    materials.eta0[1] = materials.Tref[1] * (1.0/2.0) * abs(params.Ebg)^(-1.0/materials.n[1])  # matrix viscosity
    # Material 2
    materials.Tref[2] = 1     # reference flow stress
    materials.n[2]    = 1.0
    materials.K[2]    = 174.9e8
    materials.G[2]    = 80.77e8
    materials.eta0[2] = materials.Tref[2] * (1.0/2.0) * abs(params.Ebg)^(-1.0/materials.n[2]) # inclusion viscosity
    # Scaling
    domain.xmin     /= scale.L
    domain.xmax     /= scale.L
    domain.ymin     /= scale.L
    domain.ymax     /= scale.L
    materials.gDam /= (scale.S*scale.L)
    materials.eDam /= (scale.S*scale.t)
    materials.lDam /=  scale.L
    materials.G    ./= scale.S
    materials.K    ./= scale.S
    params.Ebg      /= (1.0/scale.t)
    params.dt       /= scale.t
    params.user[1]  /=  scale.L
    # Initialize coordinates
    domain.L0 = domain.xmax-domain.xmin
    GenerateMesh!( domain )
    dx, dy = domain.dx, domain.dy
    # BC definition
    BC = BoundaryConditions()
    if params.BC_type==1
        ## Pure shear
        BC.Vx.type_W = 1;     BC.Vx.Dir_W =-0.0*params.Ebg*domain.xmin*ones(ncy+2)
        BC.Vx.type_E = 1;     BC.Vx.Dir_E =-0.0*params.Ebg*domain.xmax*ones(ncy+2)
        BC.Vx.type_S =22;     BC.Vx.Dir_S = 0.0*ones(ncx+1)
        BC.Vx.type_N =22;     BC.Vx.Dir_N = 0.0*ones(ncx+1)
        BC.Vy.type_W =22;     BC.Vy.Dir_W = 0.0*ones(ncy+1)
        BC.Vy.type_E =22;     BC.Vy.Dir_E = 0.0*ones(ncy+1)
        BC.Vy.type_S = 1;     BC.Vy.Dir_S = params.Ebg*domain.ymin*ones(ncx+2)
        BC.Vy.type_N = 1;     BC.Vy.Dir_N = params.Ebg*domain.ymax*ones(ncx+2)
    elseif params.BC_type==2
        ## Simple shear
        BC.Vx.type_W = 1;     BC.Vx.Dir_W = 0.0*ones(ncy+2)
        BC.Vx.type_E = 1;     BC.Vx.Dir_E =-0.0*ones(ncy+2)
        BC.Vx.type_S =11;     BC.Vx.Dir_S =-1.0*ones(ncx+1)
        BC.Vx.type_N =11;     BC.Vx.Dir_N = 1.0*ones(ncx+1)
        BC.Vy.type_W =22;     BC.Vy.Dir_W = 0.0*ones(ncy+1)
        BC.Vy.type_E =22;     BC.Vy.Dir_E = 0.0*ones(ncy+1)
        BC.Vy.type_S = 1;     BC.Vy.Dir_S =-0.0*ones(ncx+2)
        BC.Vy.type_N = 1;     BC.Vy.Dir_N = 0.0*ones(ncx+2)
    elseif params.BC_type==3
        ## Simple shear periodic in x
        BC.periodix  = 1
        BC.Vx.type_W = 0;     BC.Vx.Dir_W = 0.0*ones(ncy+2)
        BC.Vx.type_E = 0;     BC.Vx.Dir_E =-0.0*ones(ncy+2)
        BC.Vx.type_S =11;     BC.Vx.Dir_S =-1.0*ones(ncx+1)
        BC.Vx.type_N =11;     BC.Vx.Dir_N = 1.0*ones(ncx+1)
        BC.Vy.type_W = 0;     BC.Vy.Dir_W = 0.0*ones(ncy+1)
        BC.Vy.type_E = 0;     BC.Vy.Dir_E = 0.0*ones(ncy+1)
        BC.Vy.type_S = 1;     BC.Vy.Dir_S =-0.0*ones(ncx+2)
        BC.Vy.type_N = 1;     BC.Vy.Dir_N = 0.0*ones(ncx+2)
    end
    BC_Dam = Thermal_BC()
    BC_Dam.type_W = 2;     BC_Dam.Dir_W = 0.0*ones(ncy)
    BC_Dam.type_E = 2;     BC_Dam.Dir_E = 0.0*ones(ncy)
    BC_Dam.type_S = 2;     BC_Dam.Dir_S = 0.0*ones(ncx)
    BC_Dam.type_N = 2;     BC_Dam.Dir_N = 0.0*ones(ncx)
    BC_Dam.Pinned = zeros(Int64,ncx,ncy)
    # Allocate tables
    fields           = Fields2D()
    fields.etac      = zeros(Float64, ncx  , ncy  )
    fields.etav      = zeros(Float64, ncx+1, ncy+1)
    fields.Kc        = zeros(Float64, ncx  , ncy  )
    fields.Gc        = zeros(Float64, ncx  , ncy  )
    fields.Gv        = zeros(Float64, ncx+1, ncy+1)
    fields.We0       = zeros(Float64, ncx  , ncy  )
    fields.We        = zeros(Float64, ncx  , ncy  )
    fields.rhoc      = zeros(Float64, ncx  , ncy  )
    fields.Pc0       = zeros(Float64, ncx+0, ncy+0)
    fields.Pc        = zeros(Float64, ncx+0, ncy+0)
    fields.Vx        = zeros(Float64, ncx+1, ncy+2) # !!! GHOST ROWS
    fields.Vy        = zeros(Float64, ncx+2, ncy+1) # !!! GHOST COLUMNS
    Tau0             = Tensor2D( zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+1, ncy+1), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0) )
    Tau              = Tensor2D( zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+1, ncy+1), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0) )
    Eps              = Tensor2D( zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+1, ncy+1), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0) )
    Str0             = Tensor2D( zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+1, ncy+1), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0) )
    Str              = Tensor2D( zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+1, ncy+1), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0) )
    if BC.periodix==0
        fields.Fx    = zeros(Float64, ncx+1, ncy+0)
    else
        fields.Fx    = zeros(Float64, ncx+0, ncy+0)
    end
    fields.Fy        = zeros(Float64, ncx+0, ncy+1)
    fields.Fp        = zeros(Float64, ncx  , ncy  )
    fields.Damc      = zeros(Float64, ncx  , ncy  )
    fields.Damv      = zeros(Float64, ncx+1, ncy+1)
    fields.phiDam    = zeros(Float64, ncx,   ncy  )
    fields.phiDam0   = zeros(Float64, ncx,   ncy  )
    fields.FDam      = zeros(Float64, ncx,   ncy  )
    # Number equations
    println("Numbering")
    @time NumberingStokes!( fields, BC, domain )
    # Initialise markers
    p                = Markers()
    p.nmx            = 4            # 2 marker per cell in x
    p.nmy            = 4            # 2 marker per cell in y
    InitialiseMarkerPositions!( p, domain )
    # Over allocate markers
    p.phase          = zeros(Float64, p.nmark_max)
    p.T              = zeros(Float64, p.nmark_max)
    p.cellx          = zeros(Int64,   p.nmark_max)
    p.celly          = zeros(Int64,   p.nmark_max)
    # Marker 2 cell mapping
    fields.mpc           = zeros(Float64,(ncx  ,ncy  )) # markers per cell
    fields.mpc_th        = zeros(Float64,(nthreads(), ncx  ,ncy   )) # markers per cell
    fields.phase_perc    = zeros(Float64,(materials.nphase, ncx  ,ncy  )) # markers per cell
    fields.phase_perc_th = zeros(Float64,(nthreads(), materials.nphase, ncx  ,ncy   )) # markers per cell
    # Initial configuration
    SetMarkers!( p, params, domain )  # Define phase on markers
    SetInitialVelocity2!( fields.Vx, fields.Vy, BC, domain )
    fields.Vxc    = 0.5.*(fields.Vx[1:end-1,2:end-1] .+ fields.Vx[2:end-0,2:end-1])
    fields.Vyc    = 0.5.*(fields.Vy[2:end-1,1:end-1] .+ fields.Vy[2:end-1,2:end-0])
    # Initial phase field
    BC_Dam.Pinned = zeros(Int64,ncx,ncy)
    for i=1:ncx
        for j=1:ncy
            if (domain.xc[i] < 0.0) && (abs(domain.yc[j] - (domain.ymax-domain.ymin)/2.0) < params.user[1] )
                BC_Dam.Pinned[i,j] = 1
            end
        end
    end
    fields.phiDam[BC_Dam.Pinned.==1] .= 1.0
    @time CountMarkers3!( p, fields, materials, domain )
    @time Rheology!( fields, Eps, Tau, params, materials, BC, ncx, ncy )
    InitialDamageResidual!( fields, BC_Dam, domain, materials )
    K_Dam = InitialDamageAssembly( fields, BC_Dam, domain, materials )
    Kc    = cholesky(K_Dam)
    dd    = Kc\fields.FDam[:]
    fields.phiDam .-= dd[fields.NumP]
    InitialDamageResidual!( fields, BC_Dam, domain, materials )
    # Visualisation
    viz_directory = string( "Figures_", experiment )
    if isdir( viz_directory ) == false
        mkdir( viz_directory )
    end
    path = string( "./", viz_directory, "/" )
    anim = Plots.Animation( path, String[] )
    # Storage
    displ  = 0
    Wev    = zeros(Float64, params.nt)
    Force  = zeros(Float64, params.nt)
    Displ  = zeros(Float64, params.nt)
    resxv  = zeros(Float64, params.niter_nl)
    # TIME LOOP
    for it=1:params.nt
        @printf("-------------------------------------\n")
        @printf("---------- Time step #%04d ----------\n", it)
        @printf("-------------------------------------\n")
        # History
        fields.We0 .= fields.We
        fields.phiDam0 .= fields.phiDam
        fields.Pc0 .= fields.Pc
        Tau0.xx    .= Tau.xx
        Tau0.yy    .= Tau.yy
        Tau0.zz    .= Tau.zz
        Tau0.xy    .= Tau.xy
        Str0.xx    .= Str.xx
        Str0.yy    .= Str.yy
        Str0.zz    .= Str.zz
        Str0.xy    .= Str.xy
        # Count markers and reseed
        @time CountMarkers3!( p, fields, materials, domain )
        @time ReseedMarkers2!( p, fields.mpc, domain)
        nmark = p.nmark
        @time CountMarkers3!( p, fields, materials, domain )
        # Define dt
        if params.dt_var == 1
            maxVx     = maximum(abs.(fields.Vx))
            maxVy     = maximum(abs.(fields.Vy))
            params.dt = params.Courant*(maximum([domain.dx,domain.dy])/maximum([maxVx,maxVy]))
            @printf("dt = %2.2e\n", params.dt)
        end
        # Update cell info on markers
        @time LocateMarkers2!(p, domain)
        # Non-linear iterations
        iterations = 0
        for iter=1:params.niter_nl
            iterations+=1
            @printf("-------------------------------------------------\n")
            @printf("---------- Iteration #%04d (step #%04d)----------\n", iter, it)
            @printf("-------------------------------------------------\n")
            # Evaluate residuals
            println("Residuals")
            @time SetBCs( fields.Vx, fields.Vy, BC )
            @time StrainRate!( Eps, fields.Vx, fields.Vy, domain )
            fields.Damc .= (1.0 .- fields.phiDam ).^2 .+ 1e-6
            @time CentroidsToVertices!( fields.Damv, fields.Damc, ncx, ncy, BC )
            @time Rheology!( fields, Eps, Tau, params, materials, BC, ncx, ncy )
            @time StressVE!( Tau, Tau0, Eps, fields, params.dt )
            @time StrainEnergy!( fields.We, Tau, Eps, Str0, Str, fields.Pc, params.dt, fields, dx, dy )
            @time ResidualsCompDam!( fields, Tau, Eps, BC, domain, params )
            @time DamageResidual!( fields, BC_Dam, domain, materials, params.dt )
            println("|Fx|   = ", norm(fields.Fx)/length(fields.Fx))
            println("|Fy|   = ", norm(fields.Fy)/length(fields.Fy))
            println("|Fp|   = ", norm(fields.Fp)/length(fields.Fp))
            println("|FDam| = ", norm(fields.FDam)/length(fields.FDam))
            if ( (norm(fields.Fx)/length(fields.Fx) < params.tol_nl) & (norm(fields.Fy)/length(fields.Fy) < params.tol_nl) & (norm(fields.Fp)/length(fields.Fp) < params.tol_nl)  && (norm(fields.FDam)/length(fields.FDam) < params.tol_nl))
                println("Mechanical solver converged")
                break
            end
            resxv[iter] = norm(fields.Fx)/length(fields.Fx)

            # Assemble Stokes matrices
            println("Assembly")
            @time Kuu, Kup, Kpu = StokesAssembly( BC, fields, params.Dir_scale, domain )
            # Call solver
            println("Solver")
            @time StokesSolver!( fields, params, domain, materials, BC, Kuu, Kup, Kpu )

            K_Dam = DamageAssembly( fields, BC_Dam, domain, materials, params.dt )
            # dd = K_Dam\fields.FDam[:]
            # fields.phiDam .-= dd[fields.NumP]
            # display(Plots.spy(K_Dam))
            Kc = cholesky(K_Dam)
            dd = Kc\fields.FDam[:]
            fields.phiDam .-= dd[fields.NumP]

            # Evaluate residuals
            println("Residuals")
            @time SetBCs( fields.Vx, fields.Vy, BC )
            @time StrainRate!( Eps, fields.Vx, fields.Vy, domain )
            @time StressVE!( Tau, Tau0, Eps, fields, params.dt )
            @time ResidualsCompDam!( fields, Tau, Eps, BC, domain, params )
            @time DamageResidual!( fields, BC_Dam, domain, materials, params.dt )

            println("|Fx|   = ", norm(fields.Fx)/length(fields.Fx))
            println("|Fy|   = ", norm(fields.Fy)/length(fields.Fy))
            println("|Fp|   = ", norm(fields.Fp)/length(fields.Fp))
            println("|FDam| = ", norm(fields.FDam)/length(fields.FDam))
        end
        params.t += params.dt
        # Advect particles
        if params.advection == 1
            @time RungeKutta2!( p, fields, params, BC, domain )
            # Deform box
            if params.ALE==1
                dx, dy = PureShearBoxUpdate!( domain, BC, params )
                SetInitialVelocity2!( fields.Vx, fields.Vy, BC, domain )
            end
        end
        # Postprocessing
        # Wev[it]   = fields.Kc[1,1]#mean(Tau.II)
        Wev[it]   = mean(fields.We)
        Force[it] = sum( Tau.yy[:,end] .- fields.Pc[:,end])
        Force[it]+= sum( Tau.yy[:,1] .- fields.Pc[:,end])
        Force[it]+= sum( Tau.xx[1,:] .- fields.Pc[1,:])
        Force[it]+= sum( Tau.xx[end,:] .- fields.Pc[end,:])
        displ    += maximum(abs.(fields.Vy[:,end]).*params.dt)
        Displ[it] = displ

        DIV = maximum(fields.Vy)/(domain.ymax-domain.ymin)
        P   = -mean(fields.Kc)*DIV*params.t
        EXXd = 0.0 - 1.0/3.0*DIV
        EYYd = DIV - 1.0/3.0*DIV
        TXX = 2.0*mean(fields.Gc) * EXXd*params.t
        TYY = 2.0*mean(fields.Gc) * EYYd*params.t
        SXX = -P + TXX
        SYY = -P + TYY
        EYY = DIV*params.t
        WE  = EYY*SYY
        println("Displacement: ", displ * scale.L * 1e3, " --- Max incremental displacement: ", maximum((fields.Vy).*params.dt * scale.L * 1e3 ),  " mm")
        println("dx = ", domain.dx*scale.L*1e3, " mm --- l0 = ",materials.lDam*scale.L*1e3, " mm")
        println("min/max div: ", minimum(Eps.div)*(1.0/scale.t), " ", maximum(Eps.div)*(1.0/scale.t))
        println("bulk divergence: ",  DIV*(1.0/scale.t) )
        println("min/max P: ", minimum(fields.Pc)*scale.S, " ", maximum(fields.Pc)*scale.S)
        println("numerical pressure: ", mean(fields.Pc)*scale.S, " --- expected: ", P*scale.S)
        println("numerical Txx:      ", mean(Tau.xx)*scale.S, " --- expected: ", TXX*scale.S)
        println("numerical Tyy:      ", mean(Tau.yy)*scale.S, " --- expected: ", TYY*scale.S)
        println("numerical Sxx:      ", mean(Tau.xx.-fields.Pc)*scale.S, " --- expected: ", SXX*scale.S)
        println("numerical Syy:      ", mean(Tau.yy.-fields.Pc)*scale.S, " --- expected: ", SYY*scale.S)
        println("numerical We:       ", mean(fields.We)*scale.S, " --- expected: ", WE*scale.S)

        # Visualisation
        if params.show_figs==1 && ( mod( it, params.nout )==0 || it==1 )
            l = @layout [ [a; b] c]
            # p1 = Plots.heatmap( domain.xc, domain.yc, (Tau.xx .- 0.0*fields.Pc)*scale.S)
            # p2 = Plots.heatmap( domain.xc, domain.yc, (Tau.xy_c .- 0.0*fields.Pc)*scale.S)
            # p2 = Plots.heatmap( domain.xc, domain.yc, Eps.div*(1.0/scale.t) )
            p1 = Plots.heatmap( domain.xv, domain.yce, Array(fields.Vx*1e9*scale.L/scale.t)', aspect_ratio=1, xlims=(domain.xmin, domain.xmax), ylims=(domain.ymin, domain.ymax), c=Plots.cgrad(:roma, rev = true), title="Vx [nanom/s]" )
            p2 = Plots.heatmap( domain.xce, domain.yv, Array(fields.Vy*1e9*scale.L/scale.t)', aspect_ratio=1, xlims=(domain.xmin, domain.xmax), ylims=(domain.ymin, domain.ymax), c=Plots.cgrad(:roma, rev = true), title="Vy [nanom/s]" )
            # p1 = Plots.heatmap( domain.xc, domain.yc, Array(fields.We)'*scale.S, aspect_ratio=1, xlims=(domain.xmin, domain.xmax), ylims=(domain.ymin, domain.ymax), c=Plots.cgrad(:roma, rev = true), title="We" )
            p3 = Plots.heatmap( domain.xc, domain.yc, Array(fields.phiDam)', aspect_ratio=1, xlims=(domain.xmin, domain.xmax), ylims=(domain.ymin, domain.ymax), c=Plots.cgrad(:roma, rev = true), title="phi, niter = $iterations")
            p4 = Plots.heatmap( domain.xc, domain.yc, Array(Tau.II)'./1e6*scale.S, aspect_ratio=1, xlims=(domain.xmin, domain.xmax), ylims=(domain.ymin, domain.ymax), c=Plots.cgrad(:roma, rev = true), title="Sii [MPa], nt = $it" )
            # p4 = Plots.heatmap( domain.xc, domain.yc, Array(Str.II)', aspect_ratio=1, xlims=(domain.xmin, domain.xmax), ylims=(domain.ymin, domain.ymax), c=Plots.cgrad(:roma, rev = true), title="Strii" )

            # p3 = Plots.heatmap( domain.xc, domain.yc, Array(Tau.II*scale.S)', aspect_ratio=1, xlims=(domain.xmin, domain.xmax), ylims=(domain.ymin, domain.ymax), c=Plots.cgrad(:roma, rev = true), title="Tii" )
            # p3 = Plots.plot( Displ[1:it]*scale.L*1e3, Force[1:it]*(scale.S), markershape=:circle)
            # p3 = Plots.plot( Displ[1:it]*scale.L*1e3, Wev[1:it]*scale.S, markershape=:circle)
            # p3 = Plots.plot( [1:iterations], log10.(resxv[1:iterations]), markershape=:circle)
            display( Plots.plot( p1, p2, p3, p4, dpi=200 ) ); Plots.frame(anim)
            # display( Plots.plot( p1, p2, p3, dpi=200, layout=l ) ); Plots.frame(anim)
        end
    end
    if params.show_figs==1 Plots.gif(anim, string( path, experiment, ".gif" ), fps = 15) end
    return
end

for it=1:1
   @time main()
   @time main()
end


text:
@defstruct LookupTableLayer Layer (
    name :: AbstractString = "lookup-table",
    init :: Initializer = ConstantInitializer(0),
    lr   :: AbstractFloat = 1.0,
    (output_dim :: Int = 0, output_dim > 0),
    (tops :: Vector{Symbol} = Symbol[], length(tops) == 1),
    (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) == length(tops) == 1),
    (n_class :: Int = 0, n_class > 0)
)

@characterize_layer(InnerProductLayer,
  can_do_bp  => true,
  has_param  => true,
)

type LookupTableLayerState <: LayerState
    layer      :: LookupTableLayer
    blobs      :: Vector{Blob}
    blobs_diff :: Vector{Blob}
    parameters :: Vector{Parameter}
    frozen     :: Bool
end

function setup(backend::Backend, layer::LookupTableLayer, inputs::Vector{Blob},
            diffs::Vector{Blob})
    data_type = eltype(inputs[1])
    out_dim = layer.output_dim
    batch_size = get_num(inputs[1])

    blobs = [make_blob(backend, data_type, out_dim, batch_size)]
    blobs_diff = [make_blob(backend, data_type, out_dim, batch_size)]

    parameters = [make_parameter(backend, "lookup", data_type,
                           (out_dim, layer.n_class),
                           layer.init, layer.lr)]
    erase!(parameters[1].gradient)

    return LookupTableLayerState(layer, blobs, blobs_diff, parameters)
end

function shutdown(backend::Backend, state::LookupTableLayerState)
    map(destroy, state.blobs)
    map(destroy, state.blobs_diff)
    map(destroy, state.parameters)
end

function forward(backend::CPUBackend, state::LookupTableLayerState,
            inputs::Vector{Blob})
    output = state.blobs[1].data
    W = state.parameters[1].blob
    input = inputs[1].data

    for i in 1:length(input)
        idx = convert(Int, input[i])
        output[:, i] = W[:, idx]
    end
end

function backward(backend::Backend, state::LookupTableLayerState,
            inputs::Vector{Blob}, diffs::Vector{Blob})
    input = inputs[1].data
    ∇W = state.parameters[1].gradient
    diff = diffs[1].data

    if !state.frozen
        for i in 1:length(input)
            idx = convert(Int, input[i])
            fill!(∇W[:, idx].data, 0.)
        end

        for i in 1:length(input)
            idx = convert(Int, input[i])
            ∇W[:, idx] += diff[:, i]
        end
    end
end


text:
using Random
include("../simulation/simulation.jl")

function move_causes_collision(new_x, new_y, simulation::Simulation)
    for other_bi in 1:length(simulation.bikes)
        if new_x == simulation.bikes[other_bi].x_position && new_y == simulation.bikes[other_bi].y_position
            return true
        end
    end
    return false
end

function get_new_random_position(current)
    if bitrand()[1]
        possible_new_position =  current + rand((-1, 1))
        if possible_new_position > -1 && possible_new_position < 100 + 1
            return possible_new_position
        end
    end
    return current
end

function move_bike(bike::Bike, simulation::Simulation)
    new_x = get_new_random_position(bike.x_position)
    new_y = get_new_random_position(bike.y_position)

    collision_occured = move_causes_collision(new_x, new_y, simulation)

    if !collision_occured
        bike.x_position = new_x
        bike.y_position = new_y
    end
end


text:
module Runtests

# Run all tests with SilentNoPlot (so not plots)

using TinyModia
using Test

@time @testset verbose=true "TinyModia (with SilentNoPlot)" begin
    usePlotPackage("SilentNoPlot")
    include("include_all.jl")
    usePreviousPlotPackage()
end

end


text:
# check polylog reciprocal component
using Polylogarithms
using BenchmarkTools
using DataFrames, CSV
using PyPlot
using Printf
include("utilities.jl")

series_1 = Polylogarithms.polylog_series_1
series_2 = Polylogarithms.polylog_series_2
series_3 = Polylogarithms.polylog_series_3
polylog_reciprocal = Polylogarithms.polylog_reciprocal
L = Symbol("Li_s(z)")

# input data from Mathematica and reparse into complex numbers
for C=1:6
    # C = 1
    println("C = $C")
    filename = @sprintf("../data/polylog_test_data_c_%d.csv", C)
    data1 = CSV.read(filename; delim=",", type=String)
    #    has trouble reading in numbers like "2." so read all into strings, and parse
    data1[!,:s] = parse.(Complex{Float64}, data1[!,:s] )
    data1[!,:r] = parse.(Float64, data1[!,:r] )
    data1[!,:theta] = parse.(Float64, data1[!,:theta] )
    data1[!,:z] = parse.(Complex{Float64}, data1[!,:z] )
    data1[!,L] = parse.(Complex{Float64}, data1[!,L] )
    # check that complex numbers are being read in correctly
    z = data1[!,:z]
    z2 = data1[!,:r] .* exp.( im* data1[!,:theta] )
    error = data1[!,:z] .- z2
    maximum_abs_error = maximum(abs.(error))
    # [data1[1:10,:z] parse.(Complex{Float64}, data1[1:10,:z] ) z2[1:10] ]
    # [data1[1:10,:z]  z2[1:10] ]
    # println("max error in parsing is ", maximum_abs_error)

    # figure(1)
    # plot(real.(z), imag.(z), "o")
    # axis("equal")

    m = size(data1,1)
    Li = data1[!,L]
    s = data1[!,:s]
    r = unique( data1[!,:r])
    su = unique(s)
    if length(su) > 1
        error()
    else
        su = su[1]
        println("")
        @printf("  s = %f + %f i", real(su), imag(su))
    end
    S1 = zeros(Complex{Float64}, m)
    S2 = zeros(Complex{Float64}, m)
    n1 = zeros(Int64, m)
    error1 = zeros(Float64, m)
    rel_error1 = zeros(Float64, m)
    n2 = zeros(Int64, m)
    error2 = zeros(Float64, m)
    rel_error2 = zeros(Float64, m)

    for i=1:m
        # print(".")
        result1 = polylog_reciprocal(s[i], z[i])
        S1[i] = result1[1]
        n1[i] = result1[2]
        error1[i] = abs( S1[i] - Li[i]  )
        rel_error1[i] = error1[i]/abs( Li[i] )

        if real(s[i])>0 && abs(real(s[i]) - round(s[i])) < 1.0e-3
            rel_error2[i] = NaN
        else
            result2 = series_2(s[i], z[i])
            S2[i] = result2[1]
            n2[i] = result2[2]
            error2[i] = abs( S2[i] - Li[i]  )
            rel_error2[i] = error2[i]/abs( Li[i] )
        end
    end
    println("   max abs. error1 = $(maximum( abs.(error1) ))")
    println("   max abs. error2 = $(maximum( abs.(error2) ))")

    fig = figure(@sprintf("../data/polylog_bench_a_%02d.csv", C), figsize=(10,8))
    clf()
    title( @sprintf("s = %f + %f i", real(su), imag(su)) )
    d = 0.01
    ms = 3

    subplot(221)
    semilogy( data1[!,:r] .- d, rel_error1, "o"; markersize=ms)
    semilogy( data1[!,:r] .+ d, rel_error2, "rd"; markersize=ms)
    xlabel("|z|")
    ylabel("relative absolute error")
    plot([minimum(r), maximum(r)], [1,1]*Polylogarithms.default_accuracy)
    # xlim([0, 1.15])
    # ylim([1.0-e15, 1.0e-10])

    subplot(223)
    semilogy( data1[!,:r] .- d, n1, "o"; markersize=ms)
    semilogy( data1[!,:r] .+ d, n2, "rd"; markersize=ms)
    # θ = unique(data1[!,:theta])
    # for i=1:length(θ)
    #     k = findall(data1[!,:theta] .== θ[i])
    #     plot(data1[k,:r] .+ d, n2[k], "r-"; linewidth=0.5)
    #     text(data1[k[end],:r]+2*d, n2[k[end]], "$(θ[i]/π)"; verticalalignment="center", fontsize=9)
    # end
    # text(0.92, 65, "arg(z)/π"; verticalalignment="center")
    xlabel("|z|")
    ylabel("number of terms")
    plot([minimum(r), maximum(r)], [1,1]*Polylogarithms.default_max_iterations, "--")
    # xlim([0, 1.15])

    subplot(222)
    semilogy( data1[!,:theta] ./ π .- d, rel_error1, "o"; label="Reciprocal + Series 1", markersize=ms)
    semilogy( data1[!,:theta] ./ π .+ d, rel_error2, "rd"; label="Series 2", markersize=ms)
    xlabel("arg(z)/π")
    ylabel("relative absolute error")
    plot([-0.05,1.05], [1,1]*Polylogarithms.default_accuracy)
    legend()

    subplot(224)
    semilogy( data1[!,:theta] ./ π .- d, n1, "o"; markersize=ms)
    semilogy( data1[!,:theta] ./ π .+ d, n2, "rd"; markersize=ms)
    xlabel("arg(z)/π")
    ylabel("number of terms")
    plot([-0.05,1.05], [1,1]*Polylogarithms.default_max_iterations, "--")

    # savefig(@sprintf("plots/polylog_bench_a_%02d.svg", C))
    savefig(@sprintf("plots/polylog_bench_c_%02d.pdf", C); bbox_inches="tight")

end


text:
"""
    Context{T <: Params}
    CompressionContext = Context{CompressionParams}
    DecompressionContext = Context{DecompressionParams}

# Constructors

    Context{T}(params::T) where {T <: Params}

Create context.
- If `params` is `CompressionParams` then this is `CompressionContext`
- If `params` is `DecompressionParams` then this is `DecompressionContext`

    CompressionContext(;kwargs...)

Create compression context `kwargs` passed to CompressionParams constructor

    DecompressionContext(;kwargs...)

Create decompression context `kwargs` passed od DecompressionParams consturctor
"""
mutable struct Context{T <: Params}
    params ::T
    context ::Ptr{Lib.blosc2_context}
    function Context(params::T) where {T <: Params}
        result = new{T}(params, _create_context(params))
        finalizer(result) do c
            Lib.blosc2_free_ctx(c.context)
        end
    end
    Context{T}(params::T) where {T <: Params} = Context(params)
    Context{T}(;kwargs...) where {T <: Params} = Context(T(;kwargs...))
end

const CompressionContext = Context{CompressionParams}
const DecompressionContext = Context{DecompressionParams}

"""
    CompressionContext(::Type{T}; kwargs...)

Create `CompressionContext` for compressing vectors with element type `T`.
This is equivalent to `CompressionContext(CompressionParams(T; kwargs...))``
"""
CompressionContext(::Type{T}; kwargs...) where {T} = CompressionContext(CompressionParams(T; kwargs...))

_create_context(params::CompressionParams) = Lib.blosc2_create_cctx(make_cparams(params))
_create_context(params::DecompressionParams) = Lib.blosc2_create_dctx(make_dparams(params))


text:
path = "C:\\Users\\vsger\\Downloads\\rosalind_revp.txt"

function readdata()
	a = readlines(path)
	dna = ""
	for i in 2:length(a)
		dna *= a[i]
	end
	return dna
end

function complement(a)
	if a == "A"
		return "T"
	elseif a == "T"
		return "A"
	elseif a == "G"
		return "C"
	elseif a == "C"
		return "G"
	end
end

function isreversepallindrome(s)
	a = join(map(x->complement(x),reverse(split(s,""))))
	return s == a
end

function main()
	dna = readdata()
	i = 1
	j = 4

	while i < length(dna)
		while j <= min(12,length(dna)-i+1)

			if isreversepallindrome(dna[i:(i+j-1)])
				println(i," ",j)
			end
			j += 1
		end
		j = 4
		i += 1
	end
end
main()

max_stars_repo_name: Ved-Mahajan/JustForFun

max_issues_repo_name: Ved-Mahajan/JustForFun

max_forks_repo_name: Ved-Mahajan/JustForFun


text:
using JSON
using Interact
using DataStructures
using Plots

macro dont_manipulate(args...)
  n = length(args)
  @assert 1 <= n <= 2

  expr = args[n]
  if expr.head != :for
      error("@manipulate syntax is @manipulate for ",
            " [<variable>=<domain>,]... <expression> end")
  end

  cur_block = expr.args[2]

  if expr.args[1].head == :block
    cur_bindings = expr.args[1].args
  else
    cur_bindings = [expr.args[1]]
  end

  cur_defaults = OrderedDict()
  cur_ranges = []

  for cur_binding in cur_bindings
    @assert cur_binding.head == :(=)
    cur_symbol, cur_range = cur_binding.args

    isa(cur_range, Expr) && ( cur_range = eval(cur_range) )
    isa(cur_range, Bool) && ( cur_range = [cur_range, !cur_range] )

    isa(cur_range, Symbol) && ( cur_range = eval(cur_range) )

    if isa(cur_range, Widget)
      work_default = cur_range.output.val
      work_range = cur_range.components[:formatted_vals]

      if !isa(work_default, AbstractString)
        work_range = map(
          work_range_value -> parse(typeof(work_default), work_range_value),
          work_range
        )
      end

      push!(cur_ranges, (string(cur_symbol), work_range))
      cur_defaults[string(cur_symbol)] = work_default
    else
      push!(cur_ranges, (string(cur_symbol), [collect(cur_range)...]))
    end
  end

  plot_tree = build_tree(
    compile_tree_dict(cur_block, cur_ranges)
  )

  for cur_node in plot_tree.nodes
    isa(cur_node, TreeNode) && continue
    cur_plot = eval(cur_node.value)

    IJulia.clear_output(true)
    display(cur_plot)

    cur_node.value = JSON.parse(json(
      Plots.plotlyjs_syncplot(cur_plot).plot
    ))
  end

  plot_json = tree_to_dict(plot_tree.root)

  json_dict = OrderedDict(
    "ranges" => cur_ranges,
    "defaults" => cur_defaults,
    "plots" => plot_json
  )

  if n == 1
    JSON.print(json_dict)
  else
    file_name = args[1]
    endswith(file_name, ".json") || ( file_name *= ".json" )

    open(file_name,"w") do cur_file
      JSON.print(cur_file, json_dict)
    end
  end

  return :($(json_dict))
end

function compile_tree_dict(cur_block, cur_ranges, cur_dict=OrderedDict())
  if isempty(cur_ranges)
    return cur_block
  end

  tmp_ranges = deepcopy(cur_ranges)
  cur_symbol, cur_range = popfirst!(tmp_ranges)

  sub_dict = OrderedDict()

  for cur_value in cur_range
    tmp_block = deepcopy(cur_block)

    if isa(cur_value, Symbol)
      work_block = quote
        $(Symbol(cur_symbol)) = Symbol($(cur_value))
      end
    else
      work_block = :( $(Symbol(cur_symbol)) = $(cur_value) )
    end

    insert!(tmp_block.args, 2, work_block)
    sub_dict[cur_value] = compile_tree_dict(tmp_block, tmp_ranges, cur_dict)
  end

  return sub_dict
end

mutable struct TreeLeaf
  key::Union{Symbol,AbstractString,Number}
  value::Any
end

mutable struct TreeNode
  key::Union{Symbol,AbstractString,Number}
  children::Vector{Union{TreeNode, TreeLeaf}}
end

mutable struct TreeList
  root::Union{TreeNode, Nothing}
  nodes::Vector{Union{TreeNode, TreeLeaf}}
end

function build_node!(tree_list::TreeList, cur_key::Union{Symbol,AbstractString,Number}, nestedDict::AbstractDict)
  cur_children = []
  for (cur_key, cur_value) in nestedDict
    tmp_node = build_node!(tree_list, cur_key, cur_value)
    push!(cur_children, tmp_node)
  end

  cur_node = TreeNode(cur_key, cur_children)
  push!(tree_list.nodes, cur_node)
  return cur_node
end

function build_node!(tree_list::TreeList, cur_key::Union{Symbol,AbstractString,Number}, cur_value)
  cur_node = TreeLeaf(cur_key, cur_value)
  push!(tree_list.nodes, cur_node)
  return cur_node
end

function build_tree(nested_dict::AbstractDict)
  tree_list = TreeList(nothing, [])
  tree_list.root = build_node!(tree_list, :root, nested_dict)

  return tree_list
end

function tree_to_dict(cur_node::TreeLeaf)
  return cur_node.value
end

function tree_to_dict(cur_node::TreeNode)
  cur_dict = OrderedDict()
  for sub_node in cur_node.children
    cur_dict[sub_node.key] = tree_to_dict(sub_node)
  end
  return cur_dict
end

return


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    DirectionPartitioner(direction; tol=1e-6)

A method for partitioning spatial objects along a given `direction`
with bandwidth tolerance `tol`.
"""
struct DirectionPartitioner{T,N} <: AbstractSpatialFunctionPartitioner
  direction::SVector{N,T}
  tol::Float64
end

DirectionPartitioner(direction::NTuple{N,T}; tol=1e-6) where {T,N} =
  DirectionPartitioner{T,N}(normalize(SVector(direction)), tol)

(p::DirectionPartitioner)(x, y) = begin
  d = p.direction
  @. y = x - y
  norm(y .- (y ⋅ d) .* d) < p.tol
end


text:
# This file contains code that was formerly a part of Julia.
# License is MIT: https://julialang.org/license

using KahanSummation
using Test

@testset "cumsum_kbn" begin
    v   = [1,1e100,1,-1e100]*1000
    v2  = [1,-1e100,1,1e100]*1000

    cv  = [1,1e100,1e100,2]*1000
    cv2 = [1,-1e100,-1e100,2]*1000

    @test isequal(cumsum_kbn(v), cv)
    @test isequal(cumsum_kbn(v2), cv2)

    A = [v reverse(v) v2 reverse(v2)]

    c = cumsum_kbn(A, 1)

    @test isequal(c[:,1], cv)
    @test isequal(c[:,3], cv2)
    @test isequal(c[4,:], [2.0, 2.0, 2.0, 2.0]*1000)

    c = cumsum_kbn(A, 2)

    @test isequal(c[1,:], cv2)
    @test isequal(c[3,:], cv)
    @test isequal(c[:,4], [2.0,2.0,2.0,2.0]*1000)
end

@testset "sum_kbn" begin
    @test sum_kbn([1,1e100,1,-1e100]) === 2.0
    @test sum_kbn(Float64[]) === 0.0
    @test sum_kbn(i for i=1.0:1.0:10.0) === 55.0
    @test sum_kbn(i for i=1:1:10) === 55
    @test sum_kbn([1 2 3]) === 6
    @test sum_kbn([2+im 3-im]) === 5+0im
    @test sum_kbn([1+im 2+3im]) === 3+4im
    @test sum_kbn([7 8 9]) === sum_kbn([8 9 7])
    @test sum_kbn(i for i=1:1:10) === sum_kbn(i for i=10:-1:1)
    @test sum_kbn([-0.0]) === -0.0
    @test sum_kbn([-0.0,-0.0]) === -0.0
end


text:
#=
This example classifies CIFAR 10(https://www.cs.toronto.edu/~kriz/cifar.html)
dataset using a convolutional neural network with batch normalization

For details of applying batch normalization, see:
"Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift"
Sergey Ioffe  & Christian Szegedy
https://arxiv.org/abs/1502.03167

The trained architecture has the form of:
conv -> bn -> relu -> pool ->
conv -> bn -> relu -> pool ->
conv -> bn -> relu ->
fc -> bn  -> relu -> output

=#

using Knet, Statistics, MLDatasets

function loaddata()
    @info("Loading CIFAR 10...")
    xtrn, ytrn = CIFAR10.traindata()
    xtst, ytst = CIFAR10.testdata()
    xtrn = convert(Knet.array_type[], xtrn)
    xtst = convert(Knet.array_type[], xtst)
    #= Subtract mean of each feature
    where each channel is considered as
    a single feature following the CNN
    convention=#
    mn = mean(xtrn, dims=(1,2,4))
    xtrn = xtrn .- mn
    xtst = xtst .- mn
    @info("Loaded CIFAR 10")
    return (xtrn, ytrn), (xtst, ytst)
end

##Model definition

#=
Initialization is from
He et al., 2015,
Delving Deep into Rectifiers: Surpassing Human-Level Performance on ImageNet Classification
https://arxiv.org/abs/1502.01852
=#
kaiming(et, h, w, i, o) = et(sqrt(2 / (w * h * o))) .* randn(et, h, w, i, o)

function init_model(;et=Float32)
    # Use bnparams() to initialize gammas and betas
    w = Any[
        kaiming(et, 3, 3, 3, 16),    bnparams(et, 16),
        kaiming(et, 3, 3, 16, 32),   bnparams(et, 32),
        kaiming(et, 3, 3, 32, 64),   bnparams(et, 64),
        xavier(et, 100, 8 * 8 * 64), bnparams(et, 100),
        xavier(et, 10, 100),         zeros(et, 10, 1)
    ]
    # Initialize a moments object for each batchnorm
    m = Any[bnmoments() for i = 1:4]
    w = map(Knet.array_type[], w)
    return w, m
end

function conv_layer(w, m, x; maxpool=true)
    o = conv4(w[1], x; padding=1)
    o = batchnorm(o, m, w[2])
    o = relu.(o)
    if maxpool; o=pool(o); end
    return o
end

function lin_layer(w, m, x)
    o = w[1] * x
    o = batchnorm(o, m, w[2])
    return relu.(o)
end

function predict(w, m, x)
    o = conv_layer(w[1:2] , m[1], x)
    o = conv_layer(w[3:4] , m[2], o)
    o = conv_layer(w[5:6] , m[3], o; maxpool=false)
    o = lin_layer( w[7:8] , m[4], mat(o))
    return w[9] * o .+ w[10]
end

function loss(w, m, x, classes)
    ypred = predict(w,m, x)
    return nll(ypred, classes)
end

lossgrad = grad(loss)

# Training
function epoch!(w, m, o, xtrn, ytrn;  mbatch=64)
    data = minibatch(xtrn, ytrn, mbatch;
                   shuffle=true,
                   xtype=Knet.array_type[])
    for (x, y) in data
        g = lossgrad(w, m, x, y)
        update!(w, g, o)
    end
end

# Accuracy computation
function acc(w, m, xtst, ytst; mbatch=64)
    data = minibatch(xtst, ytst, mbatch;
                     partial=true,
                     xtype=Knet.array_type[])
    accuracy(x->predict(w,m,x), data)
end

# TODO: add command line options
function train(;optim=Momentum, epochs=5,
               lr=0.1, oparams...)
    w, m = init_model()
    o = map(_->Momentum(;lr=lr, oparams...), w)
    (xtrn, ytrn), (xtst, ytst) = loaddata()
    for epoch = 1:epochs
        println("epoch: ", epoch)
        @time epoch!(w, m, o, xtrn, ytrn)
        println("train accuracy: ", acc(w, m, xtrn, ytrn))
        println("test accuracy: ", acc(w, m, xtst, ytst))
    end
end

endswith(string(PROGRAM_FILE), "cnn_batchnorm.jl") && train()

nothing


text:
"Information container about the cells to run in a reactive call and any cells that will err."
struct CellTopology
	"Cells that form a directed acyclic graph, in topological order."
	runnable::Array{Cell, 1}
	"Cells that are in a directed cycle, with corresponding `ReactivityError`s."
	errable::Dict{Cell, ReactivityError}
end

"Return a `CellTopology` that lists the cells to be evaluated in a single reactive run, in topological order. Includes the given roots."
function topological_order(notebook::Notebook, roots::Array{Cell, 1}; allow_multiple_defs=false)::CellTopology
	entries = Cell[]
	exits = Cell[]
	errable = Dict{Cell, ReactivityError}()

	function dfs(cell::Cell)
		if cell in exits
			return
		elseif haskey(errable, cell)
			return
		elseif length(entries) > 0 && entries[end] == cell
			return # a cell referencing itself is legal
		elseif cell in entries
			currently_in = setdiff(entries, exits)
			cycle = currently_in[findfirst(isequal(cell), currently_in):end]
			for cell in cycle
				errable[cell] = CyclicReferenceError(cycle...)
			end
			return
		end

		push!(entries, cell)
		assigners = where_assigned(notebook, cell.symstate.assignments)
		if !allow_multiple_defs && length(assigners) > 1
			for c in assigners
				errable[c] = MultipleDefinitionsError(c, assigners)
			end
		end
		referencers = where_referenced(notebook, cell.symstate.assignments) |> Iterators.reverse
		for c in (allow_multiple_defs ? referencers : union(assigners, referencers))
			if c != cell
				dfs(c)
			end
		end
		push!(exits, cell)
	end

	# we first move cells to the front if they call `import` or `using`
    # we use MergeSort because it is a stable sort: leaves cells in order if they are in the same category
    prelim_order_1 = sort(roots, alg=MergeSort, by=(c -> isempty(c.module_usings)))
	# reversing because our search returns reversed order
	prelim_order_2 = Iterators.reverse(prelim_order_1)
	dfs.(prelim_order_2)
	ordered = reverse(exits)
	CellTopology(setdiff(ordered, keys(errable)), errable)
end

function disjoint(a::Set, b::Set)
	!any(x in a for x in b)
end

"Return the cells that reference any of the given symbols. Recurses down functions calls, but not down cells."
function where_referenced(notebook::Notebook, symbols::Set{Symbol})::Array{Cell, 1}
	return filter(notebook.cells) do cell
		if !disjoint(symbols, cell.symstate.references)
			return true
		end
        for func in cell.symstate.funccalls
            if haskey(notebook.combined_funcdefs, func)
                if !disjoint(symbols, notebook.combined_funcdefs[func].references)
                    return true
                end
            end
		end
		return false
	end
end

"Return the cells that assign to any of the given symbols. Recurses down functions calls, but not down cells."
function where_assigned(notebook::Notebook, symbols::Set{Symbol})::Array{Cell, 1}
	return filter(notebook.cells) do cell
		if !disjoint(symbols, cell.symstate.assignments)
			return true
		end
        for func in cell.symstate.funccalls
            if haskey(notebook.combined_funcdefs, func)
                if !disjoint(symbols, notebook.combined_funcdefs[func].assignments)
                    return true
                end
            end
		end
		return false
	end
end

"Update the combined collection of function definitions, where multiple specialisations of a function are combined into a single `SymbolsState`."
function update_funcdefs!(notebook::Notebook)
	notebook.combined_funcdefs = union((c.symstate.funcdefs for c in notebook.cells)...)
end

"Return all functions called by a cell, and all functions called by those functions, et cetera."
function all_indirect_calls(notebook::Notebook, symstate::SymbolsState, found::Set{Vector{Symbol}}=Set{Vector{Symbol}}())::Set{Vector{Symbol}}
	for func in symstate.funccalls
		if func in found
			# done
		else
			push!(found, func)
            if haskey(notebook.combined_funcdefs, func)
                inner_symstate = notebook.combined_funcdefs[func]
                all_indirect_calls(notebook, inner_symstate, found)
            end
		end
	end

	return found
end


text:
using Test, LinearAlgebra, OrderedCollections
using WaveProp
using WaveProp.Nystrom
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Mesh

@testset "Single scattering" begin
    # construct interior solution
    pde  = Helmholtz(dim=2,k=1)
    xout = SVector(3,3)
    u    = (x)   -> SingleLayerKernel(pde)(xout,x)
    dudn = (x,n) -> DoubleLayerKernel(pde)(xout,x,n)
    Geometry.clear!()
    geo  = Geometry.Kite()
    Ω    = Domain(geo)
    M    = meshgen(Ω)
    Γ    = boundary(Ω)
    # generate a Nystrom mesh with trapezoidal quadrature
    etype2qrule = OrderedDict(E => TrapezoidalP(100) for E in etypes(M))
    mesh = NystromMesh(M,Γ,etype2qrule)
    γ₀u  = γ₀(u,mesh)
    γ₁u  = γ₁(dudn,mesh)
    S    = SingleLayerOperator(pde,mesh)
    D    = DoubleLayerOperator(pde,mesh)
    Smk  = Nystrom.assemble_mk(S)
    Dmk  = Nystrom.assemble_mk(D)
    e0   = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
    e1   = WaveProp.Nystrom.error_interior_green_identity(Smk,Dmk,γ₀u,γ₁u)
    @test norm(e0,Inf) > 1e-5
    @test norm(e1,Inf) < 1e-5
end

@testset "Multiple scattering" begin
    # construct interior solution
    pde  = Helmholtz(dim=2,k=1)
    xout = SVector(-10,0)
    u    = (x)   -> SingleLayerKernel(pde)(xout,x)
    dudn = (x,n) -> DoubleLayerKernel(pde)(xout,x,n)
    geo1  = Geometry.Kite()
    geo2  = Circle(center=(10,0))
    Ω   = Domain([geo1,geo2])
    M  = meshgen(Ω,h=100)
    Γ    = boundary(Ω)
    # generate a Nystrom mesh with trapezoidal quadrature
    mesh  = NystromMesh(M,Γ;quad_rule=TrapezoidalP(100))
    γ₀u   = γ₀(u,mesh)
    γ₁u   = γ₁(dudn,mesh)
    S     = SingleLayerOperator(pde,mesh)
    D     = DoubleLayerOperator(pde,mesh)
    Smk   = Nystrom.assemble_mk(S)
    Dmk   = Nystrom.assemble_mk(D)
    e0    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
    e1 = WaveProp.Nystrom.error_interior_green_identity(Smk,Dmk,γ₀u,γ₁u)
    @test norm(e0,Inf) > 1e-5
    @test norm(e1,Inf) < 1e-5
end

@testset "Change of variables" begin
    # construct interior solution
    pde  = Helmholtz(dim=2,k=1)
    xout = SVector(-10,0)
    u    = (x)   -> SingleLayerKernel(pde)(xout,x)
    dudn = (x,n) -> DoubleLayerKernel(pde)(xout,x,n)
    geo  = ParametricEntity(ReferenceLine()) do u
        SVector(cos(2π*u[1]),sin(2π*u[1]))
    end
    bdy  = ClosedEntity([geo])
    Ω    = Domain(bdy)
    M    = meshgen(Ω)
    Γ    = boundary(Ω)
    # generate a Nystrom mesh with trapezoidal quadrature
    mesh  = NystromMesh(M,Γ;quad_rule=TrapezoidalP(40))
    γ₀u   = γ₀(u,mesh)
    γ₁u   = γ₁(dudn,mesh)
    S     = SingleLayerOperator(pde,mesh)
    D     = DoubleLayerOperator(pde,mesh)
    Smk   = Nystrom.assemble_mk(S)
    Dmk   = Nystrom.assemble_mk(D)
    e0    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
    e1    = WaveProp.Nystrom.error_interior_green_identity(Smk,Dmk,γ₀u,γ₁u)
    @test norm(e0,Inf) > 1e-5
    @test norm(e1,Inf) < 1e-5
    norm(e0,Inf)
    norm(e1,Inf)
end

@testset "Composite surface" begin
    # construct interior solution
    cov  = KressP()
    qstd = TrapezoidalP(10)
    qrule = SingularQuadratureRule(qstd,cov)
    pde  = Helmholtz(dim=2,k=1)
    xout = SVector(-2,0)
    u    = (x)   -> SingleLayerKernel(pde)(xout,x)
    dudn = (x,n) -> DoubleLayerKernel(pde)(xout,x,n)
    bdy  = Circle()
    Ω    = Domain(bdy)
    M  = meshgen(Ω;h=0.05)
    Γ    = boundary(Ω)
    # generate a Nystrom mesh with trapezoidal quadrature
    mesh  = NystromMesh(M,Γ;quad_rule=qrule)
    γ₀u   = γ₀(u,mesh)
    γ₁u   = γ₁(dudn,mesh)
    S     = SingleLayerOperator(pde,mesh)
    D     = DoubleLayerOperator(pde,mesh)
    Smk   = Nystrom.assemble_mk(S)
    Dmk   = Nystrom.assemble_mk(D)
    e0    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
    e1    = WaveProp.Nystrom.error_interior_green_identity(Smk,Dmk,γ₀u,γ₁u)
    @test norm(e0,Inf) > 1e-5
    @test norm(e1,Inf) < 1e-5
    norm(e0,Inf)
    norm(e1,Inf)
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
## This file contains utilities for scraping logs from buildbot
export BuildbotDatasource

struct BuildbotDatasource
    # Base URL this buildbot instance's API endpoints are located at.
    # Typically https://build.julialang.org/api/v2
    api_base::String

    # Mapping of builder name to ID, filled in by `load_builders!()`
    builders::Dict{String,Int}

    # Default constructor method to initialize with defaults
    # Automatically loads its builders upon construction, they
    # can be refreshed by a future `load_builders!()` invocation.
    function BuildbotDatasource(api_base::String = "https://build.julialang.org/api/v2")
        ds = new(api_base, Dict{String,Int}())
        load_builders!(ds)
        return ds
    end
end

"""
    load_builders!(ds::BuildbotDatasource)

Load the mapping of builder name (e.g. `package_linux64`) to its builder ID (an integer).
This should be fairly stable, but can change if the buildbot configuration gets changed.
This will be automatically called when a new `BuildbotDatasource` is constructed, but you
can manually refresh it by calling `load_builders!(ds)` directly.
"""
function load_builders!(ds::BuildbotDatasource)
    builders_json = get_json("$(ds.api_base)/builders")

    # We want to grab all builders that are a `package_*` or `tester_*` builder
    builders = filter(builders_json[:builders]) do builder
        return startswith(builder.name, "package_") ||
               startswith(builder.name, "tester_")
    end

    # Empty any previous `builders` that existed within the `BuildbotDatasource`
    empty!(ds.builders)

    # Load in our newly-gotten list of builders
    for builder in builders
        ds.builders[builder.name] = builder.builderid
    end
    return ds
end

"""
    get_builder_builds_list(ds::BuildbotDatasource, builder_name::String)
    get_builder_builds_list(ds::BuildbotDatasource, builder_id::Int)

Given a builder name or ID, return a sorted list of all completed builds.
"""
function get_builder_builds_list(ds::BuildbotDatasource, builder_name::String)
    return get_builder_builds_list(ds, ds.builders[builder_name])
end
function get_builder_builds_list(ds::BuildbotDatasource, builder_id::Int64)
    builders_json = get_json("$(ds.api_base)/builders/$(builder_id)/builds")
    # Only look at complete builds
    return sort(Int64[b.buildid for b in builders_json.builds if b.complete])
end

"""
    get_all_builder_builds(ds::BuildbotDatasource)

Get list of all build IDs that are stored within this `BuildbotDatasource`.
"""
function get_all_builder_builds(ds::BuildbotDatasource)
    builds = Dict{String,Vector{Int64}}()
    for (builder_name, builder_id) in ds.builders
        builds[builder_name] = get_builder_builds_list(ds, builder_id)
    end
    return builds
end

"""
    get_build_steps_list(ds::BuildbotDatasource, build_id::Int64)

Given a build ID, return the list of steps that we should download logs for.
"""
function get_build_steps_list(ds::BuildbotDatasource, build_id::Int64)
    steps_json = get_json("$(ds.api_base)/builds/$(build_id)/steps")
    return sort(Int64[s.number for s in steps_json.steps if s.complete])
end

function get_log_contents(ds::BuildbotDatasource, build_id::Int64, step::Int64)
    logs_json = get_json("$(ds.api_base)/builds/$(build_id)/steps/$(step)/logs")
    log_ids = Int64[l.logid for l in logs_json.logs if l.complete]

    content = ""
    for log_id in log_ids
        contents_json = get_json("$(ds.api_base)/logs/$(log_id)/contents")
        for chunk in contents_json.logchunks
            # Buildbot tags all lines with what stream it came from;
            #  - i for stdin
            #  - o for stdout
            #  - e for stderr
            #  - h for header (buildbot-generated textual content)
            # For now, we're just going to strip all these out.
            content = string(content, join([l[2:end] for l in split(chunk.content, "\n")], "\n"))
        end
    end
    return content
end

function download_logs(ds::BuildbotDatasource,
                       output_dir::String;
                       tail::Int64 = 100,
                       verbose::Bool = false)
    # We will store logs partitioned by `$builder_name/$build/$step.log`
    if verbose
        @info("Enumerating builders...")
    end
    builder_builds = get_all_builder_builds(ds)
    if verbose
        @info("Found $(length(keys(builder_builds))) builders")
    end

    for builder_name in keys(builder_builds)
        # Limit to only the last `tail` logs for each builder
        for build_id in builder_builds[builder_name][max(end-tail+1, 1):end]
            build_logdir = joinpath(output_dir, builder_name, string(build_id))
            mkpath(build_logdir)
            for step_id in get_build_steps_list(ds, build_id)
                log_filename = joinpath(build_logdir, "$(step_id).log")
                if isfile(log_filename)
                    continue
                end
                if verbose
                    @info("Downloading #$(build_id).$(step_id)")
                end

                data = get_log_contents(ds, build_id, step_id)
                open(log_filename, "w") do io
                    write(io, data)
                end
            end
        end
    end
end


text:
using POMDPModels
using POMDPs
using POMDPModelTools
using StaticArrays
using Statistics
using Random
using DiscreteValueIteration
using GlobalApproximationValueIteration
using Flux
using Test

Random.seed!(1234)

@testset "all" begin

    @testset "integration" begin
        include("test_with_linear_gfa.jl")
        include("test_with_nonlinear_gfa.jl")
    end
end


text:
using Manopt, Manifolds, Test, Random

@testset "Hessian access functions" begin
    M = Euclidean(2)
    F(M, p) = 1
    gradF(M, p) = zero(2)
    gradF!(M, X, p) = copyto!(M, X, p, zeros(2))
    HessF(M, p, X) = X
    HessF!(M, Y, p, X) = copyto!(M, Y, p, X)
    precon = (M, p, X) -> X
    P1 = HessianProblem(M, F, gradF, HessF, precon)
    P2 = HessianProblem(M, F, gradF!, HessF!, precon; evaluation=MutatingEvaluation())

    p = zeros(2)
    X = ones(2)
    @test get_hessian(P1, p, X) == get_hessian(P2, p, X)
    Y1 = similar(X)
    Y2 = similar(X)
    get_hessian!(P1, Y1, p, X)
    get_hessian!(P2, Y2, p, X)
    @test Y1 == Y2
end


text: noodds(values::Array{Int,1})::Array{Int,1} = filter(iseven,values)

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# Attributes

export Attribute
"""
    struct Attribute

Holds a C pointer `adios2_attribute *`.
"""
struct Attribute
    ptr::Ptr{Cvoid}
    adios::Adios
    Attribute(ptr::Ptr{Cvoid}, adios::Adios) = new(ptr, adios)
end

function Base.show(io::IO, attribute::Attribute)
    nm = name(attribute)
    T = type(attribute)
    isval = is_value(attribute)
    sz = size(attribute)
    d = data(attribute)
    print(io, "Attribute(name=$nm,type=$T,is_value=$isval,size=$sz,data=$d)")
    return
end

function Base.show(io::IO, ::MIME"text/plain", attribute::Attribute)
    nm = name(attribute)
    return print(io, "Attribute($nm)")
end

export name
"""
    attr_name = name(attribute::Attribute)
    attr_name::Union{Nothing,String}

Retrieve attribute name.
"""
function name(attribute::Attribute)
    size = Ref{Csize_t}()
    err = ccall((:adios2_attribute_name, libadios2_c), Cint,
                (Ptr{Cchar}, Ref{Csize_t}, Ptr{Cvoid}), C_NULL, size,
                attribute.ptr)
    Error(err) ≠ error_none && return nothing
    name = Array{Cchar}(undef, size[])
    err = ccall((:adios2_attribute_name, libadios2_c), Cint,
                (Ptr{Cchar}, Ref{Csize_t}, Ptr{Cvoid}), name, size,
                attribute.ptr)
    Error(err) ≠ error_none && return nothing
    return unsafe_string(pointer(name), size[])
end

export type
"""
    attr_type = type(attribute::Attribute)
    attr_type::Union{Nothing,Type}

Retrieve attribute type.
"""
function type(attribute::Attribute)
    type = Ref{Cint}()
    err = ccall((:adios2_attribute_type, libadios2_c), Cint,
                (Ref{Cint}, Ptr{Cvoid}), type, attribute.ptr)
    Error(err) ≠ error_none && return nothing
    return julia_type(AType(type[]))
end

export is_value
"""
    attr_is_value = is_value(attribute::Attribute)
    attr_is_value::Union{Nothing,Bool}

Retrieve attribute type.
"""
function is_value(attribute::Attribute)
    is_value = Ref{Cint}()
    err = ccall((:adios2_attribute_is_value, libadios2_c), Cint,
                (Ref{Cint}, Ptr{Cvoid}), is_value, attribute.ptr)
    Error(err) ≠ error_none && return nothing
    return Bool(is_value[])
end

"""
    attr_size = size(attribute::Attribute)
    attr_size::Union{Nothing,Int}

Retrieve attribute size.
"""
function Base.size(attribute::Attribute)
    size = Ref{Csize_t}()
    err = ccall((:adios2_attribute_size, libadios2_c), Cint,
                (Ref{Csize_t}, Ptr{Cvoid}), size, attribute.ptr)
    Error(err) ≠ error_none && return nothing
    return Int(size[])
end

export data
"""
    attr_data = data(attribute::Attribute)
    attr_data::Union{Nothing,AdiosType,Vector{<:AdiosType}}

Retrieve attribute Data.
"""
function data(attribute::Attribute)
    T = type(attribute)
    T ≡ nothing && return nothing
    @assert T != Union{}
    isval = is_value(attribute)
    isval ≡ nothing && return nothing
    sz = size(attribute)
    sz ≡ nothing && return nothing
    tp = type(attribute)
    tp ≡ nothing && return nothing
    if tp ≡ String
        if isval
            buffer = fill(Cchar(0), string_array_element_max_size)
            out_sz = Ref{Csize_t}()
            err = ccall((:adios2_attribute_data, libadios2_c), Cint,
                        (Ptr{Cchar}, Ref{Csize_t}, Ptr{Cvoid}), buffer, out_sz,
                        attribute.ptr)
            @assert out_sz[] == sz
            Error(err) ≠ error_none && return nothing
            data = unsafe_string(pointer(buffer))
            return data
        else
            arrays = [Array{Cchar}(undef, string_array_element_max_size)
                      for i in 1:sz]
            buffers = [pointer(array) for array in arrays]::Vector{Ptr{Cchar}}
            out_sz = Ref{Csize_t}()
            err = ccall((:adios2_attribute_data, libadios2_c), Cint,
                        (Ptr{Ptr{Cchar}}, Ref{Csize_t}, Ptr{Cvoid}), buffers,
                        out_sz, attribute.ptr)
            @assert out_sz[] == sz
            Error(err) ≠ error_none && return nothing
            data = unsafe_string.(buffers)::Vector{String}
            # Use `arrays` again to ensure it is not GCed too early
            buffers .= pointer.(arrays)
            return data
        end
    else
        data = Array{T}(undef, sz)
        out_sz = Ref{Csize_t}()
        err = ccall((:adios2_attribute_data, libadios2_c), Cint,
                    (Ptr{Cvoid}, Ref{Csize_t}, Ptr{Cvoid}), data, out_sz,
                    attribute.ptr)
        @assert out_sz[] == sz
        Error(err) ≠ error_none && return nothing
        isval && return data[1]
        return data
    end
end


text:
module ClimaSimulations

using Oceananigans.Utils:
    prettytime,
    TimeInterval,
    IterationInterval,
    WallTimeInterval

using Oceananigans.Simulations:
    Simulation,
    run!,
    time_step!,
    initialize_simulation!,
    evaluate_callbacks!,
    evaluate_output_writers!,
    @stopwatch,
    Callback,
    iteration_limit_exceeded,
    stop_time_exceeded,
    wall_time_limit_exceeded

end # module


text:
using Test
using DataFrames
using Deldir
using Gadfly
using GeometricalPredicates

include("../src/triangle.jl")

x = [0.2, 0.3, 0.4, 0.4, 0.5, 0.6]
y = [0.0, 0.2, 0.0, 0.4, 0.2, 0.0]

del, vor, sum = deldir(x, y)

function plotexample()
    x, y = edges(del)
    println(del)
    println("-----------------------------")
    println(sum)
    Gadfly.plot(x = x, y = y, Geom.path, Coord.cartesian(fixed = true))
end

@testset "Proper generators are generated" begin
    xsize = size(x)
    ysize = size(y)
    d, v, s = deldir(x, y)
    generators = emptygenerators(s)
    collectedgenerators = collect(generators)
    @test xsize == ysize
    @test size(keys(collectedgenerators)) == xsize
    nonemptygenerators = filter(g -> g == [], values(collectedgenerators))
    @test size(nonemptygenerators, 1) == 0
end

@testset "Generators are filled properly" begin
    expectedfilledgeneratos =
        Dict([
            (1, [2, 3]),
            (2, [1, 3, 4, 5]),
            (3, [1, 2, 5, 6]),
            (4, [2, 5]),
            (5, [4, 2, 6, 3]),
            (6, [3, 5])
        ])
    d, v, s = deldir(x, y)
    generators = emptygenerators(s)
    generatorstogenerated!(generators, d)
    collectedgenerators = collect(generators)
    cexpectedgenerators = collect(expectedfilledgeneratos)
    @test size(keys(collectedgenerators)) == size(keys(cexpectedgenerators))
    @test keys(collectedgenerators) == keys(cexpectedgenerators)
    i = 1
    while i <= size(keys(collectedgenerators), 1)
        sortedgenerators = sort(generators[i])
        sortedexpectedgenerators = sort(expectedfilledgeneratos[i])
        @test sortedgenerators == sortedexpectedgenerators
        i += 1
    end
 end

 @testset "Generators are transformed to triangles build with indexes" begin
    expectedtriangles = [
        [1, 2, 3],
        [2, 4, 5],
        [3, 5, 6],
        [2, 4, 5]
    ]
    d, v, s = deldir(x, y)
    generators = emptygenerators(s)
    generatorstogenerated!(generators, d)
    triangles = buildindextriangles(generators)
    @test size(triangles) == size(expectedtriangles)
    for expectedtriangle in expectedtriangles
        @test expectedtriangle in triangles
    end
end

@testset "Index triangles are transofmed to value triangles" begin
    expectedtriangles = [
        [Point2D(0.2, 0.0), Point2D(0.3, 0.2), Point2D(0.4, 0.0)],
        [Point2D(0.3, 0.2), Point2D(0.4, 0.4), Point2D(0.5, 0.2)],
        [Point2D(0.4, 0.0), Point2D(0.5, 0.2), Point2D(0.6, 0.0)],
        [Point2D(0.3, 0.2), Point2D(0.4, 0.4), Point2D(0.5, 0.2)]
    ]
    d, v, s = deldir(x, y)
    generators = emptygenerators(s)
    generatorstogenerated!(generators, d)
    triangles = buildindextriangles(generators)
    triangles = buildvaluetriangles(triangles, s)
    println(typeof(triangles))
    @test size(triangles) == size(expectedtriangles)
    for expectedtriangle in expectedtriangles
        @test expectedtriangle in triangles
    end
end

@testset "Build triangles" begin
    expectedtriangles = [
        [Point2D(0.2, 0.0), Point2D(0.3, 0.2), Point2D(0.4, 0.0)],
        [Point2D(0.3, 0.2), Point2D(0.4, 0.4), Point2D(0.5, 0.2)],
        [Point2D(0.4, 0.0), Point2D(0.5, 0.2), Point2D(0.6, 0.0)],
        [Point2D(0.3, 0.2), Point2D(0.4, 0.4), Point2D(0.5, 0.2)]
    ]
    d, v, s = deldir(x, y)
    tr = triangles(d, s)
    @test size(tr) == size(expectedtriangles)
    for expectedtriangle in expectedtriangles
        @test expectedtriangle in  tr
    end
end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
mutable struct PipelineLayoutCreateInfo
    mLayouts::Vector{vk.VkDescriptorSetLayout}
    mRanges::Vector{vk.VkPushConstantRange}

    mHandleRef::Ref{vk.VkPipelineLayoutCreateInfo}

    function PipelineLayoutCreateInfo()
        this = new()
        this.mLayouts = Vector{vk.VkDescriptorSetLayout}()
        this.mRanges = Vector{vk.VkPushConstantRange}()
        return this
    end
end

function addSetLayout(this::PipelineLayoutCreateInfo, layout::vk.VkDescriptorSetLayout)
    push!(this.mLayouts, layout)
end

function addPushConstantRange(this::PipelineLayoutCreateInfo, range::vk.VkPushConstantRange)
    push!(this.mRanges, range)
end

function handleRef(this::PipelineLayoutCreateInfo)::Ref{vk.VkPipelineLayoutCreateInfo}
    this.mHandleRef = Ref(vk.VkPipelineLayoutCreateInfo(
        vk.VK_STRUCTURE_TYPE_PIPELINE_LAYOUT_CREATE_INFO, #sType::VkStructureType
        C_NULL, #pNext::Ptr{Cvoid}
        0, #flags::VkPipelineLayoutCreateFlags
        length(this.mLayouts), #setLayoutCount::UInt32
        pointer(this.mLayouts), #pSetLayouts::Ptr{VkDescriptorSetLayout}
        length(this.mRanges), #pushConstantRangeCount::UInt32
        pointer(this.mRanges) #pPushConstantRanges::Ptr{VkPushConstantRange}
    ))
    return this.mHandleRef
end

function getRanges(this::PipelineLayoutCreateInfo)::Vector{vk.VkPushConstantRange}
    return this.mRanges
end


text:
using ConvexTests, SDPAFamily
using MathOptInterface, Dualization
const MOI = MathOptInterface
using GenericLinearAlgebra

@info "Starting SDPA-GMP tests"

presolve_exclusions = [r"lp_max_atom", r"lp_min_atom", r"affine_Partial_transpose"]
do_tests("SDPA_GMP", () -> SDPAFamily.Optimizer(presolve=true, silent=true); exclude = presolve_exclusions ∪ [r"mip", r"exp", r"sdp_Complex_Semidefinite_constraint"], variant="presolve=true", T = BigFloat, description="Tests with SDPA-GMP via SDPAFamily.jl. The problems `lp_max_atom`, `lp_min_atom` and `affine_Partial_transpose` are excluded due to excessive presolve time (> 15 minutes each).", last = false, first = true)

@info "Starting presolve=false tests"

do_tests("SDPA_GMP", () -> SDPAFamily.Optimizer(presolve=false, silent=true); exclude = [r"mip", r"exp", r"sdp_Complex_Semidefinite_constraint"], variant="presolve=false", first = false, T = BigFloat, last = false)

# Dualization does not work with BigFloats yet: https://github.com/JuliaOpt/Dualization.jl/issues/82.
# function dual_opt(;kwargs...)
#     () -> DualOptimizer{BigFloat}(SDPAFamily.Optimizer(; kwargs...))
# end

# @info "Starting dual tests"

# do_tests("SDPA_GMP", dual_opt(presolve=true, silent=true); exclude = presolve_exclusions ∪ [r"mip", r"exp", r"sdp_Complex_Semidefinite_constraint"], variant="presolve=true (dualized)", T = BigFloat, description="Tests with SDPA-GMP via SDPAFamily.jl. The problems `lp_max_atom` and `lp_min_atom` are excluded due to excessive presolve time (~1000 seconds each).", first = false, last = false)

# do_tests("SDPA_GMP", dual_opt(presolve=false, silent=true); exclude =[r"mip", r"exp", r"sdp_Complex_Semidefinite_constraint"], variant="presolve=false (dualized)", T = BigFloat, description="Tests with SDPA-GMP via SDPAFamily.jl.", first = false, last = true)