text:
#= -------------------------------------------------------------------------
# @file problema1_sin_disipador.jl
#
# @date 11/23/17 20:37:44
# @author Martin Noblia
# @email mnoblia@disroot.org
#
# @brief
# Simulacion de las ecuaciones difenciales del problema 1 sin disipador
# @detail
#
#  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
---------------------------------------------------------------------------=#
# imports
using DifferentialEquations
using PyPlot
#=------------------------------------------------------------------------------
                           constants
------------------------------------------------------------------------------=#
const  C_j   =  0.5
const  C_c   =  6.8
const  R_jc  =  1.0
const  R_ca  =  30.0
const  T_a   =  30.0
const  T_j_max = 150.0
const  P_d  = (T_j_max - T_a) / (R_jc + R_ca)
a_11 = (-1 / C_c) * ((1 / R_jc) + (1 / R_ca))
a_12 = (1 / (R_jc * C_c))
a_21 = (1 / (C_j * R_jc))
a_22 = (-1/ (C_j * R_jc))
#=------------------------------------------------------------------------------
                        system matrix
------------------------------------------------------------------------------=#
A = [a_11 a_12
     a_21 a_22]

T₀ = [0.0, 0.0] # initial condition
tspan = (0.0, 1500.0) # simulation time
#=------------------------------------------------------------------------------
                           inputs
------------------------------------------------------------------------------=#
step_Pd(t) = (t >=0.0) ? P_d: 0.0
pulse_Pd(t) = (P_d * sin(2.0*π*100.0*t) + P_d) / 2.0
#=------------------------------------------------------------------------------
                  function to modelate the system
------------------------------------------------------------------------------=#
f(t, T) = A * T + [0, step_Pd(t) / C_j] + [T_a / (R_ca * C_c), 0]
prob = ODEProblem(f, T₀, tspan)
sol = solve(prob)
final_value_T_c = sol[1, end]
final_value_T_j = sol[2, end]
println("el valor final de T_J es:", sol[2, end])
#=------------------------------------------------------------------------------
                           Ploting
------------------------------------------------------------------------------=#
plot(sol.t, sol[1, :], label=L"T_c", linewidth=3)
plot(sol.t, sol[2, :], label=L"T_j", linewidth=2)
axhline(final_value_T_j, color="k", linestyle="--")
text(0, final_value_T_j + 3, latexstring("T_J="*"$final_value_T_j"), fontsize=10)
xlabel(L"t\,[seg]")
ylabel(L"T\,[^{\circ}C]")
grid("on")
matplotlib["rcParams"][:update](["font.size" => 14, "font.family" => "serif"])
legend()


text:
"""
    detect_anagrams(subject, candidates)
A function that checks if a list of words
is an [Anagram](https://en.wikipedia.org/wiki/Anagram) or not of a subject word.

is the original word = subject
is list of words to be compared if they are an anagram of `subject` = candidates

```
julia> subject = "listen"
julia> candidates = ["inlets", "enlists", "google", "banana"]
julia> detect_anagrams(subject, candidates)
1-element Vector{String}:
 "inlets"
```

Contributed By:- [Soc V. E.](https://github.com/ReneganRonin)
Based on my exercism's Julia track problem solution on [Anagrams](https://exercism.io/my/solutions/0a1d39faa5694624b7750be3ebf6c974).

Instructions:

An anagram is a rearrangement of letters to form a new word. Given a word and a list of candidates, select the sublist of anagrams of the given word.
Given "listen" and a list of candidates like "enlists" "google" "inlets" "banana" the program should return a list containing "inlets".

Inspired by the [Extreme Startup game](https://github.com/rchatley/extreme_startup)
"""
function detect_anagrams(subject, candidates)
    function condition(i)
        return sort(collect(lowercase(subject))) == sort(collect(lowercase(i)))
    end
    collection = collect(
        word for word in candidates if lowercase(word) != lowercase(subject)
    )
    return filter(condition, collection)
end

ext: jl

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
# L2EpsilonLoss

using LaTeXStrings, LossFunctions, Measures, Latexify

Plots.plot(L2EpsilonInsLoss(1.5), c=:olive, linewidth=2, annotation=(0, 1.5, L"L\left( r \right) = max\{0, |r| - \epsilon \}^2", 15))
Plots.plot!(L2EpsilonInsLoss(1), c=:orangered4, linewidth=2.5)
Plots.plot!(L2EpsilonInsLoss(.5), c=:springgreen4, linewidth=2.5)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
const atomic_masses = Dict([
        ("Ar", 39.948*mu),
        ("Xe", 131.293*mu),
        ])

const atomic_numbers = Dict([
        ("Ar", 1),
        ("Xe", 1),
        ])

@derived_dimension ParticleDensity Unitful.𝐋^-3

"""
    P = Plasma()

Mutable type holding plasma properties like density, temperature, etc.
The plasma is quasineutral, cold-ion.
The show() method prints a table of the properties and derived quantities.

A name can be specified to quickly select the ion mass and charge.
If units are not given, the default is `mi` in kg, `n` in 1/m^3, `Te` in eV, `B` in Gauss.
"""
mutable struct QuasineutralPlasma
    name::String

    np::ParticleDensity
    Te::Energy

    B::BField
    E::EField
    nn::ParticleDensity

    mi::Mass
    Zi::Int
end
QuasineutralPlasma(name::String="Xe",np=1e18u"m^-3",Te=5u"eV";B=0u"Gauss",E=0u"V/m",nn=0u"m^-3") =
    QuasineutralPlasma(
    name,
    check_or_assume(np,u"m^-3"),
    check_or_assume(Te,u"eV"),
    check_or_assume(B,u"Gauss"),
    check_or_assume(E,u"V/m"),
    check_or_assume(nn,u"m^-3"),
    atomic_masses[name],
    atomic_numbers[name],
    )

function Base.show(io::IO,P::QuasineutralPlasma)
    println(io,P.name," quasineutral plasma")
    println(io,"")

    println(io,"Plasma density:")
    println(io,"n = ",P.np)
    println(io,"")

    println(io,"Field properties:")
    println(io,"B = ",P.B)
    println(io,"E = ",P.E)
    println(io,"")

    println(io,"Cold ion properties:")
    println(io,"mi = ",P.mi)
    println(io,"Zi = ",P.Zi," (charge number)")
    println(io,"c_s = ",c_s(P)," (sound speed based on electron temperature)")
    println(io,"f_ci = ",f_ci(P))
    println(io,"r_Lsi = ",r_Lsi(P)," (ion gyroradius based on sound speed)")
    println(io,"")

    println(io,"Electron properties:")
    println(io,"Te = ",P.Te)
    println(io,"c_the = ",c_the(P))
    println(io,"f_ce = ",f_ce(P))
    println(io,"r_Le = ",r_Le(P))
    println(io,"f_pe = ",f_pe(P))
    println(io,"λ_De = ",λ_De(P))
    println(io,"")

    println(io,"Cold neutral properties:")
    println(io,"nn = ",P.nn)
    println(io,"")

    println(io,"Collisional properties:")
    println(io,"ν_ei = ",ν_ei(P))
    println(io,"ν_ie = ",ν_ie(P))
    println(io,"ν_ee = ",ν_ee(P))
    println(io,"χ = ",χ(P)," (Hall parameter based on ν_ei only)")
    println(io,"")
end


text:
"""
    v = InnerProductVec(vec, dotf)

Create a new vector `v` from an existing vector `dotf` with a modified inner product given
by `inner`. The vector `vec`, which can be any type (not necessarily `Vector`) that supports
the basic vector interface required by KrylovKit, is wrapped in a custom struct
`v::InnerProductVec`. All vector space functionality such as addition and multiplication
with scalars (both out of place and in place using `mul!`, `rmul!`, `axpy!` and `axpby!`)
applied to `v` is simply forwarded to `v.vec`. The inner product between vectors
`v1 = InnerProductVec(vec1, dotf)` and `v2 = InnerProductVec(vec2, dotf)` is computed as
`dot(v1, v2) = dotf(v1.vec, v2.vec) = dotf(vec1, vec2)`. The inner product between vectors
with different `dotf` functions is not defined. Similarly, The norm of `v::InnerProductVec`
is defined as `v = sqrt(real(dot(v, v))) = sqrt(real(dotf(vec, vec)))`.

In a (linear) map applied to `v`, the original vector can be obtained as `v.vec` or simply
as `v[]`.
"""
struct InnerProductVec{F,T}
    vec::T
    dotf::F
end

Base.:-(v::InnerProductVec) = InnerProductVec(-v.vec, v.dotf)
Base.:+(v::InnerProductVec{F}, w::InnerProductVec{F}) where {F} =
    InnerProductVec(v.vec + w.vec, v.dotf)
Base.:-(v::InnerProductVec{F}, w::InnerProductVec{F}) where {F} =
    InnerProductVec(v.vec - w.vec, v.dotf)
Base.:*(v::InnerProductVec, a::Number) = InnerProductVec(v.vec * a, v.dotf)
Base.:*(a::Number, v::InnerProductVec) = InnerProductVec(a * v.vec, v.dotf)
Base.:/(v::InnerProductVec, a::Number) = InnerProductVec(v.vec / a, v.dotf)
Base.:\(a::Number, v::InnerProductVec) = InnerProductVec(a \ v.vec, v.dotf)

Base.similar(v::InnerProductVec, ::Type{T} = eltype(v)) where {T} =
    InnerProductVec(similar(v.vec), v.dotf)

Base.getindex(v::InnerProductVec) = v.vec

function Base.copy!(w::InnerProductVec{F}, v::InnerProductVec{F}) where {F}
    copy!(w.vec, v.vec)
    return w
end

function LinearAlgebra.mul!(
    w::InnerProductVec{F},
    a::Number,
    v::InnerProductVec{F}
) where {F}
    mul!(w.vec, a, v.vec)
    return w
end

function LinearAlgebra.mul!(
    w::InnerProductVec{F},
    v::InnerProductVec{F},
    a::Number
) where {F}
    mul!(w.vec, v.vec, a)
    return w
end

function LinearAlgebra.rmul!(v::InnerProductVec, a::Number)
    rmul!(v.vec, a)
    return v
end

function LinearAlgebra.axpy!(
    a::Number,
    v::InnerProductVec{F},
    w::InnerProductVec{F}
) where {F}
    axpy!(a, v.vec, w.vec)
    return w
end
function LinearAlgebra.axpby!(
    a::Number,
    v::InnerProductVec{F},
    b,
    w::InnerProductVec{F}
) where {F}
    axpby!(a, v.vec, b, w.vec)
    return w
end

LinearAlgebra.dot(v::InnerProductVec{F}, w::InnerProductVec{F}) where {F} =
    v.dotf(v.vec, w.vec)
LinearAlgebra.norm(v::InnerProductVec) = sqrt(real(dot(v, v)))


text:
# Has been renamed to avoid it being found by the postprocessing scripts
#
config[@__DIR__] = (
    directory=@__DIR__,
    case="Al",
    n_repeat=12,
    kwargs_scf=(smearing=Smearing.Gaussian(), temperature=0.001, randomise_guess=true),
    firstmixkey="Kerker",
    #
    construct_localiser=(lattice, atoms) -> identity,
    mixings=Dict(
        "Kerker"     => (α, localiser) -> KerkerMixing(α=α, kTF=1.0),
        "None"       => (α, localiser) -> SimpleMixing(α=α),
    ),
)


text:
using MIPVerify
using Gurobi
using JuMP
using Images
using Printf
using MAT
using CSV
using NPZ

# import model
mat_mip_file = "../../provable_robustness_max_linear_regions/models/mmr+at/2019-02-14 23:20:04 dataset=mnist nn_type=fc1 p_norm=inf lmbd=1.0 gamma_rb=0.1 gamma_db=0.1 stage1hpl=10 ae_frac=0.5 epoch=100_mip.mat"
param_dict = mat_mip_file |> matread;

# import dataset
dataset = MIPVerify.read_datasets("MNIST")

# build model
n_h = size(dataset.test.images, 2)
n_w = size(dataset.test.images, 3)
n_c = size(dataset.test.images, 4)
n_out = length(unique(dataset.test.labels))

n_in = n_h * n_w * n_c

