text:
using LinearAlgebra

# Constantes
COMANDO_START = "start"

# Funcoes do terminal
function clearterminal()
    run(`clear`)
end

function print_setup_label()
    clearterminal()
    println("#################################### SETUP ####################################")
end

function print_label()
    clearterminal()
    println("############################ Análise de Estruturas ############################")
end

function print_resumo_projeto()
    println("Projeto: $(projeto)")
    println("Tipo de Estrutura: $(tipoestrutura)")
    println()
end

# variavel etapa controla em qual etapa do setup o usuário está
etapa = 1
# variavel comando controla qual comando atual o usuário digitou
comando = COMANDO_START

while comando != "fim"
    print_setup_label()

    if comando == COMANDO_START
        print("Nome do projeto: ")
        global projeto = readline()
        global comando = "tipo_estrutura"
    elseif comando == "tipo_estrutura"
        print("Tipo de Estrutura [Portico, Grelha, Viga, pLaca]: ")
        global tipoestrutura = readline()
        if startswith(lowercase(tipoestrutura), "p")
            global tipoestrutura = "PORTICO"
        elseif startswith(lowercase(tipoestrutura), "g")
            global tipoestrutura = "GRELHA"
        elseif startswith(lowercase(tipoestrutura), "v")
            global tipoestrutura = "VIGA"
        elseif startswith(lowercase(tipoestrutura), "l")
            global tipoestrutura = "PLACA"
        else
            global tipoestrutura = "INVALIDO"
        end
        global comando = "fim"
    else
        global comando = "fim"
    end
end

print_label()
print_resumo_projeto()

max_stars_repo_name: rsd-consultorias/scripts-engenharia

max_issues_repo_name: rsd-consultorias/scripts-engenharia

max_forks_repo_name: rsd-consultorias/scripts-engenharia


text:
function tss(y::AbstractVector, hasintercept::Bool, sqrtw::AbstractVector)
    if hasintercept
        m = (mean(y) / sum(sqrtw) * length(y))::Float64
        return sum(abs2(y[i] - sqrtw[i] * m) for i in eachindex(y))
    else
        return sum(abs2, y)
    end
end