n_hidden = 1024
fc1 = get_matrix_params(param_dict, "fc1", (n_in, n_hidden))
softmax = get_matrix_params(param_dict, "softmax", (n_hidden, n_out))

n1 = Sequential([
    Flatten(4),
    fc1, ReLU(interval_arithmetic),
    softmax
    ],
    "outputs/"
)

for i = 1:100

    sample_image = MIPVerify.get_image(dataset.test.images, i);
    sample_label = MIPVerify.get_label(dataset.test.labels, i) + 1;

    d = MIPVerify.find_adversarial_example(n1,
    sample_image,
    sample_label,
    invert_target_selection=true,
    GurobiSolver(OutputFlag=0),
    rebuild=true,
    cache_model=false,
    tightening_algorithm=lp,
    tightening_solver = GurobiSolver(Gurobi.Env(), OutputFlag=0),
    norm_order=Inf);

    perturbed_sample_image = getvalue(d[:PerturbedInput])

    npzwrite("outputs/inputs_for_exact_adversarial_$(@sprintf("%.0f", i)).npz", sample_image)
    npzwrite("outputs/exact_adversarial_$(@sprintf("%.0f", i)).npz", perturbed_sample_image)

end


text:
function htmlOutput(outfile::String = "")
    htmlRegressandTransform(s::String,colmin::Int64,colmax::Int64) = "$s	"
    htmlTableHeader(numberOfResults::Int64, align::String) = "
\n
"
    htmlTableFooter(numberOfResults::Int64, align::String) = "

"
    function htmlHeaderRule(headerCellStartEnd::Vector{Vector{Int64}})
        # if length(headerCellStartEnd)<2
        #     error("Invalid headerCellStartEnd: need to have at least two columns.")
        # end
        # s = ""
        # for i in headerCellStartEnd[2:end]
        #     s = s * "\\cmidrule(lr){$(i[1])-$(i[2])}" * " "
        # end
        # return s
        return ""
    end

    # toprule: just a spacer

    toprule = "

"
    # midrule: a
 with black border on bottom, and a
 spacer
    midrule = "

"
    # bottomrule: a slightly larger spacer
    bottomrule = "

"
    headerrule = htmlHeaderRule
    headercolsep = " "
    colsep = " 	 "
    linestart = "
"
    linebreak = "
"

    label_fe_yes = "Yes"
    label_fe_no = ""

    label_statistic_n = "N"
    label_statistic_r2 = "R2"
    label_statistic_f = "F"
    label_statistic_adjr2 = "Adjusted R2"
    label_statistic_r2_within = "Within-R2"
    label_statistic_p = "F-test p value"
    label_statistic_f_kp = "First-stage F statistic"
    label_statistic_p_kp = "First-stage p value"
    label_statistic_dof = "Degrees of Freedom"

    label_estimator = "Estimator"
    label_estimator_ols = "OLS"
    label_estimator_iv = "IV"
    label_estimator_nl = "NL"

    foutfile = outfile
    encapsulateRegressand = htmlRegressandTransform
    header = htmlTableHeader
    footer = htmlTableFooter
    return RenderSettings(toprule, midrule, bottomrule, headerrule, headercolsep, colsep, linestart,
        linebreak, label_fe_yes, label_fe_no,
        label_statistic_n, label_statistic_r2, label_statistic_adjr2, label_statistic_r2_within,
        label_statistic_f, label_statistic_p, label_statistic_f_kp, label_statistic_p_kp, label_statistic_dof,
        label_estimator, label_estimator_ols, label_estimator_iv, label_estimator_nl,
        foutfile, encapsulateRegressand, header, footer)
end


text:
#!/usr/bin/env julia6
using ReadGlobal, InplaceRealFFT

FFTW.set_num_threads(Threads.nthreads())

function main()

  nx, ny, nz, lx, ly, lz = getdimsize()
  dim = (nx,ny,nz)
  isfile("fftw_wisdom") && FFTW.import_wisdom("fftw_wisdom")
  dtype, padded = checkinput(ARGS[1],nx,ny,nz)
  aux1 = PaddedArray{dtype}(ARGS[1],(nx,ny,nz),padded)
  rfft!(aux1)
  aux1[1] = zero(Complex{dtype})
  irfft!(aux1)
  write(ARGS[2], aux1.data)

  FFTW.export_wisdom("fftw_wisdom")
  return 0
end

main()


text:
for DT in [:Integer, :HalfInteger]
    @eval function Base.print_matrix_row(io::IO,
        X::Union{AbstractHalfIntegerVecOrMat, AdjOrTransAbsHalfIntVecOrMat,
        HalfIntSubArray},
        A::Vector,
        i::$DT, cols::AbstractVector, sep::AbstractString)

        for (k, j) = enumerate(cols)
            k > length(A) && break
            if isassigned(X, i, j)
                x = X[i,j]
                a = Base.alignment(io, x)
                sx = sprint(show, x, context=io, sizehint=0)
            else
                a = Base.undef_ref_alignment
                sx = Base.undef_ref_str
            end
            l = repeat(" ", A[k][1]-a[1]) # pad on left and right as needed
            r = j == axes(X, 2)[end] ? "" : repeat(" ", A[k][2]-a[2])
            prettysx = Base.replace_in_print_matrix(X,i,j,sx)
            print(io, l, prettysx, r)
            if k < length(A); print(io, sep); end
        end
    end
end

function Base.print_matrix(io::IO, X::Union{AbstractHalfIntegerVecOrMat, AdjOrTransAbsHalfIntVecOrMat,
    HalfIntSubArray},
    pre::AbstractString = " ",  # pre-matrix string
    sep::AbstractString = "  ", # separator between elements
    post::AbstractString = "",  # post-matrix string
    hdots::AbstractString = "  \u2026  ",
    vdots::AbstractString = "\u22ee",
    ddots::AbstractString = "  \u22f1  ",
    hmod::Integer = 5, vmod::Integer = 5)

    hmod, vmod = Int(hmod)::Int, Int(vmod)::Int
    if !(get(io, :limit, false)::Bool)
        screenheight = screenwidth = typemax(Int)
    else
        sz = displaysize(io)::Tuple{Int,Int}
        screenheight, screenwidth = sz[1] - 4, sz[2]
    end
    screenwidth -= length(pre)::Int + length(post)::Int
    presp = repeat(" ", length(pre)::Int)  # indent each row to match pre string
    postsp = ""
    @assert textwidth(hdots) == textwidth(ddots)
    sepsize = length(sep)::Int
    rowsA, colsA = UnitRange(axes(X,1)), UnitRange(axes(X,2))
    m, n = length(rowsA)::Int, length(colsA)::Int
    # To figure out alignments, only need to look at as many rows as could
    # fit down screen. If screen has at least as many rows as A, look at A.
    # If not, then we only need to look at the first and last chunks of A,
    # each half a screen height in size.
    halfheight = div(screenheight,2)
    if m > screenheight
        rowsA = [rowsA[(0:halfheight-1) .+ firstindex(rowsA)]; rowsA[(end-div(screenheight-1,2)+1):end]]
    end
    # Similarly for columns, only necessary to get alignments for as many
    # columns as could conceivably fit across the screen
    maxpossiblecols = div(screenwidth, 1+sepsize)
    if n > maxpossiblecols
        colsA = [colsA[(0:maxpossiblecols-1) .+ firstindex(colsA)]; colsA[(end-maxpossiblecols+1):end]]
    end
    A = Base.alignment(io, X, rowsA, colsA, screenwidth, screenwidth, sepsize)
    # Nine-slicing is accomplished using print_matrix_row repeatedly
    if m <= screenheight # rows fit vertically on screen
        if n <= length(A) # rows and cols fit so just print whole matrix in one piece
            for i in rowsA
                print(io, i == first(rowsA) ? pre : presp)
                Base.print_matrix_row(io, X,A,i,colsA,sep)
                print(io, i == last(rowsA) ? post : postsp)
                if i != last(rowsA); println(io); end
            end
        else # rows fit down screen but cols don't, so need horizontal ellipsis
            c = div(screenwidth-length(hdots)::Int+1,2)+1  # what goes to right of ellipsis
            Ralign = reverse(Base.alignment(io, X, rowsA, reverse(colsA), c, c, sepsize)) # alignments for right
            c = screenwidth - sum(map(sum,Ralign)) - (length(Ralign)-1)*sepsize - length(hdots)::Int
            Lalign = Base.alignment(io, X, rowsA, colsA, c, c, sepsize) # alignments for left of ellipsis
            for i in rowsA
                print(io, i == first(rowsA) ? pre : presp)
                Base.print_matrix_row(io, X,Lalign,i,colsA[1:length(Lalign)],sep)
                print(io, (i - first(rowsA)) % hmod == 0 ? hdots : repeat(" ", length(hdots)::Int))
                Base.print_matrix_row(io, X, Ralign, i, (n - length(Ralign)) .+ colsA, sep)
                print(io, i == last(rowsA) ? post : postsp)
                if i != last(rowsA); println(io); end
            end
        end
    else # rows don't fit so will need vertical ellipsis
        if n <= length(A) # rows don't fit, cols do, so only vertical ellipsis
            for i in rowsA
                print(io, i == first(rowsA) ? pre : presp)
                Base.print_matrix_row(io, X,A,i,colsA,sep)
                print(io, i == last(rowsA) ? post : postsp)
                if i != rowsA[end] || i == rowsA[halfheight]; println(io); end
                if i == rowsA[halfheight]
                    print(io, i == first(rowsA) ? pre : presp)
                    print_matrix_vdots(io, vdots, A, sep, vmod, 1, false)
                    print(io, i == last(rowsA) ? post : postsp * '\n')
                end
            end
        else # neither rows nor cols fit, so use all 3 kinds of dots
            c = div(screenwidth-length(hdots)::Int+1,2)+1
            Ralign = reverse(Base.alignment(io, X, rowsA, reverse(colsA), c, c, sepsize))
            c = screenwidth - sum(map(sum,Ralign)) - (length(Ralign)-1)*sepsize - length(hdots)::Int
            Lalign = Base.alignment(io, X, rowsA, colsA, c, c, sepsize)
            r = mod((length(Ralign)-n+1),vmod) # where to put dots on right half
            for i in rowsA
                print(io, i == first(rowsA) ? pre : presp)
                Base.print_matrix_row(io, X,Lalign,i,colsA[1:length(Lalign)],sep)
                print(io, (i - first(rowsA)) % hmod == 0 ? hdots : repeat(" ", length(hdots)::Int))
                Base.print_matrix_row(io, X,Ralign,i,(n-length(Ralign)).+colsA,sep)
                print(io, i == last(rowsA) ? post : postsp)
                if i != rowsA[end] || i == rowsA[halfheight]; println(io); end
                if i == rowsA[halfheight]
                    print(io, i == first(rowsA) ? pre : presp)
                    print_matrix_vdots(io, vdots, Lalign, sep, vmod, 1, true)
                    print(io, ddots)
                    print_matrix_vdots(io, vdots, Ralign, sep, vmod, r, false)
                    print(io, i == last(rowsA) ? post : postsp * '\n')
                end
            end
        end
        if isempty(rowsA)
            print(io, pre)
            print(io, vdots)
            length(colsA) > 1 && print(io, "    ", ddots)
            print(io, post)
        end
    end
end

function Base.show(io::IO, X::AbstractHalfIntegerWrapper)
    show(io, parent(X))
end

function Base.show(io::IO, X::AdjOrTransAbsHalfIntVecOrMat)
    Y = wrapperop(X)(unwraphalfint(X))
    show(io, Y)
end

Base.show(io::IO, i::CartesianIndexHalfInt) = (print(io, "CartesianIndexHalfInt"); show(io, Tuple(i)))
Base.show(io::IO, r::OneTo) = print(io, "HalfIntegerArrays.OneTo(", r.stop, ")")
Base.show(io::IO, r::IdOffsetRange) = print(io,first(r), ':', last(r))

# In general collect the array before displaying
Base.show(io::IO, a::AbstractHalfIntegerArray) = show(io, collect(a))

function Base.show_nd(io::IO, a::AbstractHalfIntegerArray, print_matrix::Function, label_slices::Bool)
    limit::Bool = get(io, :limit, false)
    if isempty(a)
        return
    end
    tailinds = Base.tail(Base.tail(axes(a)))
    nd = ndims(a)-2
    for I in CartesianIndicesHalfInt(tailinds)
        idxs = Tuple(I)
        if limit
            for i = 1:nd
                ii = idxs[i]
                ind = tailinds[i]
                if length(ind) > 10
                    if ii == ind[firstindex(ind)+3] && all(d->idxs[d]==first(tailinds[d]),1:i-1)
                        for j=i+1:nd
                            szj = length(axes(a, j+2))
                            indj = tailinds[j]
                            if szj>10 && first(indj)+2 < idxs[j] <= last(indj)-3
                                @goto skip
                            end
                        end
                        #println(io, idxs)
                        print(io, "...\n\n")
                        @goto skip
                    end
                    if ind[firstindex(ind)+2] < ii <= ind[end-3]
                        @goto skip
                    end
                end
            end
        end
        if label_slices
            print(io, "[:, :, ")
            for i = 1:(nd-1); print(io, "$(idxs[i]), "); end
            println(io, idxs[end], "] =")
        end
        slice = view(a, axes(a,1), axes(a,2), idxs...)
        print_matrix(io, slice)
        print(io, idxs == map(last,tailinds) ? "" : "\n\n")
        @label skip
    end
end

Base.replace_in_print_matrix(A::AbstractMatrix,i::HalfIntegerOrInteger,j::HalfIntegerOrInteger,s::AbstractString) = s
Base.replace_in_print_matrix(A::AbstractVector,i::HalfIntegerOrInteger,j::HalfIntegerOrInteger,s::AbstractString) = s

for DT in [:Integer, :HalfIntegerOrInteger]
    @eval function Base.replace_in_print_matrix(A::AbstractHalfIntegerWrapper{<:Any,2}, i::$DT, j::$DT, s::AbstractString)
        J = to_parentindices(axes(A),(i,j))
        Base.replace_in_print_matrix(parent(A),J...,s)
    end
    @eval function Base.replace_in_print_matrix(A::AbstractHalfIntegerWrapper{<:Any,1}, i::$DT, j::$DT, s::AbstractString)
        i′ = parentindex(Base.axes1(A), i)
        ensureInt(j)
        Base.replace_in_print_matrix(parent(A),i′,unsafeInt(j),s)
    end
end

for DT in [:Integer, :HalfIntegerOrInteger]
    @eval function Base.replace_in_print_matrix(A::AdjOrTransAbsHalfIntVecOrMat, i::$DT, j::$DT, s::AbstractString)
        Base.replace_in_print_matrix(A.parent, j, i, s)
    end
end

function Base.showarg(io::IO, a::HalfIntSubArray, toplevel)
    print(io, "view(")
    Base.showarg(io, parent(a), false)
    iocompact = IOContext(io, :compact => true)
    Base.showindices(iocompact, a.indices...)
    print(io, ')')
    toplevel && print(io, " with eltype ", eltype(a))
end

function Base.showarg(io::IO, a::HalfIntArray, toplevel)
    print(io, "HalfIntArray(")
    Base.showarg(io, parent(a), false)
    if ndims(a) > 0
        iocompact = IOContext(io, :compact => true)
        Base.showindices(iocompact, axes(a)...)
    end
    print(io, ')')
    toplevel && print(io, " with eltype ", eltype(a))
end

function Base.showarg(io::IO, a::SpinMatrix, toplevel)
    print(io, "SpinMatrix(")
    Base.showarg(io, parent(a), false)
    print(io, ", ")
    print(io, a.j)
    print(io, ')')
    toplevel && print(io, " with eltype ", eltype(a))
end

function Base.showarg(io::IO, a::CartesianIndicesHalfIntParent{N}, toplevel) where {N}
    print(io, "CartesianIndicesHalfIntParent{",N,"}")
end

function Base.showarg(io::IO, a::LinearIndicesHalfIntParent{N}, toplevel) where {N}
    print(io, "LinearIndicesHalfIntParent{",N,"}")
end

function Base.showarg(io::IO, r::Iterators.Pairs{<:HalfInteger, <:Any, <:Any, T}, toplevel) where T <: Union{AbstractVector, Tuple}
    print(io, "pairs(::$T)")
end

function Base.showarg(io::IO, r::Iterators.Pairs{<:CartesianIndexHalfInt, <:Any, <:Any, T}, toplevel) where T <: AbstractArray
    print(io, "pairs(::$T)")
end

function Base.showarg(io::IO, r::Iterators.Pairs{<:CartesianIndexHalfInt, <:Any, <:Any, T}, toplevel) where T<:AbstractVector
    print(io, "pairs(IndexCartesian(), ::$T)")
end


text:
using DataStructures
entropy(s::AbstractString) = -sum(x -> x / length(s) * log2(x / length(s)), values(counter(s)))

function fibboword(n::Int64)
    # Initialize the result
    r = Array{String}(n)
    # First element
    r[1] = "0"
    # If more than 2, set the second element
    if n ≥ 2 r[2] = "1" end
    # Recursively create elements > 3
    for i in 3:n
        r[i] = r[i - 1] * r[i - 2]
    end
    return r
end

function testfibbo(n::Integer)
    fib = fibboword(n)
    for i in 1:length(fib)
        @printf("%3d%9d%12.6f\n", i, length(fib[i]), entropy(fib[i]))
    end
    return 0
end

println("  n\tlength\tentropy")
testfibbo(37)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function
makegif(
    u::Matrix,
    Xₙ::Vector,
    Tₘ::Vector,
    ϕl::Matrix = missings(2, 2),
    ϕr::Matrix = missings(2, 2),
    ϕ::Matrix = missings(2,2),
    f1::Vector = missings(2),
    f2::Vector = missings(2);
    frames_to_write::Vector = Vector(),
    title::String = "Анимация чего-то",
    name::String = "animation.gif"
)
    check(x::Array) = !any(ismissing.(x))
    check(x::Tuple) = !any(ismissing.(x))

    @info "Создадим анимацию решения $(name)"

    N, M = size(u) .- 1;
    if isempty(frames_to_write)
        h = div(M, 100);
        @info "Отрисовываем каждый $(h)-ый кадр. Всего 100 кадров"
        frames_to_write = collect(1:h:M+1);
    end

    limits = extrema(u) .* 1.05;
    @info "Устанавливаем пределы графика $(limits)"

    a = Animation()
    @showprogress "Plotting animation..." for m in frames_to_write
        pl = plot(Xₙ, u[:, m], label = "",
             title = title, size = (800, 800),
             xlabel = "X", ylabel = "u(x, t)",
             xlims = extrema(Xₙ), ylims = limits)

        if check(ϕl) && check(ϕr) && check(ϕ)

            # Рисуем вырожденные корни
            plot!(Xₙ, ϕl[:, m], label=L"\phi_l", color=:darkgoldenrod)
            plot!(Xₙ, ϕr[:, m], label=L"\phi_r", color=:darkgoldenrod)
            plot!(Xₙ, ϕ[:, m], label=L"\widetilde{\Phi}", color=:gold)

            # Плавающая по оси Y вслед за пунктиром надпись
            buff = latexstring("f_2", @sprintf("(t) = %.2f",f2[m]));
            annotate!(0.0, f2[m] * 1.25, Plots.text(buff, :left))
            # Горизонтальный пунктир
            plot!([0, f1[m]], [f2[m], f2[m]], line=:dash, color=:black, label="")
            # Красная точка слева, около подписи f2
            scatter!( [0], [f2[m]], color=:red, label="")

            # история f1, f2
            scatter!(f1[1:div(M+1, 80)+1:m], f2[1:div(M+1, 80)+1:m], color=:red, alpha = 0.3, label="")

            # Плавающая по оси X вслед за пунктиром надпись
            buff = latexstring("f_1" , @sprintf("(t) = %.2f",f1[m]));
            annotate!(f1[m] * 1.05, 0.95 * first(limits), Plots.text(buff, :left))
            # Вертикальная линия
            plot!([f1[m], f1[m]], [limits[1], f2[m]], line=:dash, color=:black, label="")
            # Красная точка снизу, около надписи f_1
            # Значение по Y умножаем на нормировочный коэффициент, чтобы точка полностью влезла в кадр
            scatter!( [f1[m]], [limits[1]*0.99], color=:red, label="")

            # красная точка на пересечении пунктиров априорной информации
            scatter!( [f1[m]], [f2[m]], color=:red, label="")
        end

        frame(a)
    end

    if endswith(name, ".gif")
        gif(a, name);
    elseif endswith(name, ".mp4")
        mp4(a, name);
    else
        @error "Неподдерживаемый формат анимации"
    end

end


text:
using LinearAlgebra, SparseArrays
using LogisticOptTools
using BenchmarkTools

const LOT = LogisticOptTools

function dense_callbacks(;n=2_000, p=2_000, intercept=false)
    A = randn(n, p)
    b = sign.(randn(n))
    x = randn(p)

    dat = LOT.LogitData(A, b)
    logit = LOT.LogisticRegressor(dat, LOT.L2Penalty(1.0))
    logit.fit_intercept = intercept

    println("Benchmark f")
    @time LOT.loss(x, dat)
    @time LOT.loss(x, logit)

    println("Benchmark ∇f")
    g = zeros(p)
    @time LOT.gradient!(g, x, dat)
    @time LOT.gradient!(g, x, logit)

    println("Benchmark ∇^2f v")
    hvec = zeros(p)
    vec = zeros(p)
    @time LOT.hessvec!(hvec, x, vec, dat)
    @time LOT.hessvec!(hvec, x, vec, logit)

    println("Benchmark diag(∇^2f)")
    @time LOT.diaghess!(hvec, x, dat)

    println("Benchmark ∇^2f")
    p = 54
    A = randn(p, p)
    b = sign.(randn(p))
    dat = LOT.LogitData(A, b)
    logit = LOT.LogisticRegressor(dat, LOT.L2Penalty(1.0))
    hess = zeros(div(p * (p+1), 2))
    x = randn(p)
    @time LOT.hessian!(hess, x, dat)
    @time LOT.hessian!(hess, x, logit)
end

function sparse_callbacks(;n=2_000, p=2_000)
    # sparse matrix
    A = sprandn(n, p, .25)
    # dense vector
    b = sign.(randn(n))

    dat = LOT.SparseLogitData(A, b)
    logit = LOT.LogisticRegressor(dat, LOT.L2Penalty(1.0))

    n = LOT.nfeatures(dat)
    x = randn(p)
    println("Benchmark f")
    LOT.loss(x, dat)
    @time LOT.loss(x, dat)

    println("Benchmark ∇f")
    g = zeros(p)
    LOT.gradient!(g, x, dat)
    @time LOT.gradient!(g, x, dat)

    println("Benchmark ∇^2f v")
    hvec = zeros(p)
    vec = zeros(p)
    LOT.hessvec!(hvec, x, vec, dat)
    @time LOT.hessvec!(hvec, x, vec, dat)

    #= println("Benchmark ∇^2f") =#
    #= p = 54 =#
    #= A = randn(T, p, p) =#
    #= b = sign.(randn(T, p)) =#
    #= dat = LOT.LogitData(A, b) =#
    #= logit = LOT.GeneralizedLinearModel(dat, LOT.L2Penalty(1.0)) =#
    #= hess = zeros(T, div(p * (p+1), 2)) =#
    #= x = randn(T, p) =#
    #= LOT.hess!(hess, x, dat) =#
    #= @time LOT.hess!(hess, x, dat) =#
    #= @time LOT.hess!(hess, x, logit) =#
end


text:
using Changepoints
using Test

# write your own tests here
include("test_misc.jl")
include("test_ChangepointSampler.jl")
include("test_PELT.jl")
include("test_BS.jl")
include("test_R_comparison.jl")


text:
@generated function _intr(::Val{fname}, out_arg, inp_args...) where {fname,}
    inp_exprs = [:( inp_args[$i] ) for i in 1:length(inp_args)]
    inp_types = [inp_args...]
    out_type = convert(LLVMType, out_arg.parameters[1])

    # create function
    param_types = LLVMType[convert.(LLVMType, inp_types)...]
    llvm_f, _ = create_function(out_type, param_types)
    mod = LLVM.parent(llvm_f)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        # call the intrinsic
        intr_typ = LLVM.FunctionType(out_type, param_types)
        intr = LLVM.Function(mod, string(fname), intr_typ)
        value = call!(builder, intr, [parameters(llvm_f)...])
        ret!(builder, value)
    end

    call_function(llvm_f, out_arg.parameters[1], Tuple{inp_args...}, Expr(:tuple, inp_exprs...))