function Fstat(coef::Vector{Float64}, matrix_vcov::AbstractMatrix{Float64}, has_intercept::Bool)
    coefF = copy(coef)
    # TODO: check I can't do better
    length(coef) == has_intercept && return NaN
    if has_intercept
        coefF = coefF[2:end]
        matrix_vcov = matrix_vcov[2:end, 2:end]
    end
    return (coefF' * (matrix_vcov \ coefF)) / length(coefF)
end

##############################################################################
##
## VcovData stores data you need to compute errors
##
##############################################################################
struct VcovData{T, N} <: RegressionModel
    modelmatrix::Matrix{Float64}       # X
    crossmodelmatrix::T                    # X'X in the simplest case. Can be Matrix but preferably Factorization
    residuals::Array{Float64, N}      # vector or matrix of residuals (matrix in the case of IV, residuals of Xendo on (Z, Xexo))
    dof_residual::Int
end
StatsBase.modelmatrix(x::VcovData) = x.modelmatrix
Vcov.crossmodelmatrix(x::VcovData) = x.crossmodelmatrix
StatsBase.residuals(x::VcovData) = x.residuals
StatsBase.dof_residual(x::VcovData) = x.dof_residual

ext: jl


text:
module HSVGP

using LinearAlgebra, Flux, Zygote, Distances, Statistics, Random, LatinHypercubeSampling, NearestNeighbors
using MPI, CSV, DataFrames, FileIO, Distributions, StatsBase
using SpecialFunctions: loggamma
using ProgressBars

include("structs.jl")
include("prediction.jl")
include("likelihoods.jl")
include("elbos.jl")
include("fit_functions.jl")
include("covariance_function.jl")
include("mpi_comm.jl")
include("fit_mpi_hsvgp.jl")
include("fit_mpi_psvgp.jl")

end

max_stars_repo_name: ayanbiswas/PRISM

max_issues_repo_name: ayanbiswas/PRISM

max_forks_repo_name: ayanbiswas/PRISM


text:
# That's not great but it will do for now,
# It would be better to reuse Base.REPL code

abstract type REPLMode end

struct NormalMode <: REPLMode end
struct ShellMode <: REPLMode end
struct HelpMode <: REPLMode end
struct PkgMode <: REPLMode end

global const REPLModes = [NormalMode(), ShellMode(), HelpMode(), PkgMode()]

prompt(c, mode::REPLMode) = string(mode, ">")
prompt(c, mode::NormalMode) = string(c.eval_in, ">")
prompt(c, mode::ShellMode) = "shell>"
prompt(c, mode::HelpMode) = "help?>"
prompt(c, mode::PkgMode) = "pkg>"

switch_key(mode::NormalMode) = GdkKeySyms.BackSpace#Action(GdkKeySyms.BackSpace, NoModifier, "")
switch_key(mode::ShellMode)  = 0x03b#Action(0x03b, NoModifier, "")
switch_key(mode::PkgMode)    = 0x05d#Action(0x05d, NoModifier, "")
switch_key(mode::HelpMode)   = 0x03f#Action(0x03f, GdkModifierType.SHIFT, "")
#FIXME add these to GdkKeySyms https://gitlab.gnome.org/GNOME/gtk/blob/master/gdk/gdkkeysyms.h

on_return(c, mode::REPLMode, cmd)    = nothing
on_return(c, mode::NormalMode, cmd)  = remotecall_fetch(RemoteGtkREPL.eval_command_remotely, worker(c), cmd, string(c.eval_in))
on_return(c, mode::HelpMode, cmd)    = remotecall_fetch(RemoteGtkREPL.eval_command_remotely, worker(c), "@doc $cmd", string(c.eval_in))
on_return(c, mode::ShellMode, cmd)   = remotecall_fetch(RemoteGtkREPL.eval_shell_remotely,   worker(c), cmd, string(c.eval_in))
on_return(c, mode::PkgMode, cmd)     = remotecall_fetch(RemoteGtkREPL.eval_command_remotely, worker(c), "Pkg.REPLMode.pkgstr(\"$cmd\")", "Main")


text:
using ClimateMachine.DGMethods.NumericalFluxes: NumericalFluxFirstOrder
import ClimateMachine.DGMethods.NumericalFluxes:
    EntropyConservative,
    numerical_volume_conservative_flux_first_order!,
    numerical_volume_fluctuation_flux_first_order!,
    ave,
    logave,
    numerical_flux_first_order!,
    numerical_flux_second_order!,
    numerical_boundary_flux_second_order!
import ClimateMachine.BalanceLaws: entropy_variables_to_state!, state_to_entropy, wavespeed

@inline gamma(ps::EarthParameterSet) = cp_d(ps) / cv_d(ps)

function numerical_volume_conservative_flux_first_order!(
    ::EntropyConservative,
    ::DryAtmosModel,
    F::Grad,
    state_1::Vars,
    aux_1::Vars,
    state_2::Vars,
    aux_2::Vars,
)
    FT = eltype(F)
    ρ_1, ρu_1, ρe_1 = state_1.ρ, state_1.ρu, state_1.ρe
    ρ_2, ρu_2, ρe_2 = state_2.ρ, state_2.ρu, state_2.ρe
    Φ_1, Φ_2 = aux_1.Φ, aux_2.Φ
    u_1 = ρu_1 / ρ_1
    u_2 = ρu_2 / ρ_2
    p_1 = pressure(ρ_1, ρu_1, ρe_1, Φ_1)
    p_2 = pressure(ρ_2, ρu_2, ρe_2, Φ_2)
    b_1 = ρ_1 / 2p_1
    b_2 = ρ_2 / 2p_2

    ρ_avg = ave(ρ_1, ρ_2)
    u_avg = ave(u_1, u_2)
    b_avg = ave(b_1, b_2)
    Φ_avg = ave(Φ_1, Φ_2)

    usq_avg = ave(dot(u_1, u_1), dot(u_2, u_2))

    ρ_log = logave(ρ_1, ρ_2)
    b_log = logave(b_1, b_2)
    α = b_avg * ρ_log / 2b_1

    γ = FT(gamma(param_set))

    Fρ = u_avg * ρ_log
    Fρu = u_avg * Fρ' + ρ_avg / 2b_avg * I
    if total_energy
        Fρe = (1 / (2 * (γ - 1) * b_log) - usq_avg / 2 + Φ_avg) * Fρ + Fρu * u_avg
    else
        Fρe = (1 / (2 * (γ - 1) * b_log) - usq_avg / 2) * Fρ + Fρu * u_avg
    end

    F.ρ += Fρ
    F.ρu += Fρu
    F.ρe += Fρe
end

function numerical_volume_fluctuation_flux_first_order!(
    ::NumericalFluxFirstOrder,
    ::DryAtmosModel,
    D::Grad,
    state_1::Vars,
    aux_1::Vars,
    state_2::Vars,
    aux_2::Vars,
)
    if fluctuation_gravity
        FT = eltype(D)
        ρ_1, ρu_1, ρe_1 = state_1.ρ, state_1.ρu, state_1.ρe
        ρ_2, ρu_2, ρe_2 = state_2.ρ, state_2.ρu, state_2.ρe
        Φ_1, Φ_2 = aux_1.Φ, aux_2.Φ
        p_1 = pressure(ρ_1, ρu_1, ρe_1, Φ_1)
        p_2 = pressure(ρ_2, ρu_2, ρe_2, Φ_2)
        b_1 = ρ_1 / 2p_1
        b_2 = ρ_2 / 2p_2

        ρ_log = logave(ρ_1, ρ_2)
        b_avg = ave(b_1, b_2)
        α = b_avg * ρ_log / 2b_1

        D.ρu -= α * (Φ_1 - Φ_2) * I
    end
end

struct CentralVolumeFlux <: NumericalFluxFirstOrder end
function numerical_volume_conservative_flux_first_order!(
    ::CentralVolumeFlux,
    m::DryAtmosModel,
    F::Grad,
    state_1::Vars,
    aux_1::Vars,
    state_2::Vars,
    aux_2::Vars,
)
    FT = eltype(F)
    F_1 = similar(F)
    flux_first_order!(m, F_1, state_1, aux_1, FT(0), EveryDirection())

    F_2 = similar(F)
    flux_first_order!(m, F_2, state_2, aux_2, FT(0), EveryDirection())

    parent(F) .= (parent(F_1) .+ parent(F_2)) ./ 2
end

struct KGVolumeFlux <: NumericalFluxFirstOrder end
function numerical_volume_conservative_flux_first_order!(
    ::KGVolumeFlux,
    ::DryAtmosModel,
    F::Grad,
    state_1::Vars,
    aux_1::Vars,
    state_2::Vars,
    aux_2::Vars,
)

    Φ_1 = aux_1.Φ
    ρ_1 = state_1.ρ
    ρu_1 = state_1.ρu
    ρe_1 = state_1.ρe
    u_1 = ρu_1 / ρ_1
    e_1 = ρe_1 / ρ_1
    p_1 = pressure(ρ_1, ρu_1, ρe_1, Φ_1)

    Φ_2 = aux_2.Φ
    ρ_2 = state_2.ρ
    ρu_2 = state_2.ρu
    ρe_2 = state_2.ρe
    u_2 = ρu_2 / ρ_2
    e_2 = ρe_2 / ρ_2
    p_2 = pressure(ρ_2, ρu_2, ρe_2, Φ_2)

    ρ_avg = ave(ρ_1, ρ_2)
    u_avg = ave(u_1, u_2)
    e_avg = ave(e_1, e_2)
    p_avg = ave(p_1, p_2)

    F.ρ = ρ_avg * u_avg
    F.ρu = p_avg * I + ρ_avg * u_avg .* u_avg'
    F.ρe = ρ_avg * u_avg * e_avg + p_avg * u_avg
end

struct EntropyConservativeWithPenalty <: NumericalFluxFirstOrder end
function numerical_flux_first_order!(
    numerical_flux::EntropyConservativeWithPenalty,
    balance_law::BalanceLaw,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}

    FT = eltype(fluxᵀn)

    numerical_flux_first_order!(
        EntropyConservative(),
        balance_law,
        fluxᵀn,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )
    fluxᵀn = parent(fluxᵀn)

    wavespeed⁻ = wavespeed(balance_law, normal_vector, state_prognostic⁻, state_auxiliary⁻, t, direction)
    wavespeed⁺ = wavespeed(balance_law, normal_vector, state_prognostic⁺, state_auxiliary⁺, t, direction)
    max_wavespeed = max.(wavespeed⁻, wavespeed⁺)
    penalty = max_wavespeed .* (parent(state_prognostic⁻) - parent(state_prognostic⁺))

    fluxᵀn .+= penalty / 2
end

Base.@kwdef struct MatrixFlux{FT} <: NumericalFluxFirstOrder
    Mcut::FT = 0
    low_mach::Bool = false
    kinetic_energy_preserving::Bool = false
end

function numerical_flux_first_order!(
    numerical_flux::MatrixFlux,
    balance_law::BalanceLaw,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}

    FT = eltype(fluxᵀn)
    numerical_flux_first_order!(
        EntropyConservative(),
        balance_law,
        fluxᵀn,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )
    fluxᵀn = parent(fluxᵀn)

    γ = FT(gamma(param_set))

    low_mach = numerical_flux.low_mach
    Mcut = numerical_flux.Mcut
    kinetic_energy_preserving = numerical_flux.kinetic_energy_preserving

    ω = FT(π) / 3
    δ = FT(π) / 5
    random_unit_vector = SVector(sin(ω) * cos(δ), cos(ω) * cos(δ), sin(δ))
    # tangent space basis
    τ1 = random_unit_vector × normal_vector
    τ2 = τ1 × normal_vector

    ρ⁻ = state_prognostic⁻.ρ
    ρu⁻ = state_prognostic⁻.ρu
    ρe⁻ = state_prognostic⁻.ρe

    Φ⁻ = state_auxiliary⁻.Φ
    u⁻ = ρu⁻ / ρ⁻
    p⁻ = pressure(ρ⁻, ρu⁻, ρe⁻, Φ⁻)
    β⁻ = ρ⁻ / 2p⁻

    Φ⁺ = state_auxiliary⁺.Φ
    ρ⁺ = state_prognostic⁺.ρ
    ρu⁺ = state_prognostic⁺.ρu
    ρe⁺ = state_prognostic⁺.ρe

    u⁺ = ρu⁺ / ρ⁺
    p⁺ = pressure(ρ⁺, ρu⁺, ρe⁺, Φ⁺)
    β⁺ = ρ⁺ / 2p⁺

    ρ_log = logave(ρ⁻, ρ⁺)
    β_log = logave(β⁻, β⁺)
    if total_energy
        Φ_avg = ave(Φ⁻, Φ⁺)
    else
        Φ_avg = 0
    end
    u_avg = ave.(u⁻, u⁺)
    p_avg = ave(ρ⁻, ρ⁺) / 2 * ave(β⁻, β⁺)
    u²_bar = 2 * sum(u_avg .^ 2) - sum(ave(u⁻ .^ 2, u⁺ .^ 2))

    h_bar = γ / (2 * β_log * (γ - 1)) + u²_bar / 2 + Φ_avg
    c_bar = sqrt(γ * p_avg / ρ_log)

    umc = u_avg - c_bar * normal_vector
    upc = u_avg + c_bar * normal_vector
    u_avgᵀn = u_avg' * normal_vector
    R = hcat(
        SVector(1, umc[1], umc[2], umc[3], h_bar - c_bar * u_avgᵀn),
        SVector(1, u_avg[1], u_avg[2], u_avg[3], u²_bar / 2 + Φ_avg),
        SVector(0, τ1[1], τ1[2], τ1[3], τ1' * u_avg),
        SVector(0, τ2[1], τ2[2], τ2[3], τ2' * u_avg),
        SVector(1, upc[1], upc[2], upc[3], h_bar + c_bar * u_avgᵀn),
    )

    if low_mach
        M = abs(u_avg' * normal_vector) / c_bar
        c_bar *= max(min(M, FT(1)), Mcut)
    end

    if kinetic_energy_preserving
        λl = abs(u_avgᵀn) + c_bar
        λr = λl
    else
        λl = abs(u_avgᵀn - c_bar)
        λr = abs(u_avgᵀn + c_bar)
    end

    Λ = SDiagonal(λl, abs(u_avgᵀn), abs(u_avgᵀn), abs(u_avgᵀn), λr)
    #Ξ = sqrt(abs((p⁺ - p⁻) / (p⁺ + p⁻)))
    #Λ = Ξ * abs(u_avgᵀn + c_bar) * I + (1 - Ξ) * ΛM

    T = SDiagonal(ρ_log / 2γ, ρ_log * (γ - 1) / γ, p_avg, p_avg, ρ_log / 2γ)

    entropy⁻ = similar(parent(state_prognostic⁻), Size(6))
    state_to_entropy_variables!(
        balance_law,
        Vars{vars_state(balance_law, Entropy(), FT)}(entropy⁻),
        state_prognostic⁻,
        state_auxiliary⁻,
    )

    entropy⁺ = similar(parent(state_prognostic⁺), Size(6))
    state_to_entropy_variables!(
        balance_law,
        Vars{vars_state(balance_law, Entropy(), FT)}(entropy⁺),
        state_prognostic⁺,
        state_auxiliary⁺,
    )

    Δentropy = parent(entropy⁺) - parent(entropy⁻)

    fluxᵀn .-= R * Λ * T * R' * Δentropy[SOneTo(5)] / 2
end

"""
    state_to_entropy_variables!(
        ::DryAtmosModel,
        entropy::Vars,
        state::Vars,
        aux::Vars,
    )

See [`BalanceLaws.state_to_entropy_variables!`](@ref)
"""
function state_to_entropy_variables!(::DryAtmosModel, entropy::Vars, state::Vars, aux::Vars)
    ρ, ρu, ρe, Φ = state.ρ, state.ρu, state.ρe, aux.Φ

    FT = eltype(state)
    γ = FT(gamma(param_set))

    p = pressure(ρ, ρu, ρe, Φ)
    s = log(p / ρ^γ)
    b = ρ / 2p
    u = ρu / ρ

    if total_energy
        entropy.ρ = (γ - s) / (γ - 1) - (dot(u, u) - 2Φ) * b
    else
        entropy.ρ = (γ - s) / (γ - 1) - (dot(u, u)) * b
    end
    entropy.ρu = 2b * u
    entropy.ρe = -2b
    entropy.Φ = 2ρ * b
end

"""
    entropy_variables_to_state!(
        ::DryAtmosModel,
        state::Vars,
        aux::Vars,
        entropy::Vars,
    )

See [`BalanceLaws.entropy_variables_to_state!`](@ref)
"""
function entropy_variables_to_state!(::DryAtmosModel, state::Vars, aux::Vars, entropy::Vars)
    FT = eltype(state)
    β = entropy
    γ = FT(gamma(param_set))

    b = -β.ρe / 2
    ρ = β.Φ / (2b)
    ρu = ρ * β.ρu / (2b)

    p = ρ / (2b)
    s = log(p / ρ^γ)
    Φ = dot(ρu, ρu) / (2 * ρ^2) - ((γ - s) / (γ - 1) - β.ρ) / (2b)

    ρe = p / (γ - 1) + dot(ρu, ρu) / (2ρ) + ρ * Φ

    state.ρ = ρ
    state.ρu = ρu
    state.ρe = ρe
    aux.Φ = Φ
end

function state_to_entropy(::DryAtmosModel, state::Vars, aux::Vars)
    FT = eltype(state)
    ρ, ρu, ρe, Φ = state.ρ, state.ρu, state.ρe, aux.Φ
    p = pressure(ρ, ρu, ρe, Φ)
    γ = FT(gamma(param_set))
    s = log(p / ρ^γ)
    η = -ρ * s / (γ - 1)
    return η
end

# Numerical Flux
#=
numerical_flux_first_order!(::Nothing, _...) = nothing
numerical_boundary_flux_first_order!(::Nothing, _...) = nothing
numerical_flux_second_order!(::Nothing, _...) = nothing
numerical_boundary_flux_second_order!(::Nothing, _...) = nothing
=#
numerical_flux_first_order!(::Nothing, ::DryAtmosModel, _...) = nothing
numerical_flux_second_order!(::Nothing, ::DryAtmosModel, _...) = nothing
numerical_boundary_flux_second_order!(::Nothing, a, ::DryAtmosModel, _...) = nothing

function wavespeed(lm::DryAtmosLinearModel, nM, state::Vars, aux::Vars, t::Real, direction)
    ref = aux.ref_state
    return soundspeed(ref.ρ, ref.p)
end

function wavespeed(::DryAtmosModel, nM, state::Vars, aux::Vars, t::Real, direction)
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe
    Φ = aux.Φ
    p = pressure(ρ, ρu, ρe, Φ)

    u = ρu / ρ
    uN = abs(dot(nM, u))
    return uN + soundspeed(ρ, p)
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
struct HSumSpace{T} <: HRepresentation{T}
    d::Int
end
hvectortype(::Type{HSumSpace{T}}) where T = SparseVector{T, Int}
Base.length(idxs::HyperPlaneIndices{T, HSumSpace{T}}) where T = idxs.rep.d
Base.isempty(idxs::HyperPlaneIndices{T, HSumSpace{T}}) where T = idxs.rep.d <= 0
startindex(idxs::HyperPlaneIndices{T, HSumSpace{T}}) where T = isempty(idxs) ? nothing : eltype(idxs)(1)
function Base.get(h::HSumSpace{T}, idx::HyperPlaneIndex) where T
    i = idx.value
    a = sparsevec(i:h.d:i+2h.d, StaticArrays.SVector(one(T), -one(T), -one(T)), 3h.d)
    return HyperPlane(a, zero(T))
end
function nextindex(h::HSumSpace, idx::HyperPlaneIndex)
    if idx.value < h.d
        return typeof(idx)(idx.value + 1)
    else
        return nothing
    end
end
@norepelem HSumSpace HalfSpace

struct Projection{S, I}
    set::S
    dimensions::I
end
fulldim(p::Projection) = length(p.dimensions)
dimension_names(p::Projection) = dimension_names(p.set)[p.dimensions]

function __pad_map(i, d, n, start)
    if i <= start
        return i
    elseif i <= start + n
        return nothing
    else
        return i - n
    end
end
function __pad_map(i, d, n)
    if n < 0
        if i <= n
            return nothing
        else
            return i - n
        end
    else
        if i <= d
            return i
        else
            return nothing
        end
    end
end
function _pad_map(i, d, padding...)
    j = __pad_map(i, d, padding...)
    return j === nothing ? nothing : (1, j)
end

# TODO should be cartesian product with FullSpace
function zeropad(p::HRep{T}, padding...) where T
    d = map_fulldim(N -> N + abs(padding[1]), FullDim(p))
    f = (i, el) -> zeropad(el, padding...)
    return similar(p, d, T, hmap(f, d, T, p)...;
                   dimension_map = i -> _pad_map(i, fulldim(p), padding...))
end

function Base.intersect(p1::Projection{S1, <:Base.OneTo}, p2::Projection{S2, <:Base.OneTo}) where {S1, S2}
    d = length(p1.dimensions)
    d == length(p2.dimensions) || throw(DimensionMismatch())
    p = intersect(zeropad(p1.set, fulldim(p2.set) - d),
                  zeropad(p2.set, fulldim(p1.set) - d, d))
    return Projection(p, Base.OneTo(d))
end

"""
    convexhull(p1::HRepresentation, p2::HRepresentation)

Returns the Balas [Theorem 3.3, B85] extended H-representation of the convex hull
of `p1` and `p2`.

[B85] Balas, E., 1985.
*Disjunctive programming and a hierarchy of relaxations for discrete optimization problems*.
SIAM Journal on Algebraic Discrete Methods, 6(3), pp.466-486.
"""
function convexhull(p1::HRepresentation, p2::HRepresentation)
    d = FullDim(p1)
    T = promote_coefficient_type((p1, p2))
    d_2 = map_fulldim(N -> 2N, d)
    d_3 = map_fulldim(N -> 3N, d)
    function f(i, el)
        if i == 3
            return zeropad(el, 1)
        elseif i == 4
            return zeropad(el, -d_3)
        else
            if i == 1
                _padded = zeropad(el, neg_fulldim(d))
                padded = zeropad(_padded, d)
                return lift(padded, false)
            else
                padded = zeropad(el, neg_fulldim(d_2))
                return complement_lift(padded, false)
            end
        end
    end
    ei = basis(hvectortype(p1), map_fulldim(N -> 1, d), 1)
    λ = HalfSpace(-ei, zero(T)) ∩ HalfSpace(ei, one(T))
    lifted = similar((p1, p2), map_fulldim(N -> 3N + 1, d), T,
        hmap(f, d, T, p1, p2, HSumSpace{T}(fulldim(p1)), λ)...;
        dimension_map = i -> (1 <= i <= fulldim(d)) ? (1, i) : nothing)
    return Projection(lifted, Base.OneTo(fulldim(p1)))
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


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

"""
    intersecttype(l1, l2)

Compute the intersection type of two lines `l1` and `l2`.

The intersection type can be one of three types:

1. intersect at one point
2. overlap at more than one point
3. do not overlap nor intersect
"""
function intersecttype(l1::Line{3,T}, l2::Line{3,T}) where {T}
  a, b = l1(0), l1(1)
  c, d = l2(0), l2(1)

  if measure(Tetrahedron(a, b, c, d)) > 0
    NoIntersection()
  elseif isapprox(abs((b - a) × (c - d)), zero(T), atol=atol(T)^2)
    OverlappingLines(l1)
  else
    CrossingLines(intersectpoint(l1, l2))
  end
end

function intersecttype(l1::Line{2,T}, l2::Line{2,T}) where {T}
  a, b = l1(0), l1(1)
  c, d = l2(0), l2(1)

  if !isapprox(abs((b - a) × (c - d)), zero(T), atol=atol(T)^2)
    CrossingLines(intersectpoint(l1, l2))
  elseif isapprox(measure(Triangle(a, b, c)), zero(T), atol=atol(T)^2)
    OverlappingLines(l1)
  else
    NoIntersection()
  end
end

# compute the intersection of two lines assuming that it is a point
function intersectpoint(l1::Line{2}, l2::Line{2})
  a, b = l1(0), l1(1)
  c, d = l2(0), l2(1)

  v1  = a - b
  v2  = c - d
  v12 = a - c

  # the intersection point lies in between a and b at a fraction s
  # (https://en.wikipedia.org/wiki/Line-line_intersection#Formulas)
  s = (v12[1] * v2[2] - v12[2] * v2[1]) / (v1[1] * v2[2] - v1[2] * v2[1])

  a - s*v1
end


text:
using QMTK
using QMTK.Consts.Pauli
using Compat.Test

@testset "Region Iterator" begin

    chain = Chain(Fixed, 5)
    t = kronsum(bonds(chain, 1)) do i, j
        :(σ₁[$i] ⊗ σ₂[$j])
    end

    for i = 1:3
        rhs = rand(Bit, 5)
        h = LocalHamiltonian(Region{1}, σ₁ ⊗ σ₂)
        itr = RegionIterator(h, chain, rhs)

        rhs_idx = convert(Int, rhs) + 1
        for (val, lhs) in RegionIterator(h, chain, rhs)
            lhs_idx = convert(Int, lhs) + 1
            @test val == t[lhs_idx, rhs_idx]
        end
    end

    ###

    square = Square(Fixed, 4, 4)
    t = kronsum(bonds(square, 1)) do i, j
        :(σ₁[$i] ⊗ σ₂[$j])
    end

    for i=1:3
        rhs = rand(Bit, 4, 4)
        h = LocalHamiltonian(Region{1}, σ₁ ⊗ σ₂)
        itr = RegionIterator(h, square, rhs)

        rhs_idx = convert(Int, rhs) + 1
        for (val, lhs) in RegionIterator(h, square, rhs)
            lhs_idx = convert(Int, lhs) + 1
            @test val == t[lhs_idx, rhs_idx]
        end
    end

    chain = Chain(Periodic, 5)
    t = kronsum(bonds(chain, 1)) do i, j
        :(σ₁[$i] ⊗ σ₂[$j])
    end

    for i = 1:3
        rhs = rand(Bit, 5)
        h = LocalHamiltonian(Region{1}, σ₁ ⊗ σ₂)
        itr = RegionIterator(h, chain, rhs)

        rhs_idx = convert(Int, rhs) + 1
        for (val, lhs) in RegionIterator(h, chain, rhs)
            lhs_idx = convert(Int, lhs) + 1
            @test val == t[lhs_idx, rhs_idx]
        end
    end

    ###

    square = Square(Periodic, 4, 4)
    t = kronsum(bonds(square, 1)) do i, j
        :(σ₁[$i] ⊗ σ₂[$j])
    end

    for i=1:3
        rhs = rand(Bit, 4, 4)
        h = LocalHamiltonian(Region{1}, σ₁ ⊗ σ₂)
        itr = RegionIterator(h, square, rhs)

        rhs_idx = convert(Int, rhs) + 1
        for (val, lhs) in RegionIterator(h, square, rhs)
            lhs_idx = convert(Int, lhs) + 1
            @test val == t[lhs_idx, rhs_idx]
        end
    end

end


text:
# module HWM14
# isapple = Sys.isapple()
# islinux = Sys.islinux()

# export hwm_dir
hwm_dir  = @__DIR__
const _hwm14_so_filename =
            string(
                    @__DIR__, "/hwm14",
                    isapple ? "_macos" : islinux ? "_linux" : "_win",
                    ".so",
                )

_check_hwm14() = _check_so(_hwm14_so_filename)
_check_hwm14()

# export hwm14
# need cd to the directory including data files
# this_pwd = pwd()
# need_cd_hwm = this_pwd != hwm_dir
# if need_cd_hwm cd(hwm_dir)   end
# hwm14(iyd, sec, alt, glat, glon,
#         hrl_32, fbar, f10p7, aaps_hwm,
#         W, )
# if need_cd_hwm cd(this_pwd) end
function _init_hwm()
    this_pwd = pwd()
    if this_pwd != @__DIR__
        # @warn "work directory should change to $(@__DIR__) before call hwm"
        cd(@__DIR__)
        @warn "work directory (hwm14) have changed\n from: $this_pwd\n   to: $(@__DIR__)"
        # c = true
    end
    this_pwd
end

hwm14(IYD::Int32, SEC::Float32,
    ALT::Float32, GLAT::Float32, GLONG::Float32,
    STL::Float32, F107A::Float32, F107::Float32,
    AP::Array{Float32,1},
    W::Array{Float32,1},
) = ccall((:hwm14_, _hwm14_so_filename),
    Nothing,
    (Ref{Int32},   # IYD
        Ref{Float32}, # SEC
        Ref{Float32}, # ALT
        Ref{Float32}, # GLAT
        Ref{Float32}, # GLONG
        Ref{Float32}, # STL
        Ref{Float32}, # F107A
        Ref{Float32}, # F107
        Ptr{Float32}, #  AP
        Ptr{Float32},),
    IYD, SEC,
    ALT, GLAT, GLONG,
    STL, F107A, F107,
    AP,
    W,
)

function test_hwm()
    aaps_hwm  = Array{Float32,1}(undef, 2)
    aaps_hwm .= 1.0

    W = Array{Float32,1}(undef, 2)
    iyd = Int32(19199)
    sec = Float32(12 * 3600.0)
    alt  = Float32(30.0)
    glat = Float32(15.0)
    glon = Float32(36.0)
    hrl_32 = Float32(0.0)
    fbar = Float32(180.0)
    f10p7 = Float32(180.0)

    tvn0 = 1.0

    this_pwd = _init_hwm()
    hwm14(iyd, sec, alt, glat, glon,
            hrl_32, fbar, f10p7, aaps_hwm,
            W, )

    cd(this_pwd)

    v = 100.0 * W[1] * tvn0
    u = 100.0 * W[2] * tvn0
    # print(u, "\t", v)
    u, v
    # @test isapprox(u,  -2843.706703186035)
    # @test isapprox(v, -34.654366970062256)
end

# end module HWM14
# end


text:
"""
Package: DataWrangler

normalize\\[!\\]([x]; method)

Normalize values in vector `x` either in-place or returning a copy of `x` with the imputed values

# Parameters

`x`: Vector of type Real with missing values to be imputed

`method`: There are four valid values: "z-score", "min-max", "softmax", "sigmoid"

# Description

The normalization is applied ignoring any missing values in the Array, and it takes place in all the available dimensions of the array. If normalization is required in just some specific dimensions the function `mapslices` can be used to select those dimensions.

# Example
```julia
x = [1.,2,3,4,5]
normalize!(x; method = "min-max")

println(x)
[0.0, 0.25, 0.5, 0.75, 1.0]
```
"""
function normalize!(x::AbstractArray{<:Union{Missing,T}};
                 method::String = "z-score") where T<:Real

    @assert !(T <: Integer) " Method `$(typeof(x))` cannot be normalized in-place, use `normalize(x)` instead"

    @assert method ∈ ["z-score","min-max","softmax","sigmoid"]

    smx = skipmissing(x)

    if method == "z-score"
        μ = mean(smx)
        σ = std(smx)
        for i in CartesianIndices(x)
            x[i] = (x[i] - μ) / σ
        end
    end

    if method == "min-max"
        m,M = extrema(smx)
        Mm = M-m
        for i in CartesianIndices(x)
            x[i] = (x[i] - m) / Mm
        end
    end

    if method == "softmax"
        Σe = sum(exp.(smx))
        for i in CartesianIndices(x)
            x[i] = exp(x[i]) / Σe
        end
    end

    if method == "sigmoid"
        for i in CartesianIndices(x)
            x[i] = 1/(1+exp(-x[i]))
        end
    end

end

"""
Package: DataWrangler

Check 'normalize! for further information.
"""
function normalize(x::AbstractArray{<:Union{Missing,T}};
                method::String = "z-score") where T<:Real

    P = any(ismissing,x) ? Union{Missing, Base.promote_op(/,T,T)} : Base.promote_op(/,T,T)
    cx = Array{P}(undef,size(x))
    cx[:] = x
    normalize!(cx; method)
    return cx

end


text:
using MAT
using Dates
using GLPK
using JuMP
using Plots
using Plots.PlotMeasures
include("../utils.jl")

function get_wave_t_in_range(param, min, max)
    indices_to_keep = findall(>=(min), param)
    param = param[indices_to_keep]

    indices_to_keep = findall(<(max), param)
    param = param[indices_to_keep]
    param, indices_to_keep
end

function fit_poly_inf_norm(x,y,n)
	m = Model(GLPK.Optimizer)

    @variable(m, alpha[1:n+1])
	@variable(m, t >= 0.)
	@constraint(m, [i=1:length(y)], sum(alpha[j]*x[i]^(j-1) for j in 1:n+1) - y[i] <= t);
	@constraint(m, [i=1:length(y)], -t <= sum(alpha[j]*x[i]^(j-1) for j in 1:n+1) - y[i]);
	@objective(m, Min, t)

	optimize!(m) # added new 2020-04-08

    return JuMP.value.(alpha)
end

f2(x, args) = (args[2] * x) .+ args[1]

const MATLAB_EPOCH = Dates.DateTime(-1,12,31)

date2num(d::Dates.DateTime) = Dates.value(d-MATLAB_EPOCH)/(1000*60*60*24)
num2date(n::Number) =  MATLAB_EPOCH + Dates.Millisecond(round(Int64, n*1000*60*60*24))

function plot_bathy(xs, profile; show_points=false, x_last_positive=NaN, x_first_negative=NaN, value_last_positive=NaN, value_first_negative=NaN, x_at_zero=NaN, title=NaN)
    p1 = plot(xs, profile, label="profile_t0")
    if show_points
        scatter!([x_last_positive, x_first_negative], [value_last_positive, value_first_negative], label="A, B")
        scatter!([x_at_zero], [0], label="Shoreline")
    end
    display(plot(p1, title=title))
end

function find_x_at_zero(profile, xs, i)
    t0_index_last_positive = findlast(profile.>0)
    if t0_index_last_positive == nothing
        println(i, ": t0_index_last_positive")
        t0_index_last_positive = 1
    end
    # if t0_index_last_positive == nothing t0_index_last_positive = 1 end

    t0_index_first_negative = findfirst(profile.<0)
    if t0_index_first_negative == nothing
        println(i, ": t0_index_first_negative")
        t0_index_first_negative = 1
    end

    t0_x_last_positive = xs[t0_index_last_positive]
    t0_x_first_negative = xs[t0_index_first_negative]
    t0_value_last_positive = profile[t0_index_last_positive]
    t0_value_first_negative = profile[t0_index_first_negative]

    t0_coeffs = fit_poly_inf_norm([t0_value_last_positive, t0_value_first_negative], [t0_x_last_positive, t0_x_first_negative], 1)
    t0_x_at_zero = f2([0], t0_coeffs)[1]
    t0_x_at_zero
end

dir_path = get_dataset_path(;dataset_name="grand_popo")
kalman_bathy_vars = matread(string(dir_path, "Kalman_bathy.mat"))
bathy_xb = kalman_bathy_vars["xb"][1, :]# 1×213 Matrix{Float64}
# bathy_tn = kalman_bathy_vars["tn"]      # 1×1214 Matrix{Float64}
bathy_hbK = kalman_bathy_vars["hbK"]    # 213×1214 Matrix{Float64}

min_target_matdate = minimum(filter(!isnan, kalman_bathy_vars["tn"]))
max_target_matdate = maximum(filter(!isnan, kalman_bathy_vars["tn"]))

wave_vars = matread(string(dir_path, "Wave_GPP_ERA.mat"))
wave_t = wave_vars["t"]
wave_Tm = wave_vars["Tm"]
wave_hs = wave_vars["hs"]
wave_dir = wave_vars["dir"]

wave_t, indices = get_wave_t_in_range(wave_t, min_target_matdate, max_target_matdate)
wave_Tm = wave_Tm[indices]
wave_hs = wave_hs[indices]
wave_dir = wave_dir[indices]

bathy_dates = map(num2date, collect(min_target_matdate:max_target_matdate))  # collect() instead of using bathy_tn because bathy_tn includes NaN's, even though bathy_hbK has no NaN
wave_dates = map(num2date, wave_t)

changes = []

for i in 2:length(bathy_dates)
    _d0 = bathy_dates[i-1]
    _d1 = bathy_dates[i]

    _t0 = date2num(_d0)
    _t1 = date2num(_d1)

    _bathy_t0 = bathy_hbK[:, i-1]
    _bathy_t1 = bathy_hbK[:, i]

    _wave_t = deepcopy(wave_t)
    _wave_Tm = deepcopy(wave_Tm)
    _wave_dir = deepcopy(wave_dir)
    _wave_hs = deepcopy(wave_hs)

    _wave_t, _indices = get_wave_t_in_range(_wave_t, _t0, _t1)
    _wave_Tm = _wave_Tm[_indices]
    _wave_hs = _wave_hs[_indices]
    _wave_dir = _wave_dir[_indices]

    t0_x_at_zero = find_x_at_zero(_bathy_t0, bathy_xb, i)
    t1_x_at_zero = find_x_at_zero(_bathy_t1, bathy_xb, i)
    push!(changes, t1_x_at_zero - t0_x_at_zero)

end

p1 = plot(bathy_dates[2:end], changes, ylabel="Shoreline change", legend=false)
p2 = plot(wave_dates, wave_Tm, ylabel="Tm", legend=false)
p3 = plot(wave_dates, wave_hs, ylabel="Hs", legend=false)
p4 = plot(wave_dates, wave_dir, ylabel="Dir", legend=false)

full_plot = plot(p1, p2, p3, p4, layout=(4, 1), legend=false)#, title="Grand Popo shoreline change")
plot!(size=(1920, 1080), left_margin=10mm) # 1920, 1080  # 1280, 720
savefig("/Users/mn/JuliaProjects/cgp-morphcast/output/grand_popo/wavesandshorelinechange.pdf")

max_stars_repo_name: mahmoud-al-najar/cgp-morphcast

max_issues_repo_name: mahmoud-al-najar/cgp-morphcast

max_forks_repo_name: mahmoud-al-najar/cgp-morphcast


text: println("Loading...")

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text:
"""
    Panda <: Environment

    6-degree-of-freedom robotic arm designed by Franka Emika
"""
struct Panda end

function panda(;
    representation=:minimal,
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=0.8,
    damper=10.0,
    spring=0.0,
    info=nothing,
    model_type=:end_effector,
    seed=1,
    contact=true,
    limits=true,
    vis=Visualizer(),
    name=:robot,
    # infeasible_control=false,
    opts_step=SolverOptions(rtol=3.0e-4, btol=3.0e-4, undercut=1.5),
    opts_grad=SolverOptions(rtol=3.0e-4, btol=3.0e-4, undercut=1.5),
    T=Float64)

    mechanism = get_mechanism(:panda,
        timestep=timestep,
        gravity=gravity,
        friction_coefficient=friction_coefficient,
        damper=damper,
        spring=spring,
        contact=contact,
        limits=limits,
        model_type=model_type)

    initialize!(mechanism, :panda)

    if representation == :minimal
        nx = minimal_dimension(mechanism)
    elseif representation == :maximal
        nx = maximal_dimension(mechanism)
    end
    nu = input_dimension(mechanism)
    no = nx

    aspace = BoxSpace(nu,
        low=(-1.0e8 * ones(nu)),
        high=(1.0e8 * ones(nu)))
    ospace = BoxSpace(no,
        low=(-Inf * ones(no)),
        high=(Inf * ones(no)))

    rng = MersenneTwister(seed)
    z = get_maximal_state(mechanism)
    x = representation == :minimal ? maximal_to_minimal(mechanism, z) : z
    fx = zeros(nx, nx)
    fu = zeros(nx, nu)

    u_prev = zeros(nu)
    control_mask = I(nu)
    control_scaling = Diagonal(ones(nu))

    build_robot(mechanism, vis=vis, name=name)

    TYPES = [Panda, T, typeof(mechanism), typeof(aspace), typeof(ospace), typeof(info)]
    Environment{TYPES...}(mechanism, representation, aspace, ospace,
        x, fx, fu,
        u_prev, control_mask' * control_scaling,
        nx, nu, no,
        info,
        [rng], vis,
        opts_step, opts_grad)
end


text:
#===============================================================================

    thermal_consistency.jl

    72°C thermal consistency test

===============================================================================#

import DataStructures.OrderedDict
import Memento: debug

#===============================================================================
    constants >>
===============================================================================#

## preset values
# const MIN_FLUORESCENCE_VAL = 8e5 ## not used
const MIN_TM_VAL = 77
const MAX_TM_VAL = 81
const MAX_DELTA_TM_VAL = 2

#===============================================================================
    function >>
===============================================================================#

## called by dispatch()
##
## NB default values supplied to mc_analysis()
## are the same as for melting curve experiments
function act(
    ::Type{Val{thermal_consistency}},
    ## remove MySql dependency
    # db_conn ::MySQL.MySQLHandle,
    # exp_id ::Integer,
    # stage_id ::Integer,
    # calib_info ::Union{Integer,OrderedDict};
    req                 ::Associative;
    dcv                 ::Bool          = DEFAULT_CAL_DCV, ## if true, perform multi-channel deconvolution
    auto_span_smooth    ::Bool          = DEFAULT_MC_AUTO_SPAN_SMOOTH,
    span_smooth_default ::Real          = DEFAULT_MC_SPAN_SMOOTH_DEFAULT,
    span_smooth_factor  ::Real          = DEFAULT_MC_SPAN_SMOOTH_FACTOR,
    max_temperature     ::Real          = DEFAULT_MC_MAX_TEMPERATURE, ## maximum temperature to analyze
    out_format          ::OutputFormat  = pre_json_output,
    reporting           ::Function      = roundoff(JSON_DIGITS) ## reporting function
)
    debug(logger, "at act(::Type{Val{thermal_consistency}})")
    #
    ## required fields
    @get_calibration_data_from_req(thermal_consistency)
    @parse_raw_data_from_req(thermal_consistency)
    #
    ## keyword arguments
    kwargs = MC_FIELD_DEFS |>
        sift(x -> x.key in keys(req)) |>
        mold() do x
			if isa(req[x.key], Int64) && Int_T != Int64
				x.name => Int_T(req[x.key])
			else
				x.name => req[x.key]
			end
        end
    #
    ## create container for data and parameter values
    interface = McInput(
            parsed_raw_data...,
            calibration_data;
            calibration_args = CalibrationParameters(dcv = dcv),
            auto_span_smooth = auto_span_smooth,
            span_smooth_default = span_smooth_default,
            span_smooth_factor = span_smooth_factor,
            max_temperature = max_temperature,
            out_format = full_output,
            reporting = reporting,
            kwargs...)
    #
    ## analyse data analyses melting curve
    mc_w72c =
        mc_analysis(interface)
    #
    ## process the data from only one channel
    ## PROBLEM >> this does not seem appropriate for dual channel analysis
    mc_tm = map(
        field(:peaks_filtered),
        mc_w72c.peak_output[:, CHANNELS[1]]) ## mc_matrix
    min_Tm = max_temperature + 1
    max_Tm = 0
    tm_check_vec = map(mc_tm) do Ta
        if size(Ta, 1) == 0
            TmCheck1w((NaN_T, false), NaN_T)
        else
            top1_Tm = Ta[1,1]
            (top1_Tm < min_Tm) && (min_Tm = top1_Tm)
            (top1_Tm > max_Tm) && (max_Tm = top1_Tm)
            TmCheck1w(
                (top1_Tm, MIN_TM_VAL <= top1_Tm <= MAX_TM_VAL),
                Ta[1,2])
        end ## if size
    end ## do Ta
    #
    ## return values
    delta_Tm_val = max_Tm - min_Tm
    full_out() =
        ThermalConsistencyOutput(
            tm_check_vec,
            (delta_Tm_val, delta_Tm_val .<= MAX_DELTA_TM_VAL),
            true)
    pre_json_out() =
        OrderedDict(
            :tm_check => tm_check_vec,
            :delta_Tm => (reporting(delta_Tm_val),
                            delta_Tm_val .<= MAX_DELTA_TM_VAL),
            :valid    => true)
    (out_format == full_output) && return full_out()
    ## else
    return pre_json_out() |> out(out_format)
end ## act(::Type{Val{thermal_consistency}})

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
using Flux
using OceanParameterizations
using WindMixing
using LinearAlgebra
using JLD2

BLAS.set_num_threads(32)

train_files = ["wind_-5e-4_cooling_3e-8_new"]

PATH = pwd()
OUTPUT_PATH = joinpath(PATH, "training_output")
# OUTPUT_PATH = joinpath("C:\\Users\\xinle\\Documents\\OceanParameterizations.jl")

FILE_NAME_uw = "uw_NN_training_test"
FILE_NAME_vw = "vw_NN_training_test"
FILE_NAME_wT = "wT_NN_training_test"

FILE_PATH_uw = joinpath(OUTPUT_PATH, "$(FILE_NAME_uw).jld2")
FILE_PATH_vw = joinpath(OUTPUT_PATH, "$(FILE_NAME_vw).jld2")
FILE_PATH_wT = joinpath(OUTPUT_PATH, "$(FILE_NAME_wT).jld2")

EXTRACTED_FILE_PATH = joinpath(PATH, "extracted_training_output")
EXTRACTED_FILE_PATH_uw = joinpath(EXTRACTED_FILE_PATH, "$(FILE_NAME_uw)_extracted.jld2")
EXTRACTED_FILE_PATH_vw = joinpath(EXTRACTED_FILE_PATH, "$(FILE_NAME_vw)_extracted.jld2")
EXTRACTED_FILE_PATH_wT = joinpath(EXTRACTED_FILE_PATH, "$(FILE_NAME_wT)_extracted.jld2")

train_parameters = Dict("ν₀" => 1f-4, "ν₋" => 0.1f0, "Riᶜ" => 0.25f0, "ΔRi" => 1f-1, "Pr" => 1f0, "κ" => 10f0,
                        "modified_pacanowski_philander" => true, "convective_adjustment" => false,
                        "smooth_profile" => false, "smooth_NN" => false, "smooth_Ri" => false, "train_gradient" => true,
                        "zero_weights" => true, "unscaled" => false)

𝒟train = WindMixing.data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)

N_inputs = 96
hidden_units = 400
N_outputs = 31

uw_NN = Chain(Dense(N_inputs, hidden_units, relu), Dense(hidden_units, N_outputs))
vw_NN = Chain(Dense(N_inputs, hidden_units, relu), Dense(hidden_units, N_outputs))
wT_NN = Chain(Dense(N_inputs, hidden_units, relu), Dense(hidden_units, N_outputs))

# weights, re = Flux.destructure(Chain(Dense(N_inputs, hidden_units, relu), Dense(hidden_units, N_outputs)))

# uw_NN = re(weights ./ 1f5)
# vw_NN = re(weights ./ 1f5)
# wT_NN = re(weights ./ 1f5)

# train_optimizers = [ADAM(0.01), Descent()]
# train_epochs = [50,100]

train_optimizers = [ADAM(5e-4), ADAM(2e-4)]
train_epochs = [200, 200]

write_metadata_NN_training(FILE_PATH_uw, train_files, train_parameters, train_epochs, train_optimizers, uw_NN, "uw")
uw_weights = train_NN(uw_NN, 𝒟train, train_optimizers, train_epochs, FILE_PATH_uw, "uw",
                      modified_pacanowski_philander=train_parameters["modified_pacanowski_philander"],
                      convective_adjustment=train_parameters["convective_adjustment"],
                      ν₀=train_parameters["ν₀"], ν₋=train_parameters["ν₋"], ΔRi=train_parameters["ΔRi"], Riᶜ=train_parameters["Riᶜ"],
                      κ=train_parameters["κ"],
                      smooth_profile=train_parameters["smooth_profile"], smooth_NN=train_parameters["smooth_NN"], smooth_Ri=train_parameters["smooth_Ri"],
                      train_gradient=train_parameters["train_gradient"],
                      zero_weights = train_parameters["zero_weights"])

extract_NN(FILE_PATH_uw, EXTRACTED_FILE_PATH_uw, "NN")

write_metadata_NN_training(FILE_PATH_vw, train_files, train_parameters, train_epochs, train_optimizers, vw_NN, "vw")
vw_weights = train_NN(vw_NN, 𝒟train.uvT_scaled, 𝒟train.vw.scaled, train_optimizers, train_epochs, FILE_PATH_vw, "vw")

extract_NN(FILE_PATH_vw, EXTRACTED_FILE_PATH_vw, "NN")

write_metadata_NN_training(FILE_PATH_wT, train_files, train_parameters, train_epochs, train_optimizers, wT_NN, "wT")
wT_weights = train_NN(wT_NN, 𝒟train.uvT_scaled, 𝒟train.wT.scaled, train_optimizers, train_epochs, FILE_PATH_wT, "wT")

extract_NN(FILE_PATH_wT, EXTRACTED_FILE_PATH_wT, "NN")


text:
module SimpleNL

import Base: ==, >=, <=, getindex, setindex!, add_sum, copyto!, ndims, size, copyto!, convert, string, print, show, log2, max, csch, acoth, tanh, asecd, ^, *, cospi, csc, abs, hypot, /, rem2pi, -, acos, min, exp10, tan, acsc, log10, log1p, \, acotd, sin, -, exp2, acot, acosh, atan, mod, cotd, atanh, sind, rem, mod2pi, secd, atan, sinh, +, cot, cosh, acsch, sinpi, atand, cos, acosd, rad2deg, +, tand, cscd, asinh, cbrt, asin, transpose, inv, acscd, expm1, log, deg2rad, asec, asind, cosd, abs2, sech, sqrt, asech, sec, exp, coth
import SpecialFunctions: erfi, bessely, besselj, loggamma, erfcinv, hankelh2, hankelh1, erfcx, besselk, beta, invdigamma, bessely1, besselj1, dawson, airyaiprime, erf, digamma, gamma, airyai, airybi, erfinv, bessely0, erfc, trigamma, besseli, polygamma, logbeta, airybiprime, besselj0
import Requires: @require

export Constant, Variable, Parameter, Field, Gradient, Jacobian, Hessian, SparseHessian, SparseJacobian, SparseGradient, variable, parameter, objective, constraint,optimize!,instantiate!, num_variables, num_constraints, value, dual, objective_value, setvalue, set_lower_bound, set_upper_bound, lower_bound, upper_bound, set_optimizer, function_evaluator, gradient_evaluator, sparse_gradient_evaluator, hessian_evaluator, sparse_hessian_evaluator, field_evaluator, jacobian_evaluator, sparse_jacobian_evaluator, nlp_evaluator

abstract type Expression end
abstract type Gradient end
abstract type Hessian end
abstract type Field end
abstract type Entry end

mutable struct MyRef{T}
    x::T
end
@inline getindex(ref::MyRef{T}) where T = ref.x
@inline setindex!(ref::MyRef{T},val) where T = ref.x = val

include("functionlist.jl")
include("expression.jl")
include("gradient.jl")
include("hessian.jl")
include("field.jl")
include("sourcesink.jl")
include("jacobian.jl")
include("lagrangianhessian.jl")
include("sparse.jl")
include("rules.jl")
include("template.jl")
include("utils.jl")
include("evaluator.jl")
include("model.jl")
include("print.jl")

function __init__()
    @require Ipopt="b6b21f68-93f8-5de0-b562-5493be1d77c9" @eval begin
        include("Interfaces/ipopt.jl")
    end
    @require MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6" @eval begin
        include("Interfaces/madnlp.jl")
    end
end

end


text:
# ps4000enumerateunits.jl

export ps4000enumerateunits

function ps4000enumerateunits()
  count = Ref{UInt16}()
  serials = Array(UInt8, 32)
  serials[1] = 1  # make sure not null
  seriallength = Ref{UInt16}(32)
  ps_status = ccall((:ps4000EnumerateUnits, ps4000driver),
          PicoStatus,
          (Ref{UInt16},Ptr{UInt8},Ref{UInt16}),
          count,serials,seriallength)
  checkps_status(ps_status)
  serials[end] = 0
  return (count[], serials, seriallength[])
end


text:
# Note: This test suite has many hard coded values in it.
# While that is normally a bad idea, we need to guarantee reproducible results
# when using the same seed.

global SEED = 0 # Do not change this
struct MyAlphabet <: Alphabet end

@testset "Random LongSequences" begin
function test_sampler(sampler, seed, elements, firstten, T)
    rng = StableRNG(seed)
    sampled1 = [rand(rng, sampler) for i in 1:1000]
    sampled2 = rand(StableRNG(seed), sampler, 1000)
    @test sampled1 == sampled2
    @test Set(sampled1) == Set(T[convert(T, i) for i in elements])
    @test eltype(sampler) == T
    @test eltype(firstten) == T
    @test rand(StableRNG(seed), sampler, 10) == firstten
end

function test_isseq(seq, alphabettype, len)
    @test seq isa LongSequence{alphabettype}
    @test length(seq) == len
end

@testset "SamplerUniform" begin
    # Cannot instantiate with empty collection
    @test_throws ArgumentError sampler = SamplerUniform{DNA}(DNA[])

    # DNA sampler with DNA array
    sampler = SamplerUniform{DNA}([DNA_A, DNA_C, DNA_G, DNA_W])
    firstten = DNA[DNA_C, DNA_W, DNA_G, DNA_C, DNA_C, DNA_C, DNA_C, DNA_A, DNA_A, DNA_A]
    test_sampler(sampler, SEED, sampler.elems, firstten, DNA)
    # Now RNA sampler from DNA array
    sampler = SamplerUniform{RNA}([DNA_A, DNA_C, DNA_G, DNA_W])
    firstten = RNA[RNA_C, RNA_W, RNA_G, RNA_C, RNA_C, RNA_C, RNA_C, RNA_A, RNA_A, RNA_A]
    test_sampler(sampler, SEED, sampler.elems, firstten, RNA)

    # Cannot make AA sampler from DNA
    @test_throws MethodError s = SamplerUniform{AminoAcid}([DNA_A, DNA_C, DNA_G, DNA_W])

    # Automatically infer eltype
    sampler1 = SamplerUniform{DNA}([DNA_A, DNA_C])
    sampler2 = SamplerUniform([DNA_A, DNA_C])
    @test typeof(sampler2) == SamplerUniform{DNA}
    @test eltype(sampler1) == eltype(sampler2)

    # Can also infer abstract eltype
    sampler3 = SamplerUniform([DNA_A, RNA_A])
    @test eltype(sampler3) == typejoin(DNA, RNA)

    sampler = SamplerUniform{RNA}([DNA_A, DNA_C, DNA_G, DNA_W])
    @test typeof(rand(sampler)) == RNA

end # SamplerUniform

@testset "SamplerWeighted" begin
    # Must have one less weight than elements
    @test_throws ArgumentError SamplerWeighted{DNA}([DNA_A], [0.5])
    @test_throws ArgumentError SamplerWeighted{DNA}([DNA_A, DNA_C], [0.5, 0.5])
    @test_throws ArgumentError SamplerWeighted{DNA}([DNA_A], [0.5, 0.5])

    # Weights cannot exceed one
    @test_throws ArgumentError SamplerWeighted{DNA}([DNA_A, DNA_C], [1.1])
    @test_throws ArgumentError SamplerWeighted{DNA}([DNA_A, DNA_C, DNA_G], [1.00001, 0.0])

    # Weights cannot be negative
    @test_throws ArgumentError SamplerWeighted{DNA}([DNA_A, DNA_C, DNA_G], [0.5, -0.001])
    @test_throws ArgumentError SamplerWeighted{DNA}([DNA_A, DNA_C, DNA_A], [1.1, -0.2])

    # Weights always sum to one after instantiation
    s = SamplerWeighted{DNA}([DNA_A, DNA_C, DNA_G, DNA_W], [0.03, 0.2, 0.7])
    @test sum(s.probs) == 1.0
    s = SamplerWeighted{DNA}([DNA_A], [])
    @test sum(s.probs) == 1.0
    s = SamplerWeighted{DNA}([DNA_A, DNA_C], [1.0])
    @test sum(s.probs) == 1.0
    s = SamplerWeighted{DNA}([DNA_A, DNA_C, DNA_G], [0.03, 0.20000001])
    @test sum(s.probs) == 1.0

    sampler = SamplerWeighted{DNA}([DNA_N, DNA_C, DNA_W, DNA_T], [0.1, 0.2, 0.3])
    firstten = DNA[DNA_C, DNA_N, DNA_T, DNA_T, DNA_T, DNA_N, DNA_W, DNA_C, DNA_W, DNA_N]
    test_sampler(sampler, 0, sampler.elems, firstten, DNA)

    sampler = SamplerWeighted{RNA}([DNA_N, DNA_C, DNA_W, DNA_T], [0.15, 0.5, 0.2])
    firstten = RNA[RNA_C, RNA_N, RNA_U, RNA_W, RNA_U, RNA_N, RNA_C, RNA_C, RNA_C, RNA_N]
    test_sampler(sampler, 0, sampler.elems, firstten, RNA)

    @test_throws MethodError s = SamplerWeighted{AminoAcid}([DNA_A, DNA_C], [0.1])

    # Casual test to see that it can use the global RNG automatically
    sampler = SamplerWeighted{RNA}([DNA_N, DNA_C, DNA_W, DNA_T], [0.15, 0.5, 0.2])
    @test typeof(rand(sampler)) == RNA

end # SamplerWeighted

@testset "randseq Sampler" begin
    # Cannot instantiate < 0-length seq
    @test_throws ArgumentError randseq(DNAAlphabet{2}(), SamplerUniform(dna"ACG"), -1)
    @test_throws ArgumentError randseq(AminoAcidAlphabet(), SamplerUniform(aa"VTW"), -1)
    @test_throws ArgumentError randseq(RNAAlphabet{4}(), SamplerWeighted(dna"ACG", ones(2)/3), -1)

    # CAN make empty sequence with mismatching alphabets
    # Or with matching alphabets
    @test randseq(DNAAlphabet{2}(), SamplerUniform(aa"AV"), 0) == LongSequence{DNAAlphabet{2}}()
    @test randseq(DNAAlphabet{2}(), SamplerUniform(rna"U"), 2) == LongSequence{DNAAlphabet{2}}(dna"TT")

    # Cannot make nonzero sequence with mismatching alphabets
    @test_throws MethodError randseq(RNAAlphabet{4}(), SamplerUniform(aa"AV", 1))
    @test_throws MethodError randseq(AminoAcidAlphabet(), SamplerUniform(dna"ACG", 10))

    # A few samplings
    sampler = SamplerUniform(aa"TVVWYAEDK")
    @test randseq(StableRNG(SEED), AminoAcidAlphabet(), sampler, 10) == aa"VEVTKTTADY"

    sampler = SamplerUniform(dna"TGAWYKN")
    @test randseq(StableRNG(SEED), DNAAlphabet{4}(), sampler, 10) == dna"GWNGGGKTTT"

    sampler = SamplerWeighted(rna"UGCMKYN", [0.1, 0.05, 0.2, 0.15, 0.15, 0.2])
    @test randseq(StableRNG(SEED), DNAAlphabet{4}(), sampler, 10) == dna"CTNYNTKCMT"

    # Casual tests to see that it can use the global RNG automatically
    sampler = SamplerUniform(aa"TVVWYAEDK")
    seq = randseq(AminoAcidAlphabet(), sampler, 25)
    test_isseq(seq, AminoAcidAlphabet, 25)

    sampler = SamplerWeighted(rna"UGCMKYN", [0.1, 0.05, 0.2, 0.15, 0.15, 0.2])
    seq = randseq(RNAAlphabet{4}(), sampler, 25)
    test_isseq(seq, RNAAlphabet{4}, 25)

    sampler = SamplerUniform(dna"AGC")
    seq = randseq(DNAAlphabet{2}(), sampler, 25)
    test_isseq(seq, DNAAlphabet{2}, 25)
end # randseq Sampler

@testset "randseq" begin
    sampler = SamplerUniform(aa"ACDEFGHIKLMNPQRSTVWY")
    automatic = randseq(StableRNG(SEED), AminoAcidAlphabet(), 1000)
    manual = randseq(StableRNG(SEED), AminoAcidAlphabet(), sampler, 1000)
    @test automatic == manual
    @test Set(automatic) == Set(aa"ACDEFGHIKLMNPQRSTVWY")

    @test randseq(StableRNG(SEED), DNAAlphabet{4}(), 20) == dna"CTCTTCGTATGCCGTACCGT"
    @test randseq(StableRNG(SEED), DNAAlphabet{2}(), 20) == dna"CATCCGCTCCTAGGCCATTT"

    @test randseq(StableRNG(SEED), RNAAlphabet{4}(), 20) == rna"CUCUUCGUAUGCCGUACCGU"
    @test randseq(StableRNG(SEED), RNAAlphabet{2}(), 20) == rna"CAUCCGCUCCUAGGCCAUUU"

    # Casual tests to see that it can use the global RNG automatically
    seq = randseq(DNAAlphabet{2}(), 100)
    test_isseq(seq, DNAAlphabet{2}, 100)

    seq = randseq(RNAAlphabet{4}(), 100)
    test_isseq(seq, RNAAlphabet{4}, 100)

    seq = randseq(AminoAcidAlphabet(), 100)
    test_isseq(seq, AminoAcidAlphabet, 100)
end # randseq

@testset "Simple constructors" begin
    manual = randseq(StableRNG(SEED), AminoAcidAlphabet(), 100)
    automatic = randaaseq(StableRNG(SEED), 100)
    @test automatic == manual

    manual = randseq(StableRNG(SEED), DNAAlphabet{4}(), 100)
    automatic = randdnaseq(StableRNG(SEED), 100)
    @test automatic == manual

    manual = randseq(StableRNG(SEED), RNAAlphabet{4}(), 100)
    automatic = randrnaseq(StableRNG(SEED), 100)
    @test automatic == manual

    # Casual tests to see that it can use the global RNG automatically
    seq = randdnaseq(10)
    test_isseq(seq, DNAAlphabet{4}, 10)

    seq = randrnaseq(10)
    test_isseq(seq, RNAAlphabet{4}, 10)

    seq = randaaseq(10)
    test_isseq(seq, AminoAcidAlphabet, 10)
end # simple constructors

@testset "Custom alphabet" begin
    Base.length(A::MyAlphabet) = 6
    BioSequences.symbols(A::MyAlphabet) = (DNA_A, DNA_C, DNA_G, DNA_T, RNA_U, DNA_N)
    BioSequences.BitsPerSymbol(A::MyAlphabet) = BioSequences.BitsPerSymbol{8}()
    BioSequences.encode(A::MyAlphabet, x::DNA) = reinterpret(UInt8, x)
    BioSequences.encode(A::MyAlphabet, x::RNA) = reinterpret(UInt8, x) | 0x10
    Base.eltype(A::MyAlphabet) = NucleicAcid

    function BioSequences.decode(A::MyAlphabet, x::UInt64)
        uint = UInt8(x)
        if uint & 0x10 == 0x10
            return reinterpret(RNA, uint & 0x0f)
        else
            return reinterpret(DNA, uint)
        end
    end

    seq = randseq(MyAlphabet(), 1000)
    @test typeof(seq) == LongSequence{MyAlphabet}
    @test length(seq) == 1000
    @test Set(seq) == Set(symbols(MyAlphabet()))
end # Custom Alphabet
end # Entire Random LongSequences testset


text:
using CreditApprovalStub
using Test

@testset "CreditApprovalStub.jl" begin

# stubs
check_background_success(first_name, last_name) = true
check_background_failure(first_name, last_name) = false

# testing
let first_name = "John", last_name = "Doe", email = "jdoe@julia-is-awesome.com"
    @test open_account(first_name, last_name, email, checker = check_background_success) == :success
    @test open_account(first_name, last_name, email, checker = check_background_failure) == :failure
end

end


text:
using ADSeismic
using ADCMEKit
using ADCME
using PyPlot
using DelimitedFiles
using JLD2
using Statistics
mpi_init()

M = N = Int(round(sqrt(mpi_size())))
n = 1000÷M
@info M
param = MPIAcousticPropagatorParams(NX = n*M, NY = n*N, n = n,
     Rcoef=0.2, DELTAX=10, DELTAY=10, DELTAT=0.05)
compute_PML_Params!(param)

srci = [n*param.M÷5]
srcj = [n*param.N÷2]
srcv = hcat([Ricker(param, 100.0, 500.0)]...)
src = MPIAcousticSource(param, srci,srcj,srcv)

rcvj = Array(20:n*param.M-19)
rcvi = ones(Int64, length(rcvj))*((n*param.N)÷5)
receiver =  MPIAcousticReceiver(param, rcvi, rcvj)

c = 1000*ones(n, n)
center = n*M÷2
width = n*M ÷ 8
for i = 1:n
    for j = 1:n
        ii = i + (param.II-1) * n
        jj = j + (param.JJ-1) * n
        x = i*param.DELTAX + (param.II-1) * n * param.DELTAX
        y = j*param.DELTAY + (param.JJ-1) * n * param.DELTAY
        if center -width <= ii <= center + width && center -width <=jj<=center +width
            c[i, j] = 2000.
        end
    end
end

propagator = MPIAcousticPropagatorSolver(param, src, c)
MPISimulatedObservation!(propagator, receiver)

sess = Session()

timer = []

for i = 1:4
    if mpi_rank()==0
        @info i, mpi_size()
    end
    init(sess)
    d = @timed U = run(sess, propagator.u)
    push!(timer, d[2])
end
t = mean(timer[2:end])

if mpi_rank()==0
    open("data/benchmark.txt", "a") do io
        writedlm(io, [mpi_size() t])
    end
end

mpi_finalize()


text:
using X86Assembly
reload("X86Assembly")

mutable struct GenerativeIOTree <: IO
  stack::Vector{UInt8}
  masks::Vector{UInt8}
  idx::Int
  hadImmediate::Bool
end
GenerativeIOTree() = GenerativeIOTree(UInt8[], UInt8[], 1, false)

struct AbstractImmediate
  size::Int
end
Base.rem(x::AbstractImmediate, ::Type{T}) where {T} = T(0)

struct AbstractBit
  # (-1%UInt8) indicates concrete
  byte_pos::UInt8
  bit_pos::UInt8
  set::Bool
end

struct AbstractByte
  tree::Union{GenerativeIOTree, Void}
  # 1 = low bit
  abstract_bits::NTuple{8, AbstractBit}
end

concrete_bit(x::Integer) = (@assert x == 0 || x == 1; AbstractBit(-1%UInt8,-1%UInt8,x!=0))
isconcrete(x::AbstractBit) = x.byte_pos == -1%UInt8 && x.bit_pos == -1%UInt8
function Base.convert(::Type{AbstractByte}, x::UInt8)
  AbstractByte(nothing, tuple((concrete_bit((x >> i) & 1) for i = 0:7)...))
end
Base.convert(::Type{AbstractByte}, x::Integer) = AbstractByte(convert(UInt8, x))

X86Assembly.representative(x::Integer) = x
function X86Assembly.representative(x::AbstractByte)
  val::UInt8 = 0
  for i = 1:8
    val |= (UInt8(x.abstract_bits[i].set) << (i-1))
  end
  val
end
X86Assembly.representative(x::X86Assembly.Operand) =
  X86Assembly.Operand{UInt8, UInt32}(x.deref, X86Assembly.representative(x.base), X86Assembly.representative(x.scale),
                                              X86Assembly.representative(x.index), X86Assembly.representative(x.disp))
X86Assembly.representative(x::X86Assembly.AVX512Operand) =
  X86Assembly.AVX512Operand(
                X86Assembly.representative(x.mask),
                X86Assembly.representative(x.zero_combining),
                X86Assembly.representative(x.sae),
                X86Assembly.representative(x.er),
                X86Assembly.representative(x.rounding_mode))

function _step!(tree::GenerativeIOTree)
  # addition with carries over ignored bits
  mask = UInt16(tree.masks[end]) | UInt16(0xff) << 8
  val = UInt16(tree.stack[end] + 1)
  while (val & mask) != 0
      if (val & (UInt16(0xff) << 8)) != 0
          pop!(tree.stack)
          pop!(tree.masks)
          !isempty(tree.stack) && _step!(tree)
          return
      end
      carry = val & mask
      val &= ~mask
      val += carry << 1
  end
  tree.stack[end] = val
end

function step!(tree::GenerativeIOTree)
  tree.idx = 1
  tree.hadImmediate = false
  _step!(tree)
end

function Base.read(tree::GenerativeIOTree, ::Type{UInt8})
  @assert !tree.hadImmediate
  tree.idx += 1
  if tree.idx - 1 <= length(tree.stack)
    return AbstractByte(tree, tuple((AbstractBit(tree.idx - 1, i, ((tree.stack[tree.idx - 1] >> (i-1))&1) != 0) for i = 1:8)...))
  end
  push!(tree.stack, 1)
  push!(tree.masks, 0xff)
  bits = tuple((AbstractBit(tree.idx - 1, i, ((tree.stack[end] >> (i-1))&1) != 0) for i = 1:8)...)
  return AbstractByte(tree, bits)
end

import Base: &, <<, |, >>, ==, ~

function (&)(x::AbstractByte, mask::UInt8)
  AbstractByte(x.tree, tuple((((mask >> (i-1)) & 1) == 1 ? x.abstract_bits[i] : concrete_bit(0) for i = 1:8)...))
end
(&)(mask::UInt8, x::AbstractByte) = x & mask
function select_bit(x, y, i)
    xbit = x.abstract_bits[i]; ybit = y.abstract_bits[i]
    (isconcrete(xbit) && xbit.set) && return xbit
    (isconcrete(ybit) && ybit.set) && return ybit
    @assert !(!isconcrete(xbit) && !isconcrete(ybit))
    isconcrete(xbit) ? ybit : xbit
end
function (|)(x::AbstractByte, y::AbstractByte)
  xtree = x.tree
  ytree = y.tree
  (xtree === nothing) && (xtree = ytree)
  (ytree === nothing) && (ytree = xtree)
  @assert xtree === ytree
  AbstractByte(xtree, tuple((select_bit(x, y, i) for i = 1:8)...))
end
function (|)(x::AbstractByte, mask::UInt8)
  AbstractByte(x.tree, tuple((((mask >> (i-1)) & 1) == 1 ? concrete_bit(1) : x.abstract_bits[i] for i = 1:8)...))
end
function (<<)(x::AbstractByte, shift::Integer)
  AbstractByte(x.tree, tuple(((i - shift >= 1) ? x.abstract_bits[i - shift] : concrete_bit(0) for i = 1:8)...))
end
function (>>)(x::AbstractByte, shift::Integer)
  AbstractByte(x.tree, tuple(((i + shift <= 8) ? x.abstract_bits[i + shift] : concrete_bit(0) for i = 1:8)...))
end
~(x::AbstractBit) = AbstractBit(x.byte_pos, x.bit_pos, !x.set)
~(x::AbstractByte) = AbstractByte(x.tree, tuple((~x.abstract_bits[i] for i = 1:8)...))

const IntOrAbstract = Union{Integer, AbstractByte}
function X86Assembly.Operand(x::Bool, base::IntOrAbstract, scale::IntOrAbstract, index::IntOrAbstract, disp::T) where T
    X86Assembly.Operand{AbstractByte, T}(x, base, scale, index, disp)
end

function concretize!(x::AbstractByte)
  for bit in x.abstract_bits
    isconcrete(bit) && continue
    x.tree.masks[bit.byte_pos] &= ~(1 << (bit.bit_pos-1))
  end
  X86Assembly.representative(x)
end
X86Assembly.concretize!(x::AbstractByte) = concretize!(x)

function Base.isless(a::AbstractByte, b::Integer)
  Base.isless(concretize!(a), b)
end
function Base.isless(a::Integer, b::AbstractByte)
  Base.isless(a, concretize!(b))
end
function Base.getindex(a::Dict, b::AbstractByte)
  Base.getindex(a::Dict, concretize!(b))
end
function Base.haskey(a::Dict, b::AbstractByte)
  Base.haskey(a::Dict, concretize!(b))
end
function (==)(x::AbstractByte, y::Integer)
  concretize!(x) == y
end
(==)(x::Integer, y::AbstractByte) = y == x
(::Type{UInt16})(x::AbstractByte) = UInt16(concretize!(x))

function X86Assembly.readImmediate(tree::GenerativeIOTree, T)
  tree.hadImmediate = true
  AbstractImmediate(sizeof(T))
end

Base.eof(tree::GenerativeIOTree) = false

isdone(tree::GenerativeIOTree) = length(tree.stack) >= 1 && tree.stack[1] == 0xff

struct prefix_status{T}
    nprefix_bytes::UInt8
    first_f2f3::UInt8
    last_f2f3::UInt8
    last_valid_seg_prefix::UInt8
    has_lock::Bool
    has_66h::Bool
    has_67h::Bool
    rex::T
    first_non_prefix_byte::T
end

using ProgressMeter

f2f3_choices = UInt8[0x0, 0xF2, 0xF3]
function run!()
  mode = X86Assembly.processor_modes(true, true)
  n = 0
  p = Progress(1616850)
  for last_f2f3 in f2f3_choices
      for first_f2f3 in (last_f2f3 == 0 ? (UInt8(0),) : UInt8[0xF2, 0xF3])
          for last_valid_seg_prefix in UInt8[0x0, 0x64, 0x65]
              for has_lock in (true, false),
                  has_66h in (true, false),
                  has_67h in (true, false),
                  rex in [collect(0x40:0x48); 0]
                    prefix_bytes = UInt8[]
                    first_f2f3 != 0 && push!(prefix_bytes, first_f2f3)
                    (last_f2f3 != first_f2f3) && push!(prefix_bytes, last_f2f3)
                    (last_valid_seg_prefix != 0) && push!(prefix_bytes, last_valid_seg_prefix)
                    has_lock && push!(prefix_bytes, 0xF0)
                    has_66h && push!(prefix_bytes, 0x66)
                    has_67h && push!(prefix_bytes, 0x67)
                    rex != 0 && push!(prefix_bytes, rex)
                    tree = GenerativeIOTree()
                    while !isdone(tree)
                      first_non_prefix_byte = read(tree, UInt8)
                      prefix = X86Assembly.prefix_status(0x0, first_f2f3,
                        last_f2f3, last_valid_seg_prefix, has_lock, has_66h,
                        has_67h, AbstractByte(rex), first_non_prefix_byte)
                      data = X86Assembly._parseOpc(tree, prefix, mode)
                      # Invalid encoding
                      hadImmediate = tree.hadImmediate ? "..." : ""
                      if data == nothing
                          println("$([prefix_bytes; tree.stack])$hadImmediate: Invalid")
                      else
                          print("$([prefix_bytes; tree.stack])$hadImmediate: $data")
                      end
                      step!(tree)
                      next!(p)
                      n += 1
                    end
              end
          end
      end
  end
  println("Number of instruction variants: $n")
end
run!()


text:
using Thebes, Test, Luxor

p1 = Point3D(0, 0, 0)
p2 = Point3D(10, 10, 10)
p3 = Point3D(20, 20, 20)

@test p1 < p2
@test p3 > p2

@test distance(p1, p2) ≈ 17.320508075688775
@test distance(p1, p3) ≈  34.64101615137755

@test between(p1, p2, 2.0) ≈ Point3D(20.0, 20.0, 20.0)
@test between(p1, p2, 1.0) ≈ Point3D(10.0, 10.0, 10.0)
@test between(p1, p2, 0.0) ≈ Point3D(0.0, 0.0, 0.0)

@test midpoint(p1, p2) ≈ Point3D(5.0, 5.0, 5.0)
@test midpoint(p1, p3) ≈ Point3D(10.0, 10.0, 10.0)
@test midpoint(p2, p3) ≈ Point3D(15.0, 15.0, 15.0)
@test midpoint(Point3D(0, 0, 0), Point3D(10, 0, 0)) ≈ Point3D(5.0, 0.0, 0.0)
@test midpoint(Point3D(0, 0, 0), Point3D(0, 10, 0)) ≈ Point3D(0.0, 5.0, 0.0)
@test midpoint(Point3D(0, 0, 0), Point3D(0, 0, 10)) ≈ Point3D(0.0, 0.0, 5.0)

Drawing(1, 1, :png)
@test pin(Point3D(1000, 1000, 1000)) == nothing

pt = pin(Point3D(-1000, -1000, -1000))
@test isapprox(pt.x, 0)
@test isapprox(pt.y, 0)
finish()


text:
function movements!(r, count, side)

    flag = true
    for i = 1:count
        if ismarker(r)
            flag = false
            return flag
        end
        move!(r,side)
    end
end

function corner(side :: HorizonSide)
    side = HorizonSide((Int(side)+1)%4)
    return side
end

function spiral(r)

    count = 1
    side :: HorizonSide = Nord

    while ismarker(r) == false

        flag = movements!(r, count, side)
        side = corner(side)
        if flag == true
            break
        end

        flag = movements!(r,count, side)
        side = corner(side)
        if flag == true
            break
        end

        count = count + 1
    end;

end

spiral(r)
show(r)

max_stars_repo_name: AntonDudDm/mirea-progs

max_issues_repo_name: AntonDudDm/mirea-progs

max_forks_repo_name: AntonDudDm/mirea-progs


text:
# This file is part of GenericSchur.jl, released under the MIT "Expat" license

# The method in this file is derived from LAPACK's zlacon.
# LAPACK is released under a BSD license, and is
# Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.

# Hager's one-norm estimator, with modifications by N.J. Higham
"""
    norm1est!(applyA!,applyAH!,y::Vector) => γ

Estimate the 1-norm of a linear operator `A` expressed as functions which
apply `A` and `adjoint(A)` to a vector such as `y`.

cf. N.J. Higham, SIAM J. Sci. Stat. Comp. 11, 804 (1990)
"""
function norm1est!(applyA!, applyAH!, y::AbstractVector{Ty}; maxiter=5) where Ty
    n = length(y)
    x = fill(one(Ty)/n,n)
    y .= zero(Ty)
    applyA!(x)
    (n == 1) && return abs(x[1])
    γ = norm(x,1)
    tiny = safemin(real(Ty))
    for i=1:n
        absxi = abs(x[i])
        if absxi > tiny
            x[i] /= absxi
        else
            x[i] = one(Ty)
        end
    end
    applyAH!(x)
    ax0,j0 = _findamax(x)
    for iter = 1:maxiter
        x .= zero(Ty)
        x[j0] = one(Ty)
        applyA!(x)
        y .= x
        oldγ = γ
        γ = norm(y,1)
        if γ <= oldγ
            break
        end
        for i=1:n
            absxi = abs(x[i])
            if absxi > tiny
                x[i] /= absxi
            else
                x[i] = one(Ty)
            end
        end
        applyAH!(x)
        jlast = j0
        ax0, j0 = _findamax(x)
        if abs(x[jlast]) == ax0
            break
        end
    end
    # alternative estimate for tricky cases (see Higham 1990)
    v = x # reuse workspace
    isign = 1
    for i in 1:n
        v[i] = isign * (1+(i-1)/(n-1))
        isign = -isign
    end
    applyA!(v)
    t = 2*norm(v,1) / (3*n)
    return max(t,γ)
end

function _findamax(x::AbstractVector{T}) where T
    ax0 = abs(x[1])
    i0 = 1
    for i=2:length(x)
        ax = abs(x[i])
        if ax > ax0
            ax0 = ax
            i0 = i
        end
    end
    return ax0,i0
end

function _findamax(x::AbstractVector{T}) where T <: Complex
    ax0 = abs2(x[1])
    i0 = 1
    for i=2:length(x)
        ax = abs2(x[i])
        if ax > ax0
            ax0 = ax
            i0 = i
        end
    end
    return sqrt(ax0),i0
end


text:
module SpikeGLM
using StatsBase, Reexport
export bsplines, raisedcosines
@reexport using GLM
include("KroneckerProduct.jl")
include("HCatMatrix.jl")
include("basis.jl")

type SpikeGLMPredChol{T,M<:AbstractMatrix,C} <: GLM.LinPred
    X::M                           # model matrix
    beta0::Vector{T}               # base vector for coefficients
    delbeta::Vector{T}             # coefficient increment
    scratchbeta::Vector{T}
    chol::C
end
function SpikeGLMPredChol{T}(X::AbstractMatrix{T})
    SpikeGLMPredChol(X, zeros(T, size(X, 2)), zeros(T, size(X, 2)), zeros(T, size(X, 2)), cholfact(X'X))
end

function GLM.delbeta!{T}(p::SpikeGLMPredChol{T}, r::Vector{T}, wt::Vector{T})
	cholfact!(weightmul!(p.chol.factors, p.X, wt))
	A_ldiv_B!(p.chol, At_mul_B!(p.delbeta, p.X, (r.*wt)))
	p
end

function StatsBase.fit{T<:FloatingPoint,V<:AbstractVector}(::Type{GeneralizedLinearModel},
                                                           X::AbstractMatrix{T}, y::V, d::UnivariateDistribution,
                                                           l::Link=canonicallink(d);
                                                           dofit::Bool=true,
                                                           wts::V=fill!(similar(y), one(eltype(y))),
                                                           offset::V=similar(y, 0), fitargs...)
    size(X, 1) == size(y, 1) || throw(DimensionMismatch("number of rows in X and y must match"))
    n = length(y)
    length(wts) == n || throw(DimensionMismatch("length(wts) does not match length(y)"))
    length(offset) == n || length(offset) == 0 || throw(DimensionMismatch("length(offset) does not match length(y)"))
    wts = T <: Float64 ? copy(wts) : convert(typeof(y), wts)
    off = T <: Float64 ? copy(offset) : convert(Vector{T}, offset)
    mu = GLM.mustart(d, y, wts)
    eta = GLM.linkfun!(l, similar(mu), mu)
    if !isempty(off)
        @inbounds @simd for i = 1:length(eta)
            eta[i] -= off[i]
        end
    end
    rr = GlmResp{typeof(y),typeof(d),typeof(l)}(y, d, l, eta, mu, offset, wts)
    res = GeneralizedLinearModel(rr, SpikeGLMPredChol(X), false)
    dofit ? fit(res; fitargs...) : res
end

end # module


text:
module Laplacian

using SparseArrays
using LinearAlgebra

################################
# Sparse Matrix * Dense Vector #
################################

# https://github.com/JuliaLang/julia/issues/4707

# creates a sparse matrix from its diagonals and generates
# 1D derivatives via finite difference on a staggered grid
function ddx_spdiags(m)
    B = ones(m)*[-1, 1]'
    d = [0, 1]
    n = m + 1
    p = length(d)
    len = zeros(Int, p+1, 1)

    for k = 1:p
        len[k+1] = len[k] + length(max(1,1-d[k]):min(m,n-d[k]))
    end

    a = zeros(Int, len[p+1], 3)

    # Append new d[k]-th diagonal to compact form
    for k = 1:p
        i = max(1,1-d[k]):min(m,n-d[k])
        a[(len[k]+1):len[k+1],:] = [i i.+d[k] B[i.+(m>=n)*d[k],k]]
    end

    return sparse(a[:,1], a[:,2], a[:,3], m, n)
end

function laplace_sparse_matvec(n1, n2, n3)
    I_n1 = sparse(1.0I, n1, n1)
    I_n2 = sparse(1.0I, n2, n2)
    I_n3 = sparse(1.0I, n3, n3)
    D1 = kron(I_n3, kron(I_n2, ddx_spdiags(n1)))
    D2 = kron(I_n3, kron(ddx_spdiags(n2), I_n1))
    D3 = kron(ddx_spdiags(n3), kron(I_n2, I_n1))
    D = [D1 D2 D3] # divergence from faces to cell-centers
    return D*D'
end

function perf_laplace_sparse_matvec(N)
    @assert isinteger(cbrt(N)) "input must have integer cbrt"
    n = Int(cbrt(N)) # makes output sizes match up with other perf_laplace benchmarks
    return laplace_sparse_matvec(n, n, n)
end

#####################
# Iterative Methods #
#####################

function perf_laplace_iter_devec(N)
    u = zeros(N, N)
    u[1, :] .= 1
    Niter = 2^10
    dx2 = 0.1*0.1
    dy2 = dx2
    uout = copy(u)
    for iter = 1:Niter
        for i = 2:N-1
            for j = 2:N-1
                uout[i,j] = ((u[i-1, j] + u[i+1, j])*dy2 + (u[i, j-1] + u[i, j+1])*dx2) * (1 / (2*(dx2+dy2)))
            end
        end
        u, uout = uout, u
    end
    return u
end

function perf_laplace_iter_vec(N)
    u = zeros(N, N)
    u[1,:] .= 1
    Niter = 2^10
    dx2 = dy2 = 0.1*0.1
    for i = 1:Niter
        u[2:N-1, 2:N-1] .= ((u[1:N-2, 2:N-1] .+ u[3:N, 2:N-1]).*dy2 .+ (u[2:N-1, 1:N-2] .+ u[2:N-1, 3:N]).*dx2) .* (1 / (2*(dx2+dy2)))
    end
    return u
end

function perf_laplace_iter_sub(N)
    u = zeros(N, N)
    u[1,:] .= 1
    Niter = 2^10
    dx2 = dy2 = 0.1*0.1
    u0 = view(u, 2:N-1, 2:N-1)
    u1 = view(u, 1:N-2, 2:N-1)
    u2 = view(u, 3:N,   2:N-1)
    u3 = view(u, 2:N-1, 1:N-2)
    u4 = view(u, 2:N-1, 3:N)
    for i = 1:Niter
        u0 .= ((u1 .+ u2).*dy2 .+ (u3 .+ u4).*dx2) .* (1 / (2*(dx2+dy2)))
    end
    return u
end

end # module


text:
"""
    _get_highlights_from_PDF(target::String; concatenate::Bool=true) -> Vector{String}

Extract all highlights from a passed PDF or all PDFs found recursively in the
passed directory.

# Arguments
- `target::String`: $(TARGET_PDF_ARGUMENT)

# Keywords
- `concatenate::Bool=true`: $(CONCATENATE_KEYWORD)

# Returns
- `Vector{String}`: the highlights

# Throws
- Exceptions from: [`get_highlights_comments_pages`](@ref)

# Example
```jldoctest; output = false
using PDFHighlights

path_to_pdf_dir = joinpath(pathof(PDFHighlights) |> dirname |> dirname, "test", "pdf")
path_to_pdf = joinpath(path_to_pdf_dir, "TestPDF.pdf")

PDFHighlights.Internal.PDF._get_highlights_from_PDF(path_to_pdf_dir) ==
PDFHighlights.Internal.PDF._get_highlights_from_PDF(path_to_pdf) ==
String[
    "Highlight 1",
    "Highlight 2 Highlight 3",
    "Highlight 4",
    "Highhighlight 5",
    "6th Highhigh light-",
    "High light 7",
    "8th Highlight-",
]

# output

true
```
"""
function _get_highlights_from_PDF(target::String; concatenate::Bool=true)::Vector{String}
    highlights, _, _ = get_highlights_comments_pages(target; concatenate)
    return highlights
end


text:
printstyled("RandSeis\n", color=:light_green)

fs_range = exp10.(range(-6, stop=4, length=50))
fc_range = exp10.(range(-4, stop=2, length=20))

printstyled("  getbandcode\n", color=:light_green)
for fs in fs_range
  for fc in fc_range
    getbandcode(fs, fc=fc)
  end
end

printstyled("  rand_misc\n", color=:light_green)
for i = 1:100
  D = rand_misc(1000)
end

# Similarly for the yp2 codes
printstyled("  iccodes_and_units\n", color=:light_green)
for i = 1:1000
  cha, u = RandSeis.iccodes_and_units('s', true)
  @test isa(cha, String)
  @test isa(u, String)

  cha, u = RandSeis.iccodes_and_units('s', false)
  @test isa(cha, String)
  @test isa(u, String)
end

# check that rand_t produces sane gaps
printstyled("  rand_t\n", color=:light_green)
fs = 100.0
nx = 100

printstyled("    controlled gaps\n", color=:light_green)
t = RandSeis.rand_t(fs, nx, 10, 1)
@test size(t, 1) == 12
t = RandSeis.rand_t(fs, nx, 0, 1)
@test size(t, 1) == 2
t = RandSeis.rand_t(100.0, 1000, 4, 200)
@test size(t, 1) == 6
@test t[:,1] == [1, 200, 400, 600, 800, 1000]

printstyled("    gap < Δ/2 + 1\n", color=:light_green)
for i in 1:1000
  fs = rand(RandSeis.fs_vals)
  nx = round(Int64, rand(1200:7200)*fs)
  t = RandSeis.rand_t(fs, nx, 0, 1)
  δt = div(round(Int64, sμ/fs), 2) + 1
  gaps = t[2:end-1, 2]

  if length(gaps) > 0
    @test minimum(gaps) ≥ δt
  end
  @test minimum(diff(t[:,1])) > 0
end

printstyled("  rand_resp\n", color=:light_green)
R = RandSeis.rand_resp(1.0, 8)
@test length(R.z) == length(R.p) == 8

printstyled("  namestrip\n", color=:light_green)
str = String(0x00:0xff)
S = randSeisData(3)
S.name[2] = str

for key in keys(bad_chars)
  test_str = namestrip(str, key)
  @test length(test_str) == 256 - (32 + length(bad_chars[key]))
end
redirect_stdout(out) do
  test_str = namestrip(str, "Nonexistent List")
end
namestrip!(S)
@test length(S.name[2]) == 210

printstyled("  repop_id!\n", color=:light_green)
S = randSeisData()
S.id[end] = deepcopy(S.id[1])
id = deepcopy(S.id)
RandSeis.repop_id!(S)
@test id != S.id

printstyled("  randSeis*\n", color=:light_green)
for i = 1:10
  randSeisChannel()
  randSeisData()
  randSeisHdr()
  randSeisSrc()
  randSeisEvent()
end

printstyled("    keywords\n", color=:light_green)
printstyled("      a0\n", color=:light_green)
try
  randSeisData(a0=true)
catch err
  println("a0 = true threw error ", err)
end

printstyled("      c\n", color=:light_green)
for i = 1:10
  S = randSeisData(10, c=1.0)
  @test maximum(S.fs) == 0.0
end

printstyled("      fc\n", color=:light_green)
C = randSeisChannel(fc = 2.0)
resp_a0!(C.resp)
r1 = fctoresp(2.0f0, 1.0f0)
r2 = fctoresp(2.0f0, Float32(1.0/sqrt(2)))
resp_a0!(r1)
resp_a0!(r2)
r = C.resp
for f in (:a0, :f0, :z, :p)
  @test isapprox(getfield(r, f), getfield(r1, f), rtol=eps(Float32)) || isapprox(getfield(r, f), getfield(r2, f), rtol=eps(Float32))
end

printstyled("      s\n", color=:light_green)
for i = 1:10
  S = randSeisData(10, s=1.0)
  @test maximum(S.fs) > 0.0
end


text:
export AltChain, altchain

type AltChain <: Model
  S2 :: Int
  L :: Int
  Jz :: Float64
  Jxy :: Float64
  delta :: Float64
  z :: Tuple{Float64, Float64}
  xy :: Tuple{Complex128, Complex128}
  coeff :: BondCoeffs
  function AltChain(S2::Integer, L::Integer, Jz::Real, Jxy::Real, delta::Real)
    c = bond_coeffs(S2)
    z = (Jz*(1.0+delta), Jz*(1.0-delta))
    xy = (Jxy*(1.0+delta), Jxy*(1.0-delta))
    new(S2, L, Jz, Jxy, delta, z, xy, c)
  end
end
altchain(S2::Integer, L::Integer, Jz::Real, Jxy::Real, delta::Real) = AltChain(S2, L, Jz, Jxy, delta)
altchain(S2::Integer, L::Integer, Jz::Real; Jxy::Real = Jz, delta::Real=0.0) = AltChain(S2, L, Jz, Jxy, delta)
num_sites(chain::AltChain) = chain.L
num_bonds(chain::AltChain) = chain.L

function interall(chain::AltChain, t::Real=0.0)
  const io = open("interall.def", "w")
  header_interall(io, chain)

  ## twisted bond
  bond(io, 0, 1, chain.z[1], chain.xy[1]*cis(t), chain.coeff)

  @inbounds for i in 1:(chain.L-1)
    j = (i+1)%chain.L
    parity = (i%2)+1
    bond(io, i, j, chain.z[parity], chain.xy[parity], chain.coeff)
  end
  close(io)
end

max_stars_repo_name: yomichi/HPhi-samples

max_issues_repo_name: yomichi/HPhi-samples

max_forks_repo_name: yomichi/HPhi-samples