end

struct GCNIntrinsic
    jlname::Symbol
    rocname::Symbol
    isbroken::Bool # please don't laugh...
    isinverted::Bool
    # FIXME: Input/output types
    inp_args::Tuple
    out_arg::Type
    roclib::Symbol
    suffix::Symbol
end
GCNIntrinsic(jlname, rocname=jlname; isbroken=false, isinverted=false,
             inp_args=(), out_arg=(), roclib=:ocml, suffix=fntypes[first(inp_args)]) =
    GCNIntrinsic(jlname, rocname, isbroken, isinverted, inp_args, out_arg, roclib, suffix)

const MATH_INTRINSICS = GCNIntrinsic[]

for jltype in (
        #= TODO: Float16 Broken due to being i16 in Julia=#
        Float32, Float64)
    append!(MATH_INTRINSICS, GCNIntrinsic.((
        :sin, :cos, :tan, :asin, :acos, :atan, :atan2,
        :sinh, :cosh, :tanh, :asinh, :acosh, :atanh,
        :sinpi, :cospi, :tanpi, :sincospi,
        :asinpi, :acospi, :atanpi, :atan2pi,
        :sqrt, :rsqrt, :cbrt, :rcbrt, :recip,
        :log, :log2, :log10, :log1p, :logb, :ilogb,
        :exp, :exp2, :exp10, :expm1,
        :erf, :erfinv, :erfc, :erfcinv, :erfcx,
        # TODO: :brev, :clz, :ffs, :byte_perm, :popc,
        :isnormal, :nearbyint, :nextafter,
        :pow, :pown, :powr,
        :tgamma, :j0, :j1, :y0, :y1,
    ); inp_args=(jltype,), out_arg=jltype))

    push!(MATH_INTRINSICS, GCNIntrinsic(:sin_fast, :native_sin; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:cos_fast, :native_cos; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:sqrt_fast, :native_sqrt; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:rsqrt_fast, :native_rsqrt; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:recip_fast, :native_recip; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:log_fast, :native_log; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:log2_fast, :native_log2; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:log10_fast, :native_log10; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:exp_fast, :native_exp; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:exp2_fast, :native_exp2; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:exp10_fast, :native_exp10; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:abs, :fabs; inp_args=(jltype,), out_arg=jltype))
    # TODO: abs(::Union{Int32,Int64})

    # FIXME: Multi-argument functions
    #=
    push!(MATH_INTRINSICS, = map(intr->GCNIntrinsic(intr), (
        :sincos, :frexp, :ldexp, :copysign,
    )))
    =#
    #push!(MATH_INTRINSICS, GCNIntrinsic(:ldexp; inp_args=(jltype,), out_arg=(jltype, Int32), isinverted=true))
end

let jltype=Float32
    # TODO: Float64 is broken for some reason, try to re-enable on a newer LLVM
    push!(MATH_INTRINSICS, GCNIntrinsic(:isfinite; inp_args=(jltype,), out_arg=Int32))
    push!(MATH_INTRINSICS, GCNIntrinsic(:isinf; inp_args=(jltype,), out_arg=Int32))
    push!(MATH_INTRINSICS, GCNIntrinsic(:isnan; inp_args=(jltype,), out_arg=Int32))
    push!(MATH_INTRINSICS, GCNIntrinsic(:signbit; inp_args=(jltype,), out_arg=Int32))
end

for intr in MATH_INTRINSICS
    inp_vars = [gensym() for _ in 1:length(intr.inp_args)]
    inp_expr = [:($(inp_vars[idx])::$arg) for (idx,arg) in enumerate(intr.inp_args)]
    libname = Symbol("__$(intr.roclib)_$(intr.rocname)_$(intr.suffix)")
    @eval @inline function $(intr.jlname)($(inp_expr...))
        y = _intr($(Val(libname)), $(intr.out_arg), $(inp_expr...))
        return $(intr.isinverted ? :(1-y) : :y)
    end
end


text:
using Test

include("test_andes_py.jl")


text:
module BrainFlowML

    using CSV
    using DataFrames
    using Statistics
    using DSP
    using LinearAlgebra
    using IterTools

    export each_channel

    const testdata_path = abspath(joinpath(@__DIR__, "../test/testdata"))

    include("savgol.jl")

    mutable struct BioData{T}
        raw::AbstractMatrix{T}
        channels::AbstractVector{Int} # indices of the (emg/eeg) data channel columns
        sample_rate::Int
        labels::AbstractVector{Int} # labels for the gestures, 0 is no gesture
    end
    channels(b::BioData) = b.channels
    raw_data(b::BioData) = b.raw
    channel_data(b::BioData) = view(raw_data(b), :, channels(b))

    # initialize zero labels
    function BioData(raw, channels::AbstractVector{Int}, sample_rate::Int)
        labels = zeros(Int, size(raw, 1))
        return BioData(raw, channels, sample_rate, labels)
    end

    function Base.append!(b1::BioData{T}, b2::BioData{T}) where T
        @assert b1.channels == b2.channels && b1.sample_rate == b2.sample_rate
        b1.raw = vcat(b1.raw, b2.raw)
        append!(b1.labels, b2.labels)
        return nothing
    end

    struct ChannelIterator
        data::BioData
    end
    channels(itr::ChannelIterator) = channels(itr.data)
    raw_data(itr::ChannelIterator) = raw_data(itr.data)

    # never forget: https://docs.julialang.org/en/v1/manual/interfaces/
    function Base.iterate(itr::ChannelIterator, state::Int = 0)
        next_channel_state = iterate(channels(itr), state)
        isnothing(next_channel_state) && return nothing
        next_channel, state = next_channel_state
        return (view(raw_data(itr), :, next_channel), state)
    end
    Base.IteratorSize(itr::ChannelIterator) = Base.IteratorSize(channels(itr))
    Base.length(itr::ChannelIterator) = Base.length(channels(itr))
    Base.size(itr::ChannelIterator) = Base.size(channels(itr))
    Base.eltype(itr::ChannelIterator) = Array{eltype(raw_data(itr)), 1}

    each_channel(b::BioData) = ChannelIterator(b)

    each_channel(A::AbstractMatrix) = (view(A, :, n) for n in axes(A, 2))
    each_channel(A::AbstractMatrix, chans::AbstractVector) = (view(A, :, n) for n in chans)

    function detrend!(b::BioData)
        for chan in each_channel(b)
            detrend!(chan)
        end
    end

    function detrend!(A::AbstractMatrix, chans::AbstractVector = axes(A, 2))
        for chan in each_channel(A, chans)
            detrend!(chan)
        end
    end

    function detrend!(v::AbstractVector)
        v .-= mean(v)
        return nothing
    end

    # calculates the envelope function of each bio data channel
    function smooth_envelope(b::BioData)
        enveloped_data = copy(raw_data(b))
        for (index, chan) in enumerate(each_channel(b))
            enveloped_data[:, index] = smooth_envelope(chan)
        end
        return BioData(enveloped_data, 1:length(channels(b)), b.sample_rate, b.labels)
    end

    function smooth_envelope(b::AbstractVector, window_size::Int=33)
        hilbert_amplitude = abs.(DSP.Util.hilbert(b))
        v = savitzkyGolay(hilbert_amplitude, window_size, 1)
        return v
    end

    function smooth_envelope!(A::AbstractMatrix, chans::AbstractVector = axes(A, 2), window_size::Int=33)
        for chan in each_channel(A, chans)
            chan .= smooth_envelope(chan)
        end
    end

    function sum_channels(b::BioData)
        return reduce(+, each_channel(b))
    end

    # crappy empirical way to automatically label the gestures
    function label_gesture!(b::BioData; shift=0.15, gesture::Int=1)
        smooth_data = smooth_envelope(b)
        summed_data = sum_channels(smooth_data)
        threshold = maximum(summed_data) / 3.5
        labeled = summed_data .> threshold
        sample_shift = Int(round(shift * b.sample_rate))
        widen_range!(labeled, sample_shift)
        b.labels = labeled.*gesture
        return nothing
    end

    # shift the true values left and right
    function widen_range!(v::AbstractVector{Bool}, shift::Int)
        deltas = diff(v)
        len = length(v)
        starts = findall(x->x==1, deltas)
        @inbounds for start_index in starts
            new_start_index = start_index - shift + 1
            new_start_index = new_start_index < 1 ? 1 : new_start_index
            v[new_start_index:start_index] .= true
        end
        endings = findall(x->x==-1, deltas)
        @inbounds for end_index in endings
            new_end_index = end_index + shift
            new_end_index = new_end_index > len ? len : new_end_index
            v[end_index:new_end_index] .= true
        end
    end

    # slice windowed partitions out of A[window, :] and dump into matrix X[:, :, i] for use in algorithms
    # transpose(reshape(X, nsamples*nchannels, npartitions)) for use in algorithms
    function partition_samples(A::AbstractMatrix{T}, sample_size::Int, step_size::Int) where T
        nsamples = size(A, 1)
        nchannels = size(A, 2)
        npartitions = partition_length(nsamples, sample_size, step_size)
        X = zeros(T, sample_size, nchannels,  npartitions)
        @inbounds for (n, part) in enumerate(partition(1:nsamples, sample_size, step_size))
            idx = [x for x in part] # need to convert tuple to array, perhaps this stuff can be done smarter
            slice_of_A = A[idx, :]
            X[:, :, n] .= slice_of_A
        end
        return X
    end

    # get the corresponding label per partition
    function partition_labels(labels::AbstractVector{T}, sample_size::Int, step_size::Int) where T
        nsamples = length(labels)
        npartitions = partition_length(nsamples, sample_size, step_size)
        y = zeros(T, npartitions)
        for (n, p) in enumerate(partition(1:nsamples, sample_size, step_size))
            idx = p[end] # just use last labeled value as the ground truth for now
            y[n] = labels[idx]
        end
        return y
    end

    function partition_samples(b::BioData, sample_size::Int, step_size::Int)
        X = partition_samples(channel_data(b), sample_size, step_size)
        y = partition_labels(b.labels, sample_size, step_size)
        return (X, y)
    end

    function partition_length(nsamples::Int, sample_size::Int, step_size::Int)
        return Int(floor((nsamples-sample_size)/step_size)+1)
    end

    function load_gesture(filepath)
        df = DataFrame(CSV.File(filepath))
        M = Matrix(df) # using a matrix is faster than a dataframe? Could not be generated directly from CSV...
        return BioData(M, 1:8, 600)
    end

    function load_labeled_gesture(file_name, gesture::Int)
        bio_data = load_gesture(file_name)
        detrend!(bio_data)
        label_gesture!(bio_data; gesture=gesture)
        return bio_data
    end

    function load_labeled_gestures(file_names)
        (first_file, rest) = Iterators.peel(file_names)
        bio_data = load_labeled_gesture(first_file, 1)
        for (g, file_name) in enumerate(rest)
            append!(bio_data, load_labeled_gesture(file_name, g+1))
        end
        return bio_data
    end

    # example function to preprocess a brainflow matrix before a model prediction
    # assumes size(A) == (nsamples, nchannels)
    function preprocess_brainflow_data(A::AbstractMatrix, chans::AbstractVector = axes(A, 2))
        nchannels = length(chans)
        nsamples = size(A, 1)
        detrend!(A, chans)
        smooth_envelope!(A, chans)
        return reshape(view(A, :, chans), nsamples*nchannels, 1)
    end

end # module


text:
"""

    gaussian(a...; mean=0.0, std=0.01)

Return a Gaussian array with a given mean and standard deviation.  The
`a` arguments are passed to `randn`.

"""
function gaussian(a...; mean=0.0, std=0.01)
	return randn(a...) * std + mean;
end

"""

    xavier(a...)

Xavier initialization.  The `a` arguments are passed to `rand`.  See
([Glorot and Bengio 2010](http://jmlr.org/proceedings/papers/v9/glorot10a/glorot10a.pdf))
for a description.
[Caffe](http://caffe.berkeleyvision.org/doxygen/classcaffe_1_1XavierFiller.html#details)
implements this slightly differently.
[Lasagne](http://lasagne.readthedocs.org/en/latest/modules/init.html#lasagne.init.GlorotUniform)
calls it `GlorotUniform`.

"""
function xavier(a...)
    w = rand(a...)
    if ndims(w) < 2
        error("ndims=$(ndims(w)) in xavier")
    elseif ndims(w) == 2
        fanout = size(w,1)
        fanin = size(w,2)
    else
        fanout = size(w, ndims(w))
        fanin = div(length(w), fanout)
    end
    s = sqrt(2 / (fanin + fanout))
    w = 2s*w-s
end

"""

Bilinear interpolation filter weights; used for initializing deconvolution layers.

Adapted from https://github.com/shelhamer/fcn.berkeleyvision.org/blob/master/surgery.py#L33

Arguments:

`T` : Data Type

`fw`: Width upscale factor

`fh`: Height upscale factor

`IN`: Number of input filters

`ON`: Number of output filters

Example usage:

w = bilinear(Float32,2,2,128,128)

"""
function bilinear(T,fw,fh,IN,ON)
    @assert fw == fh "Filter must be square"
    @assert IN == ON "Number of input and output filters must be equal"
    f=fw; N=IN;

    sz = 2*f-f%2
    center = f-0.5
    if sz%2 == 1
        center = f-1
    end
    r = linspace(0,sz-1,sz)
    c = linspace(0,sz-1,sz)'

    kernel = (1-abs(r-center)/f)*(1-abs(c-center)/f)
    w = zeros(T,sz,sz,N,N);
    for i=1:N
        w[:,:,i,i] = kernel
    end
    return w
end


text:
function _image_uri(image::AbstractString, tag::AbstractString)
    return "$(image):$(tag)"
end

function _build_docker_image(
    image_uri::AbstractString; dockerfile_path::AbstractString=DEFAULT_DOCKERFILE
)
    cmd = ["docker", "build", "-t", image_uri, "-f", dockerfile_path]

    # Set the appropriate build context
    if dockerfile_path == DEFAULT_DOCKERFILE
        push!(cmd, DEFAULT_PROJECT_PATH)
    else
        # Remove Dockerfile from the path provided, and use that directory for build context
        push!(cmd, dirname(dockerfile_path))
    end

    run(Cmd(cmd))

    return nothing
end


text:
export fold

pre_default(t, path=(); kwargs...) = t

function fold(f, x::Tuple, pre=pre_default, path=(); kwargs...)
    new_x = pre(x, path; kwargs...)
    inner = Tuple(fold(f, steppath(pre, path)(k,v)...; kwargs...) for (k,v) in enumerate(new_x))
    return f(inner, path; kwargs...)
end

function fold(f, x::NamedTuple, pre=pre_default, path=(); kwargs...)
    new_x = pre(x, path; kwargs...)
    inner = Tuple((fold(f, steppath(pre, path)(k,v)...; kwargs...) for (k,v) in pairs(new_x)))
    return f(NamedTuple{keys(new_x)}(inner), path; kwargs...)
end

function fold(f, x, pre=pre_default, path=(); kwargs...)
    new_x = pre(x, path; kwargs...)
    return f(new_x, path; kwargs...)
end

function steppath(pre, path)
    f(k,v) = (v, pre, (path..., k))
end

#######################################

# function example_fold(x)
#     pathsize = 10
#     function pre(x, path)
#         print("↓ path = ")
#         print(rpad(path, pathsize))
#         println("value = ", x)
#         return x
#     end

#     function f(x::Union{Tuple, NamedTuple}, path)
#         print("↑ path = ")
#         print(rpad(path, pathsize))
#         println("value = ", x)
#         return x
#     end

#     function f(x, path)
#         print("↑ path = ")
#         print(rpad(path, pathsize))
#         print("value = ", x)
#         println(" ←-- LEAF")
#         return x
#     end

#     fold(f, x, pre)
# end

#######################################

# q = quote end

# function f(x, path)
#     k = last(path)
#     push!(q.args, :($k = $x))
# end

# function f(x::NamedTuple, path)
#     k = last(path)
#     if isempty(path)
#         return q
#     else
#         push!(q.args, :($k = $x))
#     end
# end

# fold(pre, post, x)

# function flat(x::NamedTuple)
#     f(x::NamedTuple, path)


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Compute the sun radiation received in a satellite surface.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export satellite_sun_radiation_earth_pointing
export satellite_sun_radiation_earth_pointing_mean

"""
    satellite_sun_radiation_earth_pointing(JD0::Number, a::Number, e::Number, i::Number, RAAN::Number, w::Number, numDays::Integer, fN_k::Function, meanAnomaly::Bool = false, step::Number = 0.1*pi/180.0)

Compute the Sun radiation on a satellite surface for an Earth-pointing mission.

# Args

* `JD0`: Initial instant for the analysis [Julian day].
* `a`: Semi-major axis of the orbit [m].
* `e`: Orbit eccentricity.
* `i`: Orbit inclination [rad].
* `w`: Argument of perigee [rad].
* `RAAN`: Right ascension of the ascending node at `JD0` [rad].
* `numDays`: Number of days for the analysis.
* `fN_k`: Function **f(s_b)** that describes the solar panel normal at each k-th
          sampling step. Notice that **s_b** is the Sun vector represented in
          the body coordinate frame.
* `meanAnomaly`: (OPTIONAL) If **true**, compute using angular steps in the mean
                 anomaly instead of in the orbit latitude (**Default**: false).
* `step`: (OPTIONAL) Mean anomaly step (**Default**: 0.1 deg).

# Returns

A matrix containing the Sun radiation [W/m²] for each position in orbit for each
day.

**NOTE**: if the Sun angle is larger than 90 deg or if the satellite is in
eclipse, then `NaN` is returned in the matrix.

# Remarks

The body reference frame is defined as:

* **Z axis** points towards the center of Earth;
* **Y axis** points towards the negative direction of orbit normal;
* **X axis** completes the right-hand reference frame.

If the **mean anomaly** is used, then the average value of the output is the
average sun radiation received by the satellite surface, because every angular
steps have a fixed time interval.

If the **mean anomaly** is used, then the angle interval is [0, 2π]. Otherwise,
the angle interval is [-π,π].

"""
function satellite_sun_radiation_earth_pointing(JD0::Number,
                                                a::Number,
                                                e::Number,
                                                i::Number,
                                                RAAN::Number,
                                                w::Number,
                                                numDays::Integer,
                                                fN_k::Function,
                                                meanAnomaly::Bool = false,
                                                step::Number = 0.1*pi/180.0)
    # Constants
    deg2rad = pi/180.0
    rad2deg = 180.0/pi
    day2sec = 24.0*60.0*60.0

    # Initialization of variables.
    theta = 0.0                   # Sun angle relative to the inertial
                                  # coordinate frame.

    days = collect(0:1:numDays-1) # Vector of the days in which the eclipse time
                                  # will be computed.

    # Angle.
    ang = (!meanAnomaly) ? collect(-pi:step:pi) : collect(0:step:2*pi)

    # Period of an orbit [rad/s].
    n = period(a, e, i, :J2)

    # Step in time
    tstep = step/n

    # Sun radiation.
    sun_radiation = zeros(length(ang),numDays)

    # Perturbations.
    #
    # RAAN rotation rate [rad/s].
    dOmega = dRAAN(a, e, i, :J2)

    # Perturbation of the argument of perigee [rad/s].
    dw = dArgPer(a, e, i, :J2)

    # DCM that rotates the orbit reference frame to the body reference frame.
    #
    # In this case, the body reference frame is defined as:
    #     _ Z axis points towards the center of Earth;
    #     _ Y axis points towards the negative direction of orbit normal;
    #     _ X axis completes the right-hand reference frame.
    # which is common for Earth pointing satellites.

    Dbo = DCM([ 0.0 1.0  0.0;
                0.0 0.0 -1.0;
               -1.0 0.0  0.0])

    # Loop for each day.
    for d in days
        # Get the sun position represented in the Inertial coordinate frame.
        s_i = sun_position_i(JD0+d)
        norm_s_i = norm(s_i)

        # Compute the new orbit parameters due to perturbations.
        w_d    = w + dw*(d*day2sec)
        RAAN_d = RAAN + dOmega*(d*day2sec)

        # Loop through the orbit.
        for k in 1:length(ang)
            # Get the satellite position vector represented in the Inertial
            # coordinate frame.

            if (!meanAnomaly)
                f = ang[k]-w_d
            else
                f = M_to_f(e, ang[k])
            end

            (r_i, rt_i) = satellite_position_i(a, e, i, RAAN_d, w_d, f)

            # Check the lighting conditions.
            lighting = satellite_lighting_condition(r_i, s_i)

            if (lighting == SAT_LIGHTING_SUNLIGHT)
                # Convert the sun vector from the Inertial coordinate frame to
                # the body coordinate frame.
                Doi = angle_to_dcm(RAAN_d, i, w_d+f, :ZXZ)
                s_b = Dbo*Doi*(s_i/norm_s_i)

                # Vector normal to the solar panel.
                N_k = fN_k(s_b)

                # Normalize N_k.
                N_k = N_k/norm(N_k)

                # Compute the sun angle.
                sun_angle_k = acos(dot(s_b,N_k))

                # If the sun angle is larger than 90 deg, then the surface is
                # not illuminated. Thus, the angle will be defined as NaN.
                if (sun_angle_k > pi/2)
                    sun_radiation[k, d+1] = NaN
                else
                    sun_radiation[k, d+1] =
                        sunRad/(4*pi*norm_s_i^2)*cos(sun_angle_k)
                end
            else
                # If the satellite is in eclipse, then the surface is not
                # illuminated. Thus, the angle will be defined as NaN.
                sun_radiation[k,d+1] = NaN
            end
        end
    end

    sun_radiation
end

"""
    satellite_sun_radiation_earth_pointing_mean(JD0::Number, a::Number, e::Number, i::Number, RAAN::Number, w::Number, numDays::Integer, fN_k::Function, step::Number = 0.1*pi/180.0)

Compute the mean Sun radiation on a satellite surface for an Earth-pointing
mission.

# Args

* `JD0`: Initial instant for the analysis [Julian day].
* `a`: Semi-major axis of the orbit [m].
* `e`: Orbit eccentricity.
* `i`: Orbit inclination [rad].
* `w`: Argument of perigee [rad].
* `RAAN`: Right ascension of the ascending node at `JD0` [rad].
* `numDays`: Number of days for the analysis.
* `fN_k`: Function **f(s_b)** that describes the solar panel normal at each k-th
          sampling step. Notice that **s_b** is the Sun vector represented in
          the body coordinate frame.
* `meanAnomaly`: (OPTIONAL) If **true**, compute using angular steps in the mean
                 anomaly instead of in the orbit latitude (**Default**: false).
* `step`: (OPTIONAL) Mean anomaly step (**Default**: 0.1 deg).

# Returns

The mean Sun radiation on a surface [W/m²].

# Remarks

For more details, see **satellite_sun_radiation_earth_pointing**.

"""
function satellite_sun_radiation_earth_pointing_mean(JD0::Number,
                                                     a::Number,
                                                     e::Number,
                                                     i::Number,
                                                     RAAN::Number,
                                                     w::Number,
                                                     numDays::Integer,
                                                     fN_k::Function,
                                                     step::Number = 0.1*pi/180.0)
    sun_radiation = satellite_sun_radiation_earth_pointing(JD0,
                                                           a,
                                                           e,
                                                           i,
                                                           RAAN,
                                                           w,
                                                           numDays,
                                                           fN_k,
                                                           true,
                                                           step)

    # Remove the NaNs.
    sun_radiation[findall(x->(isnan(x)), sun_radiation)] .= 0.0

    # Compute the mean.
    squeeze(mean(sun_radiation, dims = 1), dims = 1)
end

"""
    satellite_sun_radiation_earth_pointing(JD0::Number, a::Number, e::Number, i::Number, RAAN::Number, w::Number, numDays::Integer, N::Vector, meanAnomaly::Bool = false, step::Number = 0.1*pi/180.0)

Compute the Sun radiation on a satellite surface for an Earth-pointing mission.

# Args

* `JD0`: Initial instant for the analysis [Julian day].
* `a`: Semi-major axis of the orbit [m].
* `e`: Orbit eccentricity.
* `i`: Orbit inclination [rad].
* `w`: Argument of perigee [rad].
* `RAAN`: Right ascension of the ascending node at `JD0` [rad].
* `numDays`: Number of days for the analysis.
* `N`: Vector normal to the surface represented in the body reference frame.
* `meanAnomaly`: (OPTIONAL) If **true**, compute using angular steps in the mean
                 anomaly instead of in the orbit latitude (**Default**: false).
* `step`: (OPTIONAL) Mean anomaly step (**Default**: 0.1 deg).

# Returns

A matrix containing the Sun radiation [W/m²] for each position in orbit for each
day.

**NOTE**: if the Sun angle is larger than 90 deg or if the satellite is in
eclipse, then `NaN` is returned in the matrix.

# Remarks

The body reference frame is defined as:

* **Z axis** points towards the center of Earth;
* **Y axis** points towards the negative direction of orbit normal;
* **X axis** completes the right-hand reference frame.

If the **mean anomaly** is used, then the average value of the output is the
average sun radiation received by the satellite surface, because every angular
steps have a fixed time interval.

If the **mean anomaly** is used, then the angle interval is [0, 2π]. Otherwise,
the angle interval is [-π,π].

"""
function satellite_sun_radiation_earth_pointing(JD0::Number,
                                                a::Number,
                                                e::Number,
                                                i::Number,
                                                RAAN::Number,
                                                w::Number,
                                                numDays::Integer,
                                                N::Vector,
                                                meanAnomaly::Bool = false,
                                                step::Number = 0.1*pi/180.0)
    fN_k(x) = N
    satellite_sun_radiation_earth_pointing(JD0, a, e, i, RAAN, w, numDays, fN_k,
                                           meanAnomaly, step)
end

"""
    satellite_sun_radiation_earth_pointing_mean(JD0::Number, a::Number, e::Number, i::Number, RAAN::Number, w::Number, numDays::Integer, N::AbstractVector, step::Number = 0.1*pi/180.0)

Compute the mean Sun radiation on a satellite surface for an Earth-pointing
mission.

# Args

* `JD0`: Initial instant for the analysis [Julian day].
* `a`: Semi-major axis of the orbit [m].
* `e`: Orbit eccentricity.
* `i`: Orbit inclination [rad].
* `w`: Argument of perigee [rad].
* `RAAN`: Right ascension of the ascending node at `JD0` [rad].
* `numDays`: Number of days for the analysis.
* `N`: Vector normal to the surface represented in the body reference frame.
* `meanAnomaly`: (OPTIONAL) If **true**, compute using angular steps in the mean
                 anomaly instead of in the orbit latitude (**Default**: false).
* `step`: (OPTIONAL) Mean anomaly step (**Default**: 0.1 deg).

# Returns

The mean Sun radiation on a surface [W/m²].

# Remarks

For more details, see **satellite_sun_radiation_earth_pointing**.

"""
function satellite_sun_radiation_earth_pointing_mean(JD0::Number,
                                                     a::Number,
                                                     e::Number,
                                                     i::Number,
                                                     RAAN::Number,
                                                     w::Number,
                                                     numDays::Integer,
                                                     N::AbstractVector,
                                                     step::Number = 0.1*pi/180.0)
    sun_radiation = satellite_sun_radiation_earth_pointing(JD0,
                                                           a,
                                                           e,
                                                           i,
                                                           RAAN,
                                                           w,
                                                           numDays,
                                                           N,
                                                           true,
                                                           step)

    # Remove the NaNs.
    sun_radiation[findall(x->(isnan(x)), sun_radiation)] .= 0.0

    # Compute the mean.
    squeeze(mean(sun_radiation, dims = 1), dims = 1)
end


text:
@info "Installing FOOOF"
using PyCall
pip = pyimport("pip")
flags = split(get(ENV, "PIPFLAGS", ""))
ver = get(ENV, "FOOOFVERSION", "1")
packages = ["""fooof$(isempty(ver) ? "" : "==")$(ver)"""]

@info "Package requirements:" packages
@info "Flags for pip install:" flags
ver = isempty(ver) ? "latest" : ver
@info "FOOOF version:" ver
pip.main(["install"; flags; packages])


text:
module TreeSitter

include("native.jl")

const type_regular = TSSymbolTypeRegular
const type_anonymous = TSSymbolTypeAnonymous
const type_auxiliary = TSSymbolTypeAuxiliary

# Julia wrapper around the tree-sitter Parser type.
struct Parser
    c::Ptr{TSParser}
end

# Julia wrapper around the tree-sitter Language type.
struct Language
    c::Ptr{TSLanguage}
end

# Julia wrapper around the tree-sitter Tree type.
struct Tree
    c::Ptr{TSTree}
end

#  Julia wrapper around the tree-sitter TSNode
struct Node
    c::TSNode
end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
import Pkg: PackageSpec

name = "libcgal_julia"
version = VersionNumber(0, 16, julia_version.minor)

isyggdrasil = get(ENV, "YGGDRASIL", "") == "true"
rname = "libcgal-julia"

# Collection of sources required to build CGAL
sources = [
    isyggdrasil ?
        GitSource("https://github.com/rgcv/$rname.git",
                  "1d6eabab1bb94e83676a19f3a6c5877b17b1e11c") :
        DirectorySource(joinpath(ENV["HOME"], "src/github/rgcv/$rname"))
]

# Bash recipe for building across all platforms
jlcgaldir = ifelse(isyggdrasil, rname, ".")
script = raw"""
## pre-build setup
# exit on error
set -eu

macosflags=
case $target in
  *apple-darwin*)
    macosflags="-DCMAKE_CXX_COMPILER_ID=AppleClang"
    macosflags="$macosflags -DCMAKE_CXX_COMPILER_VERSION=10.0.0"
    macosflags="$macosflags -DCMAKE_CXX_STANDARD_COMPUTED_DEFAULT=11"
    ;;
esac
""" * """
## configure build
cmake $jlcgaldir """ * raw"""\
  -B /tmp/build \
  `# cmake specific` \
  -DCMAKE_BUILD_TYPE=Release \
  -DCMAKE_FIND_ROOT_PATH=$prefix \
  -DCMAKE_INSTALL_PREFIX=$prefix \
  -DCMAKE_TOOLCHAIN_FILE=$CMAKE_TARGET_TOOLCHAIN \
  `# tell jlcxx where julia is` \
  -DJulia_PREFIX=$prefix \
  $macosflags

## and away we go..
VERBOSE=ON cmake --build /tmp/build --config Release --target install -- -j$nproc
""" * """
install_license $jlcgaldir/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")
platforms = libjulia_platforms(julia_version)
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libcgal_julia_exact", :libcgal_julia_exact),
    LibraryProduct("libcgal_julia_inexact", :libcgal_julia_inexact),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version)),
    BuildDependency(PackageSpec(name="GMP_jll", version=v"6.1.2")),
    BuildDependency(PackageSpec(name="MPFR_jll", version=v"4.0.2")),

    Dependency(PackageSpec(name="CGAL_jll", version="5.2")),
    Dependency("libcxxwrap_julia_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"8",
               julia_compat = "$(julia_version.major).$(julia_version.minor)")

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
using LinearAlgebra
using DynamicPolynomials
using SwitchOnSafety
using Combinatorics
using SparseArrays
using JuMP, Ipopt, MosekTools,NLopt
using SpecialFunctions

include("../src/RandomTrajectories.jl")
include("../src/AlgebraicLift.jl")
include("../src/ScenarioOpti.jl")
include("../src/ProbabilisticCertificates.jl")
include("../src/WhiteBoxAnalysis.jl")

dim = 2; numMode = 2; dimIn = 1

numScen_budget = 1000

Asave = [zeros(Float64, (dim, dim)) for i=1:numMode]
Bsave = zeros(Float64, (dim, dimIn))
ratio = 1

while true

A = [10*rand(Float64, (dim, dim)).-5 for i=1:numMode]
#A = [[1 1; 0 1], [1 0; 1 1]]
#println(A)
jsrbound = white_box_jsr(A)

if jsrbound > 1 #rhomax < jsrbound*0.9 &&

B = 10*rand(Float64, (dim, dimIn)).-5
gaTrue,K = white_box_stabilization_quad(A,B)
#println(K)
if isempty(K)
    println("A: $A")
    println("B: $B")
end
jsrboundclosed = white_box_jsr([Ai+B*K for Ai in A])

println("A: $Asave")
println("B: $Bsave")
println("Ratio: $ratio")

if gaTrue > 1

    (state0_budget,state_budget) = generate_trajectories(1,A,numScen_budget)

    batchsize = numScen_budget

    K0 = zeros(dimIn,dim)

    K1,ga1,P1 = soslyap_alternating(state0_budget,state_budget;B=B,d=1,batchsize=batchsize,K0=K0,tol=1e-4)
    if ga1 > 1

        K2,ga2,P2 = soslyap_alternating(state0_budget,state_budget;B=B,d=2,batchsize=batchsize,K0=K0,tol=1e-4)

        Aclose1 = [Ai+B*K1 for Ai in A]
        jsrboundclose1 = white_box_jsr(Aclose1)
        Aclose2 = [Ai+B*K2 for Ai in A]
        jsrboundclose2 = white_box_jsr(Aclose2)

        println(repeat('*', 80))
        println("White-box stabilization: $gaTrue")
        println("JSR closed: $jsrboundclosed")
        println("JSR closed 1: $jsrboundclose1")
        println("JSR closed 2: $jsrboundclose2")

        println("JSR closed Prob1: $ga1")
        println("JSR closed Prob2: $ga2")

        println(repeat('*', 80))

        ratiotemp = ga2/ga1

        if ratiotemp < ratio && jsrboundclose2 < jsrboundclose1 && jsrboundclose1 < 1.5*jsrboundclosed && jsrboundclose2 < 0.95*jsrboundclosed
            global Asave = A
            global Bsave = B
            global ratio = ratiotemp
            println("jsr_bound1: $ga1 || K1: $K1")
            println("jsr_bound2: $ga2 || K2: $K2")
        end

        if ratio < 0.7
            println("A: $Asave")
            println("B: $Bsave")
            println("Ratio: $ratio")
            break
        end
    end
end
end

end

max_stars_repo_name: zhemingwang/DataDrivenSwitchControl

max_issues_repo_name: zhemingwang/DataDrivenSwitchControl

max_forks_repo_name: zhemingwang/DataDrivenSwitchControl


text:
"""
    gmesh2fem(name)

Extract points and elements from Gmesh file `name.msh`. This will create the two files
`name_nodes.txt` and `name_elements.txt` in the same directory.
"""
function gmesh2fem(name)
    # Read points and elements
    mshname = name * ".msh"
    @info "Reading nodes and elements from " * mshname
    points, elements = open(mshname, "r") do io
        # Points
        for i = 1:4
            readline(io)
        end
        npoint = parse(Int, readline(io))
        points = zeros(3, npoint)
        for ipoint = 1:npoint
            points[:, ipoint] = parse.(Float64, split(readline(io))[2:4])
        end
        readline(io)

        # Elements
        readline(io)
        nelement = parse(Int, readline(io))
        elements = zeros(Int, 4, nelement)
        for ielement = 1:nelement
            tmp = parse.(Int, split(readline(io)))
            elements[:, ielement] = tmp[6:9]
        end

        points, elements
    end

    # Write nodes
    nodename = name * "_nodes.txt"
    @info "Writing nodes to " * nodename
    npoint = size(points, 2)
    open(nodename, "w") do io
        for i = 1:npoint
            for j = 1:3
                write(io, "  ")
                write(io, string(points[j, i]))
            end
            write(io, "\n")
        end
    end

    # Write elements
    elementname = name * "_elements.txt"
    @info "Writing elements to " * elementname
    nelement = size(elements, 2)
    open(elementname, "w") do io
        for i = 1:nelement
            for j = 1:4
                write(io, "  ")
                write(io, string(elements[j, i]))
            end
            write(io, "\n")
        end
    end

    nothing
end


text:
"""
    ServerOutput(frames, rulset, args...; port=8080, kwargs...)

A basic Mux.jl webserver, serving a [`InteractOutput`](@ref)s to the web.

Unlike ElectronOutput, the parameter modifications are not written back
to the original rulset, and the simulations are not stored.
Each page load gets a identical initialised rulset.

### Arguments
- `init`: `AbstractArray` or `NamedTuple` of `Array`
- `ruleset::Models`: tuple of rulset wrapped in Models().

### Keyword arguments
- `port`: port number to reach the server. ie localhost:8080
- `kwargs`: keyword arguments to be passed to [`InteractOuput`](@ref).
"""
mutable struct ServerOutput{I}
    init::I
    port::Int
end
function ServerOutput(init; port=8080, ruleset, kwargs...)
    server = ServerOutput(init, port)
    function app(request)
        InteractOutput(deepcopy(server.init); ruleset=deepcopy(ruleset), kwargs...).page
    end
    WebIO.webio_serve(Mux.page("/", request -> app(request)), port)
    return server
end

function Base.show(io::IO, output::ServerOutput)
    println(io, typeof(output))
    println(io, "Port: ", output.port)
end


text:
module Quake
using Blosc, Dates, DSP, LightXML, LinearAlgebra, Printf, SeisIO, SeisIO.FastIO, Sockets
using HTTP: request, Messages.statustext
Blosc.set_compressor("lz4")
Blosc.set_num_threads(Sys.CPU_THREADS)
path = Base.source_dir()

const tracefields = (:az, :baz, :dist, :id, :loc, :fs, :gain, :misc, :name, :notes, :pha, :resp, :src, :t, :units, :x)
const loc_qual_fields = (:se, :gap, :dmin, :dmax)
const loc_qual_names = ("standardError", "azimuthalGap", "minimumDistance", "maximumDistance")

# imports
include("Quake/imports.jl")

# types for earthquake data
include("Quake/Types/EQLoc.jl")
include("Quake/Types/EQMag.jl")
include("Quake/Types/SourceTime.jl")
include("Quake/Types/SeisSrc.jl")
include("Quake/Types/SeisPha.jl")
include("Quake/Types/PhaseCat.jl")
include("Quake/Types/SeisHdr.jl")
include("Quake/Types/EventTraceData.jl")
include("Quake/Types/EventChannel.jl")
include("Quake/Types/SeisEvent.jl")

# formats
for i in ls(path*"/Quake/Formats/")
  if endswith(i, ".jl")
    include(i)
  end
end

# processing
for i in ls(path*"/Quake/Processing/")
  if endswith(i, ".jl")
    include(i)
  end
end

# utilities
for i in ls(path*"/Quake/Utils/")
  if endswith(i, ".jl")
    include(i)
  end
end

# web
for i in ls(path*"/Quake/Web/")
  if endswith(i, ".jl")
    include(i)
  end
end

# exports
include("Quake/exports.jl")

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -7.9284408379108005 <= q <= 14.377266175441108)

                     add_NL_constraint(m, :(softplus(0.4343485145630286 + -0.5735104656090533*softplus(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2])) + -0.40833038450807946*softplus(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2])) + 0.09714914815406805*softplus(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2])) + -0.19623354791767778*softplus(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2])) + -0.3237276355524821*softplus(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(x[2]))) + softplus(-0.11905941659917563 + 0.17224590981590104*softplus(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2])) + 0.8406566605626429*softplus(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2])) + 0.24719777641241691*softplus(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2])) + -0.24308200319021767*softplus(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2])) + -0.39557038452022475*softplus(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(x[2]))) + softplus(0.967203829156619 + 0.3785230261408685*softplus(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2])) + -0.23118999327456224*softplus(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2])) + -0.5197073892915105*softplus(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2])) + 0.6005003791404913*softplus(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2])) + 0.15463985115916623*softplus(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(x[2]))) + softplus(-0.18340968019959902 + 0.01980109053282808*softplus(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2])) + 0.8767310690037822*softplus(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2])) + 0.42243485881005904*softplus(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2])) + -0.1694649191201152*softplus(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2])) + 0.5278959282338329*softplus(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(x[2]))) + softplus(0.7043455967116397 + -0.48930630508951234*softplus(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2])) + 0.7375930151357304*softplus(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2])) + -0.9188673755139103*softplus(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2])) + 0.5801272291202126*softplus(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2])) + 0.06743755305487698*softplus(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(x[2]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""

```julia
Atom::DataType
```

Structure that contains the atom properties. It is mutable, so it can be edited.
Fields:

```julia
mutable struct Atom
  index::Int # The sequential index of the atoms in the file
  index_pdb::Int # The index as written in the PDB file (might be anything)
  name::String # Atom name
  resname::String # Residue name
  chain::String # Chain identifier
  resnum::Int # Number of residue as written in PDB file
  residue::Int # Sequential residue (molecule) number in file
  x::Float64 # x coordinate
  y::Float64 # y coordinate
  z::Float64 # z coordinate
  beta::Float64 # temperature factor
  occup::Float64 # occupancy
  model::Int # model number
  segname::String # Segment name (cols 73:76)
end
```

### Example

```julia-repl
julia> pdb = wget("1LBD");

julia> printatom(pdb[1])
   index name resname chain   resnum  residue        x        y        z  beta occup model segname index_pdb
       1    N     ALA     P        1        1    2.062  -13.995   21.747  0.00  1.00     1    PROT         1

julia> pdb[1].resname
"ALA"

julia> pdb[1].chain
"P"

julia> element(pdb[1])
"N"

julia> mass(pdb[1])
14.0067

```

"""
Base.@kwdef mutable struct Atom
    index::Int = 0 # The sequential index of the atoms in the file
    index_pdb::Int = 0 # The index as written in the PDB file (might be anything)
    name::String = "X"
    resname::String = "XXX"
    chain::String = "X"
    resnum::Int = 0 # Number of residue as written in PDB file
    residue::Int = 0 # Sequential residue (molecule) number in file
    x::Float64 = 0.0
    y::Float64 = 0.0
    z::Float64 = 0.0
    beta::Float64 = 0.0
    occup::Float64 = 0.0
    model::Int = 0
    segname::String = "XXXX" # Segment name (cols 73:76)
end

index(atom::Atom) = atom.index
index_pdb(atom::Atom) = atom.index_pdb
name(atom::Atom) = atom.name
resname(atom::Atom) = atom.resname
chain(atom::Atom) = atom.chain
resnum(atom::Atom) = atom.resnum
residue(atom::Atom) = atom.residue
beta(atom::Atom) = atom.beta
occup(atom::Atom) = atom.occup
model(atom::Atom) = atom.model
segname(atom::Atom) = atom.segname

const atom_title = @sprintf(
    "%8s %4s %7s %5s %8s %8s %8s %8s %8s %5s %5s %5s %7s %9s",
    "index",
    "name",
    "resname",
    "chain",
    "resnum",
    "residue",
    "x",
    "y",
    "z",
    "beta",
    "occup",
    "model",
    "segname",
    "index_pdb"
)
atom_line(atom::Atom) = @sprintf(
    "%8i %4s %7s %5s %8i %8i %8.3f %8.3f %8.3f %5.2f %5.2f %5i %7s %9i",
    atom.index,
    atom.name,
    atom.resname,
    atom.chain,
    atom.resnum,
    atom.residue,
    atom.x,
    atom.y,
    atom.z,
    atom.beta,
    atom.occup,
    atom.model,
    atom.segname,
    atom.index_pdb
)

function printatom(atom::Atom)
    println(atom_title)
    println(atom_line(atom))
end

#
# Print a formatted list of atoms
#
function print_short_atom_list(atoms::AbstractVector{Atom})
    println(atom_title)
    for i = 1:min(length(atoms), 3)
        println(atom_line(atoms[i]))
    end
    if length(atoms) > 7
        @printf("%57s\n", "⋮ ")
    end
    for i = max(4, length(atoms) - 2):length(atoms)
        println(atom_line(atoms[i]))
    end
end

function Base.show(io::IO, atom::Atom)
    println(atom_line(atom))
end

function Base.show(io::IO, ::MIME"text/plain", atoms::AbstractVector{Atom})
    println("   Array{Atoms,1} with $(length(atoms)) atoms with fields:")
    print_short_atom_list(atoms)
end

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
# ubSmoteExs.R
using DataFrames

# d = readtable("./data/people.csv", makefactors = true)

function factor_cols(dat::DataFrame)
    p = size(dat, 2)
    are_factors = falses(p)
    for j = 1:p
        if isa(dat[:, j], NullableCategoricalArray)
            are_factors[j] = true
        end
    end
    indcs = find(are_factors)
    return indcs
end

# @code_warntype factor_cols(d)

function factor_to_float(v)
    unique_cats = levels(v)         # unique categories
    sort!(unique_cats)
    cat_dictionary = Dict{Nullable{String}, Float64}()
    val = 1.0
    for k in unique_cats
        cat_dictionary[Nullable(k)] = val
        val += 1.0
    end
    n = length(v)
    out = Array{Float64, 1}(n)
    for i = 1:n
        out[i] = cat_dictionary[v[i]]
    end
    out
end

function float_to_factor(v, levels)
    sort!(levels)
    str_vect = map(x -> levels[Int(round(x))], v)
    out = CategoricalArray(str_vect)
    return out
end

# This function behaves a bit like R's scale()
# function when it's call with MARGIN = 2.
function rscale(X, center, scale)
    n, p = size(X)
    out = Array{Float64, 2}(n, p)
    for i = 1:n
        for j = 1:p
            out[i, j] = (X[i, j] - center[j])/scale[j]
        end
    end
    out
end

function column_ranges{T<:Real}(X::Array{T, 2})
    p = size(X, 2)
    ranges = Array{Float64,1}(p)

    for j = 1:p
        ranges[j] = maximum(X[:, j]) - minimum(X[:, j])
    end
    ranges
end

function smote_exs(dat::DataFrame, tgt::Symbol, pct = 200, k = 5)
    n, m = size(dat)
    T = Array{Float64, 2}(n, m-1)

    # HACK: We assume outcome var is last column
    factor_indcs = factor_cols(dat)[1:end-1]

    for j = 1:size(T, 2)
        if j ∈ factor_indcs
            T[:, j] = factor_to_float(dat[:, j])
        else
            T[:, j] = convert(Array{Float64,1}, dat[:, j])
        end
    end

    # when pct < 100, only a percentage of cases will be SMOTEd
    if pct < 100
        n_needed = round(Int, (pct/100)*n)
        idx = sample(1:n, n_needed)
        T = T[idx, :]
        pct = 100
    end
    n, p = size(T)
    # display(T)
    ranges = column_ranges(T)

    n_exs = round(Int, floor(pct/100))   # num. of artificial ex for each member of T
    xnew = Array{Float64, 2}(n_exs*n, p)

    for i = 1:n

        # the k nearest neighbors of case T[i, ]
        xd = rscale(T, T[i, :], ranges)

        for col in factor_indcs
            xd[:, col] = map(x -> x == 0.0 ? 1.0 : 0.0, xd[:, col])
        end

        dd = xd.^2 * ones(p)
        kNNs = sortperm(dd)[2:(k+1)]

        for l = 1:n_exs
            neighbor = sample(1:k)
            ex = Array{Float64, 1}(p)

            # the attribute values of generated case
            difs = T[kNNs[neighbor], :] - T[i, :]
            xnew[(i - 1)*n_exs + l, :] = T[i, :] + rand()*difs

            # For each of factor variable, sample at random the original value
            # of Person i or the value that one of Person i's nearest neighbors has.
            for col in factor_indcs
                xnew[(i - 1)*n_exs + l, col] = sample(vcat(T[kNNs[neighbor], col], T[i, col]))
            end
        end
    end

    new_cases = DataFrame()
    for j = 1:p
        if j ∈ factor_indcs
            new_cases[:, j] = float_to_factor(xnew[:, j], levels(dat[:, j]))
        else
            new_cases[:, j] = xnew[:, j]
        end
    end
    yval = String(dat[1, tgt].value)
    new_cases[:, tgt] = CategoricalArray(repeat([yval], inner = n_exs*n))
    return new_cases
end

# This version of the function is to be used when we have no factor
# variables. And it assumes input is simply a numeric matrix, where
# the last column is the outcome (or target) variable.
# NOTE: `pct` is the pctent of positive examples relative to total
# sample size to be returned.
function smote_exs{S<:Number}(dat::Array{S, 2}, tgt::Int, pct = 200, k = 5)
    if pct < 1
        warn("Percent over-sampling cannot be less than 1\n
              Setting `pct` to 1.")
        pct = 1
    end

    n, m = size(dat)
    T = Array{Float64, 2}(n, m-1)

    for j = 1:size(T, 2)
        T[:, j] = convert(Array{Float64,1}, dat[:, j])
    end

    # When pct < 100, only a percentage of cases will be SMOTEd
    if pct < 100
        n_needed = round(Int, (pct/100)*n)
        idx = sample(1:n, n_needed)
        T = T[idx, :]
        pct = 100
    end

    n, p = size(T)
    # display(T)
    ranges = column_ranges(T)

    n_exs = round(Int, pct/100)   # num. of artificial ex for each member of T
    xnew = Array{Float64, 2}(n_exs*n, p)

    for i = 1:n

        # The k nearest neighbors of case T[i, ]
        xd = rscale(T, T[i, :], ranges)

        dd = xd.^2 * ones(p)
        kNNs = sortperm(dd)[2:(k+1)]

        for l = 1:n_exs
            neighbor = sample(1:k)
            ex = Array{Float64, 1}(p)

            # The attribute values of generated case
            difs = T[kNNs[neighbor], :] - T[i, :]
            xnew[(i - 1)*n_exs + l, :] = T[i, :] + rand()*difs
        end
    end
    # Find what the minority class is in outcome
    yval = dat[1, tgt]
    new_cases = hcat(xnew, repeat([yval], inner = n_exs*n))
    return new_cases
end

# m = 150
# X = rand(m, 10)
# y = ones(m)
# X = hcat(X, y)
#
# smote_exs(X, 11)

function cases_needed{T<:Real}(y::Array{T, 1})
    n_minority = count(x -> x == 1, y)
    n = length(y)
    0.5n - n_minority
end

function pct_needed{T<:Real}(y::Array{T, 1})
    numer = cases_needed(y)
    denom = count(x -> x == 1, y)
    return 100 * numer/denom
end

# w1 = [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
#
# cases_needed(w1)
# pct_needed(w1)
#
#
#
# X = randn(100, 10)
# y = vcat(zeros(90), ones(10))
# ub_smote(X, y)


text:
using BcdV3
using StaticArrays
using MimuConstants
using LinearAlgebra
using Printf
using Dates

function crb(w, r, Q_inv)
    r_m = hcat(r...)
    Q_inv_tot = cat(Q_inv...; dims = (1,2))
    Nt =  length(Q_inv)
    Na = length(r)
    Ng = Nt  - Na
    H_m = zeros(3*Nt, 6)
    BcdV3.fill_H_constants!(H_m, Na)
    BcdV3.H!(H_m, r)
    crb = zeros(3, length(w))
    J_m = zeros(3*Nt, 3)

    BcdV3.fill_J_h_buff!(J_m, Na, Ng)
    for n = 1:length(w)
        BcdV3.J_h!(J_m, w[n], r_m, Ng)
        J_tot = hcat(J_m, H_m)
        fim_n = J_tot'*Q_inv_tot*J_tot

        fim_n_inv = inv(fim_n)
        crb[:,n] = diag(fim_n_inv)[1:3]
    end
    return crb
end
function get_callback_trace(trace, trace_p, show_every = 1)
    function callback(k::Int, dx, dlogp, log_p, r, T, b)
        trace[:f][k] = log_p
        trace[:f_abs][k] = dlogp
        trace[:f_rel][k] = abs(dlogp)/abs(log_p)

        trace[:x_abs][k] = dx

        if mod(k,show_every) == 0
            time = Dates.format(Dates.now(), "yyyy-mm-dd HH:MM:SS")
            @printf("%20s %7d logp %.16e dlogp % .8e ||dx|| %.8e\n",
                    time, k, log_p, dlogp, trace[:x_abs][k])

            # C = 0.0
            # for k = 1:8
            #     C += tr(T[k+8])
            # end
            # C /= (3*8)
            # @printf(" C_norm %.8e\n", C)
            push!(trace_p[:r], Array(r))
            push!(trace_p[:T], Array.(T))
            push!(trace_p[:b], Array.(b))
            push!(trace_p[:k], k)
        end
    end
    return callback
end

# Noise less simulations
TT = Float64
r_m = [1.0  0.0 0.0;
-1.0  0.0 0.0;
0.0  1.0 0.0;
0.0 -1.0 0.0]'.*1.0e-2 # [m]
r_m .-= r_m[:,1]
Na = size(r_m, 2)
r = [SVector{3,TT}(r_m[:,k]) for k = 1:Na]

# sig_a = get_rms_noise_acc()
# sig_g = get_rms_noise_gyro()
sig_a = 0.01
sig_g = deg2rad(1.0)

Ng = 4

N = 100 # Dyn
Ns = []
show_every = 100
tol_dlogp = 1.0e-9

Nt = Na + Ng
g_mag = 9.81
Random.seed!(1)

b = [@SVector randn(3) for k = 1:Nt]
T = [I + 1.0e-3*@SMatrix randn(3,3) for k = 1:Nt]
T[1] = I + 1.0e-3*SMatrix{3,3, Float64}(triu(randn(3,3)))

# T = [I + @SMatrix zeros(3,3) for k = 1:Nt]
# b = [@SVector zeros(3) for k = 1:Nt]
# r = [SVector{3}(r_m[:,k])*100 for k = 1:Na]

w_norm = map(x-> 10^x, LinRange(2, 4, N)) .|> deg2rad
# w = [deg2rad(1000.0) * @SVector randn(3) for _ = 1:N]
w_rot = SVector{3}([1.0, 0.0, 0.0])
w = [w_rot*w_norm[n] for n = 1:N]
w_m = hcat(Array.(w)...)
w_dot = [@SVector randn(3) for _ = 1:N]
s = [@SVector randn(3) for _ = 1:N]

y_d_true = Array{SVector{3,TT}}(undef, Nt,N)
y_d_true_2 = Array{SVector{3,TT}}(undef, Nt,N)
for n = 1:N
    for k = 1:Na
        y_d_true[k,n] = T[k]*(Ω(w[n])^2*r[k] + Ω(w_dot[n])*r[k] + s[n]) + b[k]
        y_d_true_2[k,n] = T[k]*(Ω(2*w[n])^2*r[k] + Ω(2*w_dot[n])*r[k] + 2*s[n]) + b[k]
    end
    for k in range(Na+1, length=Ng)
        y_d_true[k,n] = T[k]*w[n] + b[k]
        y_d_true_2[k,n] = T[k]*w[n]*2 + b[k]
    end
end

Q_inv = [SMatrix{3,3,Float64}(I)/( k <= Na ? sig_a^2 :  sig_g^2) for k = 1:Nt]

N_mc = 10^4
dw = zeros(3, N)
dw2 = zeros(3, N)
show_all(x) = show(stdout, "text/plain", x)
Juno.@progress for n = 1:N_mc
    y = Array{SVector{3,TT}}(undef, Nt,N)
    y2 = Array{SVector{3,TT}}(undef, Nt,N)
    for n = 1:N
        for k = 1:Na
            e = sig_a * randn(3)
            y[k,n] = y_d_true[k,n] + e
            y2[k,n] = y_d_true_2[k,n] + e
        end
        for k in range(Na+1, length=Ng)
            e = sig_g * randn(3)
            y[k,n] = y_d_true[k,n] + e
            y2[k,n] = y_d_true_2[k,n] + e
        end
    end

    w0 = deepcopy(w)
    r_hat, T_hat, b_hat, eta = bcd!(
    copy(r), copy(T), copy(b), y, Q_inv, Ns, g_mag, Val{3*(Na+Ng)}(),Val{Na}(),
    1.0e-12, 100, 1.0e-8, 200, tol_dlogp, 1.0-14, 1; warn = false, w0 = w0
    )
    w_hat_m = hcat(Array.(eta[:w])...)
    dw .+= (w_m  - w_hat_m).^2

    r_hat, T_hat, b_hat, eta = bcd!(
    copy(r), copy(T), copy(b), y2, Q_inv, Ns, g_mag, Val{3*(Na+Ng)}(),Val{Na}(),
    1.0e-12, 100, 1.0e-8, 200, tol_dlogp, 1.0-14, 1; warn = false, w0 = w0.*2
    )
    w_hat_m = hcat(Array.(eta[:w])...)
    dw2 .+= (w_m.*2  - w_hat_m).^2

    # w_hat, w_dot_hat, s_hat  = dynamics(deepcopy(w0), y,  Q_inv, r, Ng,  Val{3*Nt}(), 1.0e-8, 200)
    # w_hat_m2 = hcat(Array.(w_hat)...)
    # dw2 .+= (w_m  - w_hat_m2).^2
end
mse = dw./N_mc
rmse_w_tot = sqrt.(mse)
println("zeta  1")
@show sum(mse)/length(mse)

mse2 = dw2./N_mc
rmse_w_tot2 = sqrt.(mse2)
println("zeta  2")
@show sum(mse2)/length(mse2)
if true
    fig, ax = subplots(1,1)
    crb_w =  sqrt.(crb(w, r, Q_inv))
    for i = 1:3
        ax.plot(w_norm .|> rad2deg, rmse_w_tot[i,:] .|> rad2deg, label = string(i));
        ax.plot(w_norm .|> rad2deg, crb_w[i,:] .|> rad2deg, label = "CRB $(i)");
    end
    ax.plot(w_norm .|> rad2deg, ones(size(w_norm))*sig_g/sqrt(Ng) .|> rad2deg);
    # ax.set_ylim(sig_g |>  rad2deg |> x -> x/10, sig_g |>  rad2deg |> x -> x*10)
    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.grid(true, which = "both")
    ax.legend()
    # ax.set_ylim(0.01, 0.1)
end

nothing


text:
using Viznet
using Test

@testset "chimera" begin
    lt = ChimeraLattice(3,3)
    @test length(lt) == 72 == length(CartesianIndices(lt)) == length(LinearIndices(lt))
    @test length(bonds(lt)) == 16*9 + (2*4*3)*2
end


text:
@testset "Aqua.jl" begin
  Aqua.test_all(Octavian, ambiguities=false, project_toml_formatting=false)
  @test isempty(Test.detect_ambiguities(Octavian))
end


text:
"""
# script run

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-08-20

Runs set of tests and benchmarks.
Commits and pushes results to Git.
# Examples

```jldoctest
julia>
```
"""

using DrWatson
@quickactivate "MEngProject"
using MEngProject,
    CUDA,
    DifferentialEquations,
    PyPlot,
    NNlib,
    ImageFiltering,
    Images,
    MEngProject.LaminartKernels,
    MEngProject.LaminartInitFunc,
    MEngProject.Utils,
    BenchmarkTools,
    Test

using OrdinaryDiffEq,
    ParameterizedFunctions, LSODA, Sundials, DiffEqDevTools, Noise

CUDA.allowscalar(false)

try include("test_illusions.jl") catch err; print(err) end
CUDA.reclaim()
try include("test_noise.jl") catch err; print(err) end
CUDA.reclaim()
try include("test_parameterVar.jl")catch err; print(err) end
CUDA.reclaim()
try include("test_kernelsVar.jl")catch err; print(err) end
CUDA.reclaim()
# try include("test_eqAtEquil.jl")catch err; print(err) end
CUDA.reclaim()

# try include("bench_solver.jl")catch err; print(err) end
try include("bench_imp.jl")catch err; print(err) end
CUDA.reclaim()
# try include("bench_eq.jl")catch err; print(err) end
try include("bench_kernels_dim.jl")catch err; print(err) end
CUDA.reclaim()
try include("bench_dim.jl")catch err; print(err) end
CUDA.reclaim()

run(`git add plots/*`)
run(`git commit -m "results"`)
run(`git pull`)
run(`git push`)


text:
using AlphaZero
using Test

include("../games/tictactoe/main.jl")
include("../games/connect-four/main.jl")
include("../scripts/lib/dummy_run.jl")

include("test_game.jl")

@testset "Testing Games" begin
  test_symmetries(Tictactoe.Game)
  test_symmetries(ConnectFour.Game)
end

@testset "Dummy Runs" begin
  @test dummy_run(Tictactoe)
  #@test dummy_run(ConnectFour) # Takes a bit too long for Travis
end