text:
export start_end, slide, nearest_start_end
const μs = 1.0e-6
const sμ = 1000000.0
"""

  start_end(C::SeisChannel)

Return start and endtimes of SeisChannel in DateTime format.
"""
function start_end(C::SeisChannel)
    return u2d.(SeisIO.t_win(C.t,C.fs) * 1e-6)
end

"""

  start_end(S::SeisData)

Return start and endtimes of SeisChannel in DateTime format.
"""
function start_end(S::SeisData)
  irr = falses(S.n)
  non_ts = findall(S.fs .== 0)
  irr[non_ts] .= true
  z = zero(Int64)
  start_times = Array{DateTime}(undef,S.n)
  end_times =  Array{DateTime}(undef,S.n)
  for i = 1:S.n
    start_times[i] = u2d(S.t[i][1,2] * μs)
    end_times[i] = u2d(S.t[i][end,2] + (irr[i] ? z : S.t[i][1,2] + sum(S.t[i][2:end,2]) +
        round(Int64, (length(S.x[i])-1)/(μs*S.fs[i]))) * μs)
        # ts data: start time in μs from epoch + sum of time gaps in μs + length of trace in μs
        # non-ts data: time of last sample
  end
  return start_times,end_times
end

"""
    slide(C::SeisChannel, cc_len::Int, cc_step::Int)

Cut `C` into equal length sliding windows.

# Arguments
- `C::SeisChannel`: SeisChannel.
- `cc_len::Int`: Cross-correlation window length [s].
- `cc_step::Int`: Step between cross-correlation windows [s].

# Returns
- `A::Array`: Array of sliding windows
- `starts::Array`: Array of start times of each window, in Unix time. E.g to convert
        Unix time to date time, use u2d(starts[1]) = 2018-08-12T00:00:00
- `ends::Array`: Array of end times of each window, in Unix time.
"""
function slide(C::SeisChannel, cc_len::Int, cc_step::Int)
  window_samples = Int(cc_len * C.fs)
  su,eu = SeisIO.t_win(C.t,C.fs) * μs
  starts = Array(range(su,stop=eu,step=cc_step))
  ends = starts .+ cc_len .- 1. / C.fs
  ind = findlast(x -> x <= eu,ends)
  starts = starts[1:ind]
  ends = ends[1:ind]

  # fill array with overlapping windows
  A = Array{eltype(C.x),2}(undef, window_samples,length(starts))
  s = convert.(Int,round.((hcat(starts,ends) .- su) .* C.fs .+ 1.))
  for ii = 1:length(starts)
    A[:,ii] .= @view(C.x[s[ii,1]:s[ii,2]])
  end
  return A,starts
end

"""
    nearest_start_end(C::SeisChannel, cc_len::Int, cc_step::Int)

Return best possible start, end times for data in `C` given the `cc_step` and `cc_len`.
"""
function nearest_start_end(C::SeisChannel, cc_len::Int, cc_step::Int)
  su,eu = SeisIO.t_win(C.t,C.fs) * μs
  ideal_start = d2u(DateTime(Date(u2d(su)))) # midnight of same day
  starts = Array(range(ideal_start,stop=eu,step=cc_step))
  ends = starts .+ cc_len .- 1. / C.fs
  return starts[findfirst(x -> x >= su, starts)], ends[findlast(x -> x <= eu,ends)]
end

"""
    nearest_start_end(S::DateTime,E::DateTime, cc_len::Int, cc_step::Int)

Return best possible start, end times for given starttime `S` and endtime `E`.
"""
function nearest_start_end(S::DateTime, E::DateTime, fs::Float64, cc_len::Int, cc_step::Int)
  ideal_start = DateTime(Date(S)) # midnight of same day
  starts = Array(ideal_start:Second(cc_step):endtime)
  ends = starts .+ Second(cc_len) .- Millisecond(convert(Int,1. / fs * 1e3))
  return starts[findfirst(x -> x >= S, starts)], ends[findlast(x -> x <= E,ends)]
end


text:
"""
Compare analyzers over one or more images.
"""
function compare(
    imgs::AbstractVector{<:AbstractArray},
    labels::AbstractVector{<:AbstractString},
    algs::AbstractVector{<:AbstractXAIMethod};
    backgroundcolor=:white, #:transparent
)
    n_algs = length(algs)
    n_imgs = length(imgs)
    fig = Figure(; backgroundcolor=backgroundcolor)

    # Add supertitle
    fig[0, :] = Label(fig, "ExplainableAI.jl"; textsize=24, color=(:black, 0.6))

    for (row, img) in enumerate(imgs)
        # Run first analyzer to get class probs
        output, _ = output_and_explain(img, algs[1])
        lab = Flux.onecold(output, labels)
        println(lab)

        fig[row, 1, Left()] = Label(fig, "asd"; halign=:right, padding=(0, 5, 0, 0))
        ax = Axis(fig[row, 1]; aspect=1)
        raw = image!(ax, RGB.(img))
        hidedecorations!(ax) # remove tics

        for (col, alg) in enumerate(algs)
            ax = Axis(fig[row, col + 1]; aspect=1)

            attr = explain(img, alg)
            hm = heatmap!(ax, attr; colormap=:bwr)
            hm.colorrange = (-1, 1)
            hidedecorations!(ax)
        end
    end
    cbar = fig[1:n_imgs, n_algs + 2] = Colorbar(fig, hm; label="Activations")

    # Add title for Input
    alg_name_padding = (0, 0, 5, 0)
    alg_name_rotation = pi / 8
    fig[1, 1, Top()] = Label(
        fig,
        "Input";
        halign=:left,
        valign=:bottom,
        padding=alg_name_padding,
        rotation=alg_name_rotation,
    )

    # Add titles for Algorithms
    for (col, alg) in enumerate(algs)
        fig[1, col + 1, Top()] = Label(
            fig,
            string(typeof(alg));
            halign=:left,
            valign=:bottom,
            padding=alg_name_padding,
            rotation=alg_name_rotation,
        )
    end

    # Trim gaps
    rowgap!(fig.layout, 7)
    colgap!(fig.layout, 7)
    colgap!(fig.layout, n_algs + 1, 25) # larger space before colorbar
    trim!(fig.layout)

    return fig
end

function compare(
    img::AbstractArray,
    labels::AbstractVector{<:AbstractString},
    algs::AbstractVector{<:AbstractXAIMethod};
    kwargs...,
)
    return compare([img], labels, algs, kwargs...)
end

ext: jl


text:
import Base.var
export var

function var(x::AffineExpr)
  return sum_squares(mean(x) - x) / (x.size[1] * x.size[2])
end


text:
function _c(lt, a, b)
    LI = LinearIndices(lt)
    isconnected(lt, LI[a...], LI[b...])
end

@i function isolve_largemem(out!::T, sg::AbstractSpinglass{MaskedSquareLattice}, reg::ArrayReg{1,TT}, REG_STACK; usecuda=false) where {T,TT<:Tropical{T}}
    @invcheckoff begin
    lt ← sg.lattice
    Lx ← size(lt, 1)
    Ly ← size(lt, 2)
    nbit ← Ly + 2
    Js ← sg.Js
    hs ← sg.hs
    k ← 0
    l ← 0
    for i=1:Lx
        @safe println("Layer $i/$Lx")
        for j=1:Ly-1
            if (_c(lt, (i,j), (i,j+1)), ~)
                INC(k)
                apply_Gvb!(reg, (@const (j, j+1)), Js[k], REG_STACK)
            end
        end
        for j=1:Ly
            if (lt.mask[i,j], ~)
                INC(l)
                apply_Gh!(reg, j, hs[l], REG_STACK)
            end
        end
        if (i!=Lx, ~)
            for j=1:Ly
                # store the information in qubit `j` to ancilla `nbit-j%2`
                if (j!=Ly && _c(lt, (i,j), (i+1,j+1)), ~)
                    apply_Gcp!(reg, (@const (j,nbit-j%2)), REG_STACK)
                end
                # interact with j-1 th qubit (a)
                if (j!=1 && _c(lt, (@const (i+1,j-1)), (i,j)), ~)
                    INC(k)
                    apply_Gvb!(reg, (@const (j-1, j)), Js[k], REG_STACK)
                end
                # onsite term (b)
                if (_c(lt, (i,j), (i+1,j)), ~)
                    INC(k)
                    apply_Ghb!(reg, j, Js[k], REG_STACK)
                else
                    apply_Gcut!(reg, j, REG_STACK)
                end
                if (j!=1 && _c(lt, (i,j-1), (i+1,j)), ~)
                    INC(k)
                    # interact with cached j-1 th qubit (c)
                    apply_Gvb!(reg, (@const (j,nbit-(j-1)%2)), Js[k], REG_STACK)
                    # erease the information in previous ancilla `nbit-(j-1)%2`
                    apply_Gcut!(reg, (@const nbit-(j-1)%2), REG_STACK)
                end
            end
        end
    end
    summed ← one(TT)
    i_sum(summed, reg.state)
    NiLang.SWAP(summed.n, out!)
    l → length(hs)
    k → length(Js)
    summed → one(TT)
    end
end

function cachesize_largemem(lt::MaskedSquareLattice)
    Lx = size(lt, 1)
    Ly = size(lt, 2)
    ncache = 0
    for i=1:Lx
        i!=Lx && for j=1:Ly
            if j!=Ly && _c(lt, (i,j), (i+1,j+1))
                ncache += 1
            end
            ncache += 1
            if j!=1 && _c(lt, (i,j-1), (i+1,j))
                ncache += 1
            end
        end
    end
    return ncache
end

@i function isolve(out!::T, sg::Spinglass{MaskedSquareLattice,T}, reg::ArrayReg{1,TT}, A_STACK, B_STACK; usecuda=false) where {T,TT<:Tropical{T}}
    lt ← sg.lattice
    Lx ← size(lt, 1)
    Ly ← size(lt, 2)
    nbit ← Ly + 2
    Js ← sg.Js
    hs ← sg.hs
    k ← 0
    l ← 0
    for i=1:Lx
        @safe println("Layer $i/$Lx")
        for j=1:Ly-1
            if (_c(lt, (i,j), (i,j+1)), ~)
                INC(k)
                apply_Gvb!(reg, (@const (j, j+1)), Js[k], A_STACK)
            end
        end
        for j=1:Ly
            if (lt.mask[i,j], ~)
                INC(l)
                apply_Gh!(reg, j, hs[l], A_STACK)
            end
        end
        if (i!=Lx, ~)
            @routine begin
                for j=1:Ly
                    # store the information in qubit `j` to ancilla `nbit-j%2`
                    if (j!=Ly && _c(lt, (i,j), (i+1,j+1)), ~)
                        apply_Gcp!(reg, (@const (j,nbit-j%2)), A_STACK)
                    end
                    # interact with j-1 th qubit (a)
                    if (j!=1 && _c(lt, (i+1,j-1), (i,j)), ~)
                        INC(k)
                        apply_Gvb!(reg, (@const (j-1, j)), Js[k], A_STACK)
                    end
                    # onsite term (b)
                    if (_c(lt, (i,j), (i+1,j)), ~)
                        INC(k)
                        apply_Ghb!(reg, j, Js[k], A_STACK)
                    else
                        apply_Gcut!(reg, j, A_STACK)
                    end
                    if (j!=1 && _c(lt, (i,j-1), (i+1,j)), ~)
                        INC(k)
                        # interact with cached j-1 th qubit (c)
                        apply_Gvb!(reg, (@const (j,nbit-(j-1)%2)), Js[k], A_STACK)
                        # erease the information in previous ancilla `nbit-(j-1)%2`
                        apply_Gcut!(reg, (@const nbit-(j-1)%2), A_STACK)
                    end
                end
            end
            incstack!(B_STACK)
            store_state!(B_STACK, reg.state)
            ~@routine
            # but wait, `k` should not be uncomputed
            for j=1:Ly
                if (j!=1 && _c(lt, (i+1,j-1), (i,j)), ~)
                    INC(k)
                end
                # onsite term (b)
                if (_c(lt, (i,j), (i+1,j)), ~)
                    INC(k)
                end
                if (j!=1 && _c(lt, (i,j-1), (i+1,j)), ~)
                    INC(k)
                end
            end
            # restore state
            swap_state!(B_STACK, reg.state)
        end
    end
    summed ← one(TT)
    i_sum(summed, reg.state)
    NiLang.SWAP(summed.n, out!)
    l → length(hs)
    k → length(Js)
    summed → one(TT)
end

cachesize_A(lt::MaskedSquareLattice) = size(lt,2)*3-1
cachesize_B(lt::MaskedSquareLattice) = size(lt,1)-1


text:
#An Extended Kalman Filter (EKF) example taken from Markus Gesmann's
#Mages' blog titled "Extended Kalman filter example in R" [1] which in turn
#was taken from a blog post titled "Fun with (Extended Kalman) Filters"
#by Dominic Steinitz [2].
#The web addresses are:
#[1] http://www.magesblog.com/2015/01/extended-kalman-filter-example-in-r.html
#[2] https://idontgetoutmuch.wordpress.com/2014/09/09/fun-with-extended-kalman-filters-4/

#Here we import the required modules
using StateSpace
using Distributions
using Gadfly
using DataFrames
using Colors

################################################################################
#Section: Generate noisy Observations
#-------------------------------------------------------------------------------

#Here we set the parameters with their true values
r = 0.2 #r is the growth rate
k = 100.0 #k is the carrying capacity
p0 = 0.1 * k # p0 is the initial population
Δt = 0.1 # Δt is the change in time.

#Define the logistic growth function to set the observations
function logisticGrowth(r, p, k, t)
    k * p * exp(r*t) / (k + p * (exp(r*t) - 1))
end
logisticGrowth(state) = logisticGrowth(state[1], state[2], k, Δt)

#Set the measurement noise variance
measurement_noise_variance = 25.0

#create the noisy observations (zero mean Gaussian noise)
numObs = 100
true_values = Vector{Float64}(numObs)
population_measurements = Vector{Float64}(numObs)
for i in 1:numObs
    true_values[i] = logisticGrowth(r, p0, k, i*Δt)
    population_measurements[i] = true_values[i] + randn() * sqrt(measurement_noise_variance)
end

#Since we're going to assume that our state consists of a growth rate, r,
#as well as a population rate, p, we need measurements of the growth rate.
#The problem is that we don't actually observe the growth rate directly.
#So we will set these measurements to zero.
growth_rate_measurements = zeros(numObs)

#Then we put these measurements together
measurements = [growth_rate_measurements population_measurements]'
#End Section: Generate noisy Observations
################################################################################

################################################################################
#Section: Describe Extended Kalman Filter parameters
#-------------------------------------------------------------------------------
#The state consists of the growth rate, r, and the population number, p. We'll
#assume that the growth rate is constant throughout time. And we'll assume that
#the population update follows the logistic growth pattern. So let's create that
#process function
function process_fcn(state)
    predict_growth_rate = state[1]
    predict_population = logisticGrowth(state)
    new_state = [predict_growth_rate, predict_population]
    return new_state
end
#Here we assume that there is no evolution noise so well create a zero matrix
process_noise_mat = diagm([0.001, 0.001])

#Now we need to describe the observation model. We'll assume that we don't
#observe the growth rate but we do observe the population. So the observation
#function is:
function observation_fcn(state)
    growth_rate_observation = 0.0
    population_observation = 1.0 * state[2]
    observation = [growth_rate_observation, population_observation]
    return observation
end

#Now we need to set the observation noise. We already set the measurement noise
#for the population earlier. We'lll arbitrarily make the variance for the
#growth rate (it doesn't matter because we don't explicity observe it):
observation_noise_mat = diagm([1.0, measurement_noise_variance])

#Create instance of our EKF model
nonLinSSM = NonlinearGaussianSSM(process_fcn, process_noise_mat, observation_fcn, observation_noise_mat)
#End Section: Describe Extended Kalman Filter parameters
################################################################################

################################################################################
#Section: Set initial guess of the state
#-------------------------------------------------------------------------------
initial_guess = MvNormal([0.5, 10], diagm([1.0,20.0]))
################################################################################

################################################################################
#Section: Execute the Extended Kalman Filter
#-------------------------------------------------------------------------------
filtered_state = filter(nonLinSSM, measurements, initial_guess)
#End Section: Execute  the Extended Kalman Filter
################################################################################

################################################################################
#Section: Plot Filtered results
#-------------------------------------------------------------------------------
#Here we are plotting the filtered results with Gadfly. See the Gadfly
#documentation for information about how plotting works if you are unfamiliar.
#Website: http://gadflyjl.org/

x_data = 1:numObs
population_array = Vector{Float64}(numObs+1)
confidence_array = Vector{Float64}(numObs+1)
population_array[1] = initial_guess.μ[2]
confidence_array[1] = 2*sqrt(initial_guess.Σ.mat[2,2])
for i in x_data
    current_state = filtered_state.state[i]
    population_array[i+1] = current_state.μ[2]
    confidence_array[i+1] = 2*sqrt(current_state.Σ.mat[2,2])
end
df_fs = DataFrame(
    x = [0;x_data],
    y = population_array,
    ymin = population_array - confidence_array,
    ymax = population_array + confidence_array,
    f = "Filtered values"
    )

n = 3
getColors = distinguishable_colors(n, Color[LCHab(70, 60, 240)],
                                   transform=c -> deuteranopic(c, 0.5),
                                   lchoices=Float64[65, 70, 75, 80],
                                   cchoices=Float64[0, 50, 60, 70],
                                   hchoices=linspace(0, 330, 24))
population_state_plot = plot(
    layer(x=0:numObs, y=[p0;measurements[2,:]'], Geom.point, Theme(default_color=getColors[2])),
    layer(x=0:numObs, y=[p0;true_values], Geom.line, Theme(default_color=getColors[3])),
    layer(df_fs, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon),
    Guide.xlabel("Measurement Number"), Guide.ylabel("Population"),
    Guide.manual_color_key("Colour Key",["Filtered Estimate", "Measurements","True Value "],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Extended Kalman Filter Example")
    )
display(population_state_plot)
#End Section: Plot Filtered results
################################################################################

################################################################################
#Section: Plot Filtered results
#-------------------------------------------------------------------------------
#Here we are plotting the filtered results with Gadfly. See the Gadfly
#documentation for information about how plotting works if you are unfamiliar.
#Website: http://gadflyjl.org/

x_data = 1:numObs
growth_rate_array = Vector{Float64}(numObs+1)
confidence_array = Vector{Float64}(numObs+1)
growth_rate_array[1] = initial_guess.μ[1]
confidence_array[1] = initial_guess.Σ.mat[1,1]
for i in x_data
    current_state = filtered_state.state[i]
    growth_rate_array[i+1] = current_state.μ[1]
    confidence_array[i+1] = 2*sqrt(current_state.Σ.mat[1,1])
end
df_fs = DataFrame(
    x = [0;x_data],
    y = growth_rate_array,
    ymin = growth_rate_array - confidence_array,
    ymax = growth_rate_array + confidence_array,
    f = "Filtered values"
    )

growth_rate_state_plot = plot(
    layer(x=0:numObs, y=ones(numObs+1)*r, Geom.line, Theme(default_color=getColors[3])),
    layer(df_fs, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon),
    Guide.xlabel("Measurement Number"), Guide.ylabel("Growth Rate"),
    Guide.manual_color_key("Colour Key",["Filtered Estimate", "Measurements","True Value "],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Extended Kalman Filter Example")
    )
display(growth_rate_state_plot)
#End Section: Plot Filtered results
################################################################################


text:
module Alpine

export
    # distances.jl
    sd,
    msd,
    rmsd,
    mag,
    wrapcoords!,
    unwrapcoords!,
    distancesq,
    distances!,

    # atoms.jl
    Frame,
    Atoms,
    Atom,
    add_property!,

    # atoms_analysis.jl
    group,
    box_dims,
    unscale,
    unscale!,
    unwrap!,

    # elements.jl
    element_property,
    element_properties,

    # files.jl
    AtomFileType,
    AtomFile,
    FrameHeader,
    read_atoms,
    read_frame_headers,
    write_atoms,
    LAMMPSTRJ, TRJ,
    LAMMPSDAT, DAT,
    XYZ,

    # potentials.jl
    vdw,
    stretching,
    bending,
    torsion,
    inversion,
    atom_type_number,
    atom_type_name,
    atom_type_params,

    # analysis.jl
    density,
    get_atoms,

    # graphs.jl
    Graph,
    vertices,
    edges,
    acyclic_path_table,
    acyclic_paths,
    unique_paths,
    generate_bonds,
    unique_upto_reversal,
    unique_neighborhoods,
    neighborhood_type_ids,

    # rdf.jl
    box_upper,
    box_lower,
    box_widths,
    ndf!,
    func_ndf,
    rdf_axis,
    split_by_types,
    sort_atoms_into_cells,
    rdf,
    rdf_components,
    rdf_components_all,
    rdf_axis

include("distances.jl")
include("atoms.jl")
include("atoms_analysis.jl")

include("elements/elements.jl")

include("files/files.jl")
include("files/lammpstrj.jl")
include("files/lammpsdat.jl")
include("files/xyz.jl")

#include("potentials/potentials.jl")
#include("potentials/uff/uff.jl")

include("analysis.jl")
include("graphs.jl")
include("rdf.jl")

end

max_stars_repo_name: NickMcNutt/Alpine

max_issues_repo_name: NickMcNutt/Alpine

max_forks_repo_name: NickMcNutt/Alpine


text:
models = [LogisticModel(), LinearModel(), QuantileModel()]

#thest that loss and gradients agree
n=50
p=10
x = rand(n, p)
y = round(rand(n))
for model in models
    mylearner=  GLMLearner(model, SimpleSGD(1.0, 1.0))
    update!(mylearner, x, y)
    @test isa(predict(mylearner, x), Vector)

    analytic_grad = grad!(mylearner.m, mylearner.gr, OnlineLearning.which_weights(mylearner.optimizer, mylearner.coefs), x, y)
    numeric_grad = Calculus.gradient(coefs -> loss(mylearner.m, predict!(mylearner.m, similar(y), coefs, x), y), mylearner.coefs)

    @test_approx_eq_eps analytic_grad numeric_grad 1e-8
end

# myadagrad = GLMLearner(LogisticModel(), AdaGrad(1.0))
# update!(myadagrad, x, y)

# myglmnet = GLMNetLearner(LogisticModel(), AdaGrad(1.0), 0.1, 0.1)

# update!(myglmnet, x, y)

# update!(myglmnet, sprand(10, 10, 0.3), round(rand(10)))


text:
using GLMakie

xs = range(0, 10, length = 25)
ys = range(0, 15, length = 25)

fig = Figure()

sl = Slider(fig[2, 1], range = 0:0.1:3, startvalue = 1)

zs = lift(sl.value) do s
    convert(Float64,s).*[cos(x) * sin(y) for x in xs, y in ys]
end

titleString = lift(sl.value) do s
    "Step: "*string(s)
end

ax = Axis(fig[1, 1], title = titleString)

hm=heatmap!(xs, ys, zs,colormap = Reverse(:Spectral))
Colorbar(fig[1, 2],hm,label = "Values")
fig

max_stars_repo_name: Kevin-Mattheus-Moerman/JuliaAdventures

max_issues_repo_name: Kevin-Mattheus-Moerman/JuliaAdventures

max_forks_repo_name: Kevin-Mattheus-Moerman/JuliaAdventures


text:
push!(LOAD_PATH, "../models")
using SuperResModels
using SparseInverseProblems

@everywhere begin
    include("../ufus_parameters.jl")
    parameters = SuperResModels.Conv2dParameters(x_max, x_max, filter, filter_dx, filter_dy, sigma, sigma, sigma, sigma)
    model_static = SuperResModels.Conv2d(parameters)
    dynamic_parameters = SuperResModels.DynamicConv2dParameters(K, tau, v_max, 20)
    model_dynamic = SuperResModels.DynamicConv2d(model_static, dynamic_parameters)
    n_x = model_static.n_x
    n_y = model_static.n_y
end

# Get time of script start
now_str = string(now())
now_str = replace(now_str, ":", "-")
now_str = replace(now_str, ".", "-")

# L2 norm single particle
thetas = reshape([x_max/2; x_max/2], 2, 1)
weights = [1.0]
single_particle_norm = norm(phi(model_static, thetas, weights), lp_norm)
@everywhere single_particle_norm = $single_particle_norm
println("norm = ", single_particle_norm)

# Test dynamic
println("Testing dynamic...")
particles_m = [x_max/4]
#particles_p = [3*x_max/4]
thetas = hcat([(x_max/2 - dx) * ones(1, length(particles_m)); particles_m'; 0.0; -v_max/2])
              #[(x_max/2 + dx) * ones(1, length(particles_p)); particles_p'; 0.0; v_max/2])
weights = ones(1)
target = phi(model_dynamic, thetas, weights)
function callback(old_thetas, thetas, weights, output, old_obj_val)
    #evalute current OV
    new_obj_val,t = SparseInverseProblems.loss(SparseInverseProblems.LSLoss(), output - target)
    #println("gap = $(old_obj_val - new_obj_val)")
    if old_obj_val - new_obj_val < 1E-4
        return true
    end
    return false
end
(thetas_est,weights_est) = SparseInverseProblems.ADCG(model_dynamic, SparseInverseProblems.LSLoss(), target, 2.0, callback=callback, max_iters=200)
println("theta = ", thetas, ", weights = ", weights)
println("theta_est = ", thetas_est, ", weights_est = ", weights_est)


text:
"""
nsolsc_solution_test

Test nsolsc with the atan function. Check answers and iteration stats.
"""
function nsolsc_solution_test()
    #
    # Local convergence with forward difference derivative
    #
    sdatal = nsolsc(atan, 1.0)
    solok = (abs(sdatal.solution) < 1.e-8)
    funok = (abs(sdatal.functionval) < 1.e-8)
    hs = size(sdatal.history)
    histok = (hs[1] == 5)
    locok = funok && solok && histok
    if ~locok
        println("local FD fails")
    end
    #
    # Local convergence with analytic derivative
    #
    sdataa = nsolsc(atan, 1.0, fpatan)
    solok = (abs(sdataa.solution) < 1.e-8)
    funok = (abs(sdataa.functionval) < 1.e-8)
    hs = size(sdataa.history)
    histok = (hs[1] == 5)
    analyticok = funok && solok && histok
    if ~analyticok
        println("failure with analytic derivative ")
        println(sdataa)
    end
    #
    # Global convergence
    #
    sdatag = nsolsc(atan, 10.0; maxit = 11, armfix = true)
    solok = (abs(sdatag.solution) < 1.e-8)
    funok = (abs(sdatag.functionval) < 1.e-8)
    hs = size(sdatag.history)
    histok = (hs[1] == 12)
    globok = funok && solok && histok
    if ~globok
        println("global FD fails")
    end
    #
    # Global convergence with parab3p
    #
    sdatap3p = nsolsc(atan, 30.0; rtol = 1.e-10, maxit = 11)
    solok = (abs(sdatap3p.solution) < 1.e-8)
    funok = (abs(sdatap3p.functionval) < 1.e-8)
    hs = size(sdatap3p.history)
    histok = (hs[1] == 12)
    p3pok = funok && solok && histok
    if ~p3pok
        println("parab3p fails")
    end

    #
    # Local convergence with secant method
    #
    #sdatas=nsolsc(atan,1.0; solver="secant")
    sdatas = secant(atan, 1.0)
    solok = (abs(sdatas.solution) < 1.e-10)
    funok = (abs(sdatas.functionval) < 1.e-10)
    hs = size(sdatas.history)
    histok = (hs[1] == 6)
    secantok = funok && solok && histok
    if ~secantok
        println("secant failure")
    end
    #
    # Initialize secant method when x0=0
    #
    #zedata=nsolsc(fcos,0.0;solver="secant",rtol=1.e-9)
    zedata = secant(fcos, 0.0; solver = "secant", rtol = 1.e-9)
    solution = 7.390851333858823e-01
    solok = (abs(zedata.solution - solution) < 1.e-9)
    funok = (abs(zedata.functionval) < 1.e-9)
    hs = size(zedata.history)
    histok = (hs[1] == 7)
    zecok = funok && solok && histok
    if ~zecok
        println("local FD fixup at zero fails")
    end
    #
    # Tricky line search problem
    #
    sdatal = nsolsc(
        linatan,
        200.0;
        sham = 5,
        maxit = 20,
        armmax = 10,
        armfix = true,
        rtol = 1.e-10,
    )
    solution = -100.0
    solok = (abs(sdatal.solution - solution) < 1.e-8)
    funok = (abs(sdatal.functionval) < 1.e-8)
    hs = size(sdatal.history)
    histok = (hs[1] == 6)
    shamfastok = funok && solok && histok
    if ~shamfastok
        println("Fast Shamanskii response FAILURE")
    end
    #
    # Test linesearch failure complaints.
    #
    armfail = nsolsc(atan, 10.0; armmax = 1, armfix = true, printerr = false)
    afok = false
    if armfail.idid == false && armfail.errcode == 1
        afok = true
    else
        println("Armijo failure test FAILED.")
    end
    #
    # Test residual failure mode and no history.
    #
    resok = false
    resfail = nsolsc(atan, 10.0; maxit = 3, armfix = true, keepsolhist = false)
    if resfail.idid == false && resfail.errcode == 10
        resok = true
    else
        println("Residual failure test FAILED.")
    end
    #
    # Test stagnation mode
    #
    stagdatan =
        nsolsc(ftanx, 4.5, ftanxp; rtol = 1.e-17, atol = 1.e-17, armfix = true, maxit = 14)
    fvals = stagdatan.history
    avals = stagdatan.stats.iarm
    ifvals = stagdatan.stats.ifun
    jvals = stagdatan.stats.ijac
    stagl = (length(fvals) == 6)
    stagf = (fvals[5] < 1.e-15)
    stags = (avals[6] == 5) && (ifvals[6] == 6) && (jvals[6] == 1)
    stagok = stagl && stags && stagf
    if ~stagok
        println("Stagnation test FAILED")
    end
    #
    #
    # Test chord method
    #
    lttest = nsolsc(atan, 0.5; solver = "chord")
    fvals = lttest.history
    chordl = (length(fvals) == 11)
    ratl = fvals[11] / fvals[10]
    chordr = (abs(ratl - 0.25) < 1.e-7)
    solok = (fvals[11] < 1.e-6)
    chordok = chordl && chordr && solok
    if ~chordok
        println("Chord test failed")
    end
    #
    # Make sure nsolsc knows if the solution and the intial iterate are
    # the same
    #
    lotout = nsolsc(flot, 0.0)
    lotok = (lotout.errcode == -1)
    if ~lotok
        println("Lottery test failed")
    end
    return locok &&
           globok &&
           secantok &&
           analyticok &&
           zecok &&
           shamfastok &&
           afok &&
           resok &&
           p3pok &&
           chordok &&
           lotok
end

function flot(x)
    flot = x * exp(x)
end


text:
n = [
  37107287533902102798797998220837590246510135740250
  46376937677490009712648124896970078050417018260538
  74324986199524741059474233309513058123726617309629
  91942213363574161572522430563301811072406154908250
  23067588207539346171171980310421047513778063246676
  89261670696623633820136378418383684178734361726757
  28112879812849979408065481931592621691275889832738
  44274228917432520321923589422876796487670272189318
  47451445736001306439091167216856844588711603153276
  70386486105843025439939619828917593665686757934951
  62176457141856560629502157223196586755079324193331
  64906352462741904929101432445813822663347944758178
  92575867718337217661963751590579239728245598838407
  58203565325359399008402633568948830189458628227828
  80181199384826282014278194139940567587151170094390
  35398664372827112653829987240784473053190104293586
  86515506006295864861532075273371959191420517255829
  71693888707715466499115593487603532921714970056938
  54370070576826684624621495650076471787294438377604
  53282654108756828443191190634694037855217779295145
  36123272525000296071075082563815656710885258350721
  45876576172410976447339110607218265236877223636045
  17423706905851860660448207621209813287860733969412
  81142660418086830619328460811191061556940512689692
  51934325451728388641918047049293215058642563049483
  62467221648435076201727918039944693004732956340691
  15732444386908125794514089057706229429197107928209
  55037687525678773091862540744969844508330393682126
  18336384825330154686196124348767681297534375946515
  80386287592878490201521685554828717201219257766954
  78182833757993103614740356856449095527097864797581
  16726320100436897842553539920931837441497806860984
  48403098129077791799088218795327364475675590848030
  87086987551392711854517078544161852424320693150332
  59959406895756536782107074926966537676326235447210
  69793950679652694742597709739166693763042633987085
  41052684708299085211399427365734116182760315001271
  65378607361501080857009149939512557028198746004375
  35829035317434717326932123578154982629742552737307
  94953759765105305946966067683156574377167401875275
  88902802571733229619176668713819931811048770190271
  25267680276078003013678680992525463401061632866526
  36270218540497705585629946580636237993140746255962
  24074486908231174977792365466257246923322810917141
  91430288197103288597806669760892938638285025333403
  34413065578016127815921815005561868836468420090470
  23053081172816430487623791969842487255036638784583
  11487696932154902810424020138335124462181441773470
  63783299490636259666498587618221225225512486764533
  67720186971698544312419572409913959008952310058822
  95548255300263520781532296796249481641953868218774
  76085327132285723110424803456124867697064507995236
  37774242535411291684276865538926205024910326572967
  23701913275725675285653248258265463092207058596522
  29798860272258331913126375147341994889534765745501
  18495701454879288984856827726077713721403798879715
  38298203783031473527721580348144513491373226651381
  34829543829199918180278916522431027392251122869539
  40957953066405232632538044100059654939159879593635
  29746152185502371307642255121183693803580388584903
  41698116222072977186158236678424689157993532961922
  62467957194401269043877107275048102390895523597457
  23189706772547915061505504953922979530901129967519
  86188088225875314529584099251203829009407770775672
  11306739708304724483816533873502340845647058077308
  82959174767140363198008187129011875491310547126581
  97623331044818386269515456334926366572897563400500
  42846280183517070527831839425882145521227251250327
  55121603546981200581762165212827652751691296897789
  32238195734329339946437501907836945765883352399886
  75506164965184775180738168837861091527357929701337
  62177842752192623401942399639168044983993173312731
  32924185707147349566916674687634660915035914677504
  99518671430235219628894890102423325116913619626622
  73267460800591547471830798392868535206946944540724
  76841822524674417161514036427982273348055556214818
  97142617910342598647204516893989422179826088076852
  87783646182799346313767754307809363333018982642090
  10848802521674670883215120185883543223812876952786
  71329612474782464538636993009049310363619763878039
  62184073572399794223406235393808339651327408011116
  66627891981488087797941876876144230030984490851411
  60661826293682836764744779239180335110989069790714
  85786944089552990653640447425576083659976645795096
  66024396409905389607120198219976047599490197230297
  64913982680032973156037120041377903785566085089252
  16730939319872750275468906903707539413042652315011
  94809377245048795150954100921645863754710598436791
  78639167021187492431995700641917969777599028300699
  15368713711936614952811305876380278410754449733078
  40789923115535562561142322423255033685442488917353
  44889911501440648020369068063960672322193204149535
  41503128880339536053299340368006977710650566631954
  81234880673210146739058568557934581403627822703280
  82616570773948327592232845941706525094512325230608
  22918802058777319719839450180888072429661980811197
  77158542502016545090413245809786882778948721859617
  72107838435069186155435662884062257473692284509516
  20849603980134001723930671666823555245252804609722
  53503534226472524250874054075591789781264330331690
]

println(string(sum(n))[1:10])

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl


text:
#!/usr/bin/env julia

import Celeste.GalsimBenchmark
import Celeste.DeterministicVI: infer_source
import Celeste.DeterministicVIImagePSF: infer_source_fft

test_case_names = String[]
infer_source_callback = infer_source
for arg in ARGS
    if arg == "--fft"
        infer_source_callback = infer_source_fft
    elseif ismatch(r"^--test_case_names", arg)
        # test case name: function name in test_case_definitions.py, or CL_DESCR
        # field from FITS header
        test_case_names =
            split(replace(arg, "--test_case_names=", ""), ",")
    else
        println("Usage: run_galsim_benchmark.jl " *
                "[--fft] [--test_case_names=test1,test2,...]")
        error("Invalid argument: ", arg)
    end
end

GalsimBenchmark.main(
    test_case_names=test_case_names,
    infer_source_callback=infer_source_callback)


text:
# misc CCI-file loaders

# function load!(ds::DataSet{IVData, CCILoader}, dem, glaciermask, target_proj, outline, cropbox,
#                pp::Phys)::Tuple{SomeData,Matrix{Bool}}
#     # TODO: allow other IV sets
#     @unpack files, opts, daterange = ds

#     zone = split(split(target_proj)[2],'=')[2]
#     flx = files[1]*"utm$zone-x.asc.gz"
#     fly = files[1]*"utm$zone-y.asc.gz"

#     dsx = reconstruct(ds, files=[flx])
#     ivx = load_asciigrid!(dsx, :ivx, cropbox)
#     dsy = reconstruct(ds, files=[fly])
#     ivy = load_asciigrid!(dsy, :ivy, cropbox)
#     fill = ivx.v.==FILL
#     tmp = sqrt.(ivx.v.^2 + ivy.v.^2)
#     # m/d -> m/y
#     tmp .*= 365
#     tmp[fill] = FILL
#     iv = reconstruct(ivx, v=tmp)

#     # # check projection -> now using re-projected files
#     # proj = ""
#     # pfl = splitext(splitext(files[1])[1])[1]*".prj"
#     # if isfile(pfl)
#     #     for st in split(readstring(`gdalsrsinfo iv-2017-y.prj`), '\n')
#     #         if startswith(st, "PROJ.4 :")
#     #             proj = strip(split(st, '\'')[2])
#     #             @assert startswith(proj, "+proj=")
#     #         end
#     #     end
#     #     if proj!=target_proj

#     # end

#     # according to *.xml sheet
#     # <meanE units="m/d">0.00</meanE>
#     # <rmseE units="m/d">0.01</rmseE>
#     # <meanN units="m/d">0.00</meanN>
#     # <rmseN units="m/d">0.02</rmseN>
#     @get sigma = opts | sqrt((0.01*365)^2 + (0.02*365)^2)

#     # get/make masks
#     ivmask = get_cache!(ds, outline.gid, :ivmask,
#                         () -> maskit(iv, outline)[1],
#                         Bool, size(iv))

#     # fill in any gaps inside the glacier
#     fillgaps!(iv, ivmask, mask_unfilled_holes=true)
#     return SomeData(vals=iv, sigma=sigma), ivmask
# end

function load!(ds::DataSet{IVData, CCILoader}, dem, glaciermask, target_proj, outline, cropbox,
               pp::Phys)::Tuple{SomeData,Matrix{Bool}}
    @unpack files, opts, daterange = ds
    # figure out utm-zone
    zone = split(split(target_proj)[2],'=')[2]
    fln = files[1]*"-utm$zone.tif"
    ds = reconstruct(ds, files=[fln])
    iv = load_geotiff!(ds, Symbol("iv-$zone"), cropbox)

    @assert !any(isnan.(iv.v)) "Loaded IV has NaNs"

    # get/make masks
    ivmask = get_cache!(ds, outline.gid, :ivmask,
                        () -> maskit(iv, outline)[1],
                        Bool, size(iv))

    # I don't think this is a good idea for the rough data we got # fill in any gaps inside the glacier
    # fillgaps!(iv, ivmask, mask_unfilled_holes=true)

    # mask it where it is within threshold
    th = get(opts, :mask_threshold, (0.0, Inf))
    ivmask[.!(th[1].< iv.v .<th[2])] = false
    # mask where FILL value
    ivmask[iv.v.==FILL] = false

    if any(size(iv.v).<2)
        # if there is not enough IV data, don't use it
        return SomeData(nothing), Matrix{Bool}(0,0)
    else
        return SomeData(vals=iv, sigma=opts[:sigma]), ivmask
    end
end

function load!(ds::DataSet{DhdtData, CCILoader},
               bdot, terminus_flux,
               dem,
               glaciermask, target_proj, outline, cropbox,
               iscalving,
               pp::Phys)::SomeData
    @unpack files, opts, daterange = ds
    zone = split(split(target_proj)[2],'=')[2]
    fln = files[1]*"-utm$zone.tif"
    ds = reconstruct(ds, files=[fln])
    dhdt = SomeData(vals=load_geotiff!(ds, Symbol("dhdt-$zone"), cropbox),
                    sigma=ds.opts[:sigma], vals_prior=0.0)
    if length(dhdt.vals.v)==0
        # use SyntheticBenchLoader
        ds = DataSet{DhdtData, SyntheticBenchLoader}(ds)
        dhdt = load!(ds,bdot, terminus_flux,
                     dem,
                     glaciermask, target_proj, outline, cropbox,
                     iscalving,
                     pp)
    end
    return dhdt
end


text:
import QuadGK

mutable struct Integrate{V} <: State{V}
    value::V
end

Integrate(; unit, _type, _...) = begin
    V = valuetype(_type, value(unit))
    v = zero(V)
    Integrate{V}(v)
end

constructortags(::Val{:Integrate}) = (:unit,)

updatetags!(d, ::Val{:Integrate}; _...) = begin
    !haskey(d, :from) && (d[:from] = @q zero($(d[:_type])))
    !haskey(d, :to) && (d[:to] = @q zero($(d[:_type])))
end

genvartype(v::VarInfo, ::Val{:Integrate}; V, _...) = @q Integrate{$V}

gendefault(v::VarInfo, ::Val{:Integrate}) = nothing

genupdate(v::VarInfo, ::Val{:Integrate}, ::MainStep; kw...) = begin
    kwarg = only(v.kwargs)
    u = extractfunckwargtuple(kwarg)[3]
    @gensym s a b f i
    @q let $s = $(symstate(v)),
           $a = $C.unitfy($C.value($(gettag(v, :from))), $u),
           $b = $C.unitfy($C.value($(gettag(v, :to))), $u),
           $f = $(genfunc(v; unitfy=false))
        $i = $C.QuadGK.quadgk($f, $a, $b) |> first
        $C.store!($s, $i)
    end
end


text:
function update_remote_backup_cli(argv=ARGS)

    # read ARGS
    argset = ArgParse.ArgParseSettings()
    ArgParse.@add_arg_table! argset begin
        "--force", "-f"
            help = "The search pattern (It will be evaluated as a julia expression)"
            action = :store_true
    end

    parsed_args = ArgParse.parse_args(argv, argset)
    force = parsed_args["force"]

    # update
    update_remote_backup(;force, loglevel = :verbose)

    return nothing
end


text:
using RNLA, Test, LinearAlgebra

@testset "rsvd_test1" begin
    X = [1 2 3 4; 5 6 7 8]
    u1, s1, v1 = svd(X)
    u2, s2, v2 = rsvd(X, 2)

    @test norm(abs.(u1) - abs.(u2)) < sqrt(eps())
    @test norm(abs.(v1) - abs.(v2)) < sqrt(eps())
    @test norm(abs.(s1) - abs.(s2)) < sqrt(eps())
end

@testset "rsvd_test2" begin
    X = [1 2 3; 4 5 6; 7 8 9]
    u1, s1, v1 = svd(X)
    u2, s2, v2 = rsvd(X, 3)

    @test norm(abs.(u1) - abs.(u2)) < sqrt(eps())
    @test norm(abs.(v1) - abs.(v2)) < sqrt(eps())
    @test norm(abs.(s1) - abs.(s2)) < sqrt(eps())
end

@testset "rsvd_test3" begin
    X = [1 2; 3 4; 5 6; 7 8; 9 10; 11 12]
    u1, s1, v1 = svd(X)
    u2, s2, v2 = rsvd(X, 2)

    @test norm(abs.(u1) - abs.(u2)) < sqrt(eps())
    @test norm(abs.(v1) - abs.(v2)) < sqrt(eps())
    @test norm(abs.(s1) - abs.(s2)) < sqrt(eps())
end

@testset "rsvd_test4" begin
    m = 1000
    k = 500
    X = rand(m, m) * rand(m, k)

    u1, s1, v1 = svd(X)
    u2, s2, v2 = rsvd(X, k)

    @test norm(abs.(u1) - abs.(u2)) < sqrt(eps())
    @test norm(abs.(v1) - abs.(v2)) < sqrt(eps())
    @test norm(abs.(s1) - abs.(s2)) < sqrt(eps())
end


text:
module Dif
    export mxdiflg

    function mxdiflg(a1::Array, a2::Array)
        try
            min1, max1 = extrema( map( length, a1 ) )
            min2, max2 = extrema( map( length, a2 ) )
            max( abs(min1-max2), abs(max1-min2) )
        catch
            -1
        end
    end

end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using BackpropNeuralNet
using Base.Test

# write your own tests here
include("test_neural.jl")


text:
function probit_reg(y, X, Xfull, nmcmc)

	N = size(Xfull)[1]
	n = length(y)

	beta_save = fill(0.0, 2, nmcmc)
	p_save = fill(0.0, N, nmcmc)
	y_save = BitArray(undef, N, nmcmc)

	# Priors
	Σ_β = ScalMat(2, 2.25)
	Σ_β_inv = inv(Σ_β)
	μ_β = fill(0.0, 2)

	# Initial values
	β = μ_β

	Σtemp = PDMat(Symmetric(inv(Σ_β_inv + X' * X)))

	# Gibbs loop
	for k in 1:nmcmc

		# Sample z
		z = fill(0.0, length(y))
		μ = X * β

		z[y .== 1] = rand.(Truncated.(Normal.(μ[y .== 1], 1.0), 0.0, Inf))
		z[y .== 0] = rand.(Truncated.(Normal.(μ[y .== 0], 1.0), -Inf, 0.0))

		# Sample beta
		βtemp = Σtemp * (Σ_β_inv * μ_β + X' * z)

		β = rand(MvNormal(βtemp, Σtemp))

		beta_save[:, k] = β
		p_save[:, k] = cdf.(Normal(), Xfull * β)
		y_save[:, k] = rand.(Bernoulli.(p_save[:, k]))

	end

	return (β = beta_save, p = p_save, y = y_save)
end

function probit_pprb(ynew, loc, p_prop, β_prop, nmcmc)

	nprop = size(p_prop)[2]

	β_save = fill(0.0, 2, nmcmc)
	y_save = BitArray(undef, N, nmcmc)
	ynew_save = fill(0, nmcmc)

	init = sample(1:nprop, 1)
	p = p_prop[:, init]

	β = β_prop[:, init]
	β_save[:, 1] = β

	@inbounds for k in 1:nmcmc

		idx = sample(1:nprop, 1)
		pstar = p_prop[:, idx]

		mh1 = logpdf(Bernoulli(pstar[loc]), ynew)
		mh2 = logpdf(Bernoulli(p[loc]), ynew)

		r = exp(mh1 - mh2)
		if r > rand(Uniform())
			p = pstar
			β = β_prop[:, idx]
		end

		β_save[:, k] = β
		# p_save[:, k] = p
		y_save[:, k] = rand.(Bernoulli.(p))
		ynew_save[k] = ynew

	end

	return (y = y_save, ynew = ynew_save, β = β_save)

end

max_stars_repo_name: clint-leach/optimal-design

max_issues_repo_name: clint-leach/optimal-design

max_forks_repo_name: clint-leach/optimal-design


text:
"""
  shrink(v, c)

Soft-threshold operator. Returns sign(v)⋅max(0, |v|-c)
"""
shrink(v::T, c::T) where {T<:AbstractFloat} = v > c ? v - c : (v < -c ? v + c : zero(T))

"""
  out = max(1-c/|x|, 0) ⋅ x
"""
function shrinkL2!(out::M, x::M, c::T) where {M <: AbstractVecOrMat{T}} where T
  tmp = max(one(T) - c / norm(x), zero(T))
  tmp > zero(T) ? rmul!(copyto!(out, x), tmp) : fill!(out, zero(T))
end

# """
# Computes a proximal operation for the penalty
# \[ \lambda_1\cdot(|u|+|v|) + \lambda_2\cdot|u-v| . \]
#
# Returns (u, v) = arg min (1/2γ) * ((u - x)^2 + (v-y)^2) + λ1 * (|u|+|v|) + λ2 * |u-v|
# """
function proxL1Fused(x::T, y::T, λ1::T, λ2::T, γ::T=one(T)) where {T<:AbstractFloat}

  if λ2 > zero(T)
    t = λ2 * γ
    if x > y + 2. * t
      u = x - t
      v = y + t
    elseif y > x + 2. * t
      u = x + t
      v = y - t
    else
      u = v = (x + y) / 2.
    end
  else
    u = x
    v = y
  end

  t = λ1 * γ
  (shrink(u, t), shrink(v, t))
end

####################################
#
# helper functions
#
####################################

function A_mul_B_row(A::AbstractMatrix{T}, b::AbstractVector{T}, row::Int64) where {T<:AbstractFloat}
  n, p = size(A)
  ((p == length(b)) && (1 <= row <= n)) || throw(DimensionMismatch())

  v = zero(T)
  for i=1:p
    @inbounds v += A[row, i] * b[i]
  end
  v
end

function A_mul_B_row(A::AbstractMatrix{T}, b::SparseVector{T}, row::Int64) where {T<:AbstractFloat}
  n, p = size(A)
  ((p == length(b)) && (1 <= row <= n)) || throw(DimensionMismatch())

  nzval = SparseArrays.nonzeros(b)
  rowval = SparseArrays.nonzeroinds(b)
  v = zero(T)
  for i=1:length(nzval)
    @inbounds v += A[row, rowval[i]] * nzval[i]
  end
  v
end

function A_mul_B_row(A::AbstractMatrix{T}, b::SparseIterate{T, 1}, row::Int64) where {T<:AbstractFloat}
  n, p = size(A)
  ((p == length(b)) && (1 <= row <= n)) || throw(DimensionMismatch())

  v = zero(T)
  @inbounds for i = 1:nnz(b)
      v += A[row, b.nzval2ind[i]] * b.nzval[i]
  end
  v
end

function At_mul_B_row(A::AbstractMatrix{T}, b::StridedVector{T}, row::Int64) where {T<:AbstractFloat}
  n, p = size(A)
  ((n == length(b)) && (1 <= row <= p)) || throw(DimensionMismatch())

  v = zero(T)
  for i=1:n
    @inbounds v += A[i, row] * b[i]
  end
  v
end

function At_mul_B_row(A::AbstractMatrix{T}, b::SparseVector{T}, row::Int64) where {T<:AbstractFloat}
  n, p = size(A)
  ((n == length(b)) && (1 <= row <= p)) || throw(DimensionMismatch())

  nzval = SparseArrays.nonzeros(b)
  rowval = SparseArrays.nonzeroinds(b)
  v = zero(T)
  for i=1:length(nzval)
    @inbounds v += A[rowval[i], row] * nzval[i]
  end
  v
end

function At_mul_B_row(A::AbstractMatrix{T}, b::SparseIterate{T, 1}, row::Int64) where {T<:AbstractFloat}
  n, p = size(A)
  ((n == length(b)) && (1 <= row <= n)) || throw(DimensionMismatch())

  v = zero(T)
  @inbounds for i = 1:b.nnz
      v += A[b.nzval2ind[i], row] * b.nzval[i]
  end
  v
end

function A_mul_X_mul_B_rc(
  A::Symmetric{T},
  X::SymmetricSparseIterate{T},
  B::Symmetric{T},
  r::Int,
  c::Int
  ) where {T<:AbstractFloat}

  p = X.p
  v = zero(T)
  for j=1:nnz(X)
    ind = X.nzval2ind[j]
    ri, ci = ind2subLowerTriangular(p, ind)
    if ri == ci
      @inbounds v += A[ri, r] * B[ci, c] * X.nzval[j]
    else
      @inbounds v += (A[ri, r] * B[ci, c] + A[ci, r] * B[ri, c]) * X.nzval[j]
    end
  end
  v
end

function A_mul_X_mul_B_rc(
  A::Symmetric{T},
  X::SparseIterate{T, 1},
  B::Symmetric{T},
  r::Int,
  c::Int
  ) where {T<:AbstractFloat}

  v = zero(T)
  I = CartesianIndices((size(A, 2), size(B, 1)))
  for j=1:nnz(X)
    ind = X.nzval2ind[j]
    @inbounds v += A[I[ind][1], r] * B[I[ind][2], c] * X.nzval[j]
  end
  v
end

function A_mul_X_mul_B(
  A::Symmetric{T},
  X::Union{SymmetricSparseIterate{T},SparseIterate{T}},
  B::Symmetric{T}
  ) where {T<:AbstractFloat}

  p = size(A, 1)
  out = zeros(T, p, p)

  for c=1:p, r=1:p
    @inbounds out[r,c] = A_mul_X_mul_B_rc(A, X, B, r, c)
  end
  out
end

function A_mul_UVt_mul_B_rc(
  A::Symmetric{T},
  U::StridedMatrix{T},
  V::StridedMatrix{T},
  B::Symmetric{T},
  r::Int,
  c::Int
  ) where {T<:AbstractFloat}
  # check input dimensions

  nr, nc = size(U)
  v = zero(T)
  for j=1:nc
    v1 = zero(T)   # stores A[r, :]*U[:,j]
    v2 = zero(T)   # stores B[c, :]*U[:,j]
    for k=1:nr
      @inbounds v1 += A[k, r] * U[k, j]
      @inbounds v2 += B[k, c] * V[k, j]
    end
    v += v1 * v2
  end
  v
end

function A_mul_UVt_mul_B(
    A::Symmetric{T},
    U::StridedMatrix{T},
    V::StridedMatrix{T},
    B::Symmetric{T}
  ) where {T<:AbstractFloat}

  # check input dimensions

  p = size(A, 1)
  out = zeros(T, p, p)

  for c=1:p, r=1:p
    @inbounds out[r,c] = A_mul_UVt_mul_B_rc(A, U, V, B, r, c)
  end
  out
end

function A_mul_UUt_mul_B_rc(
  A::Symmetric{T},
  U::StridedMatrix{T},
  B::Symmetric{T},
  r::Int,
  c::Int
  ) where {T<:AbstractFloat}
  # check input dimensions

  nr, nc = size(U)
  v = zero(T)
  for j=1:nc
    v1 = zero(T)   # stores A[r, :]*U[:,j]
    v2 = zero(T)   # stores B[c, :]*U[:,j]
    for k=1:nr
      @inbounds v1 += A[k, r] * U[k, j]
      @inbounds v2 += B[k, c] * U[k, j]
    end
    v += v1 * v2
  end
  v
end

function A_mul_UUt_mul_B(A::Symmetric{T},
  U::StridedMatrix{T},
  B::Symmetric{T}
  ) where {T<:AbstractFloat}

  # check input dimensions

  p = size(A, 1)
  out = zeros(T, p, p)

  for c=1:p, r=1:p
    @inbounds out[r,c] = A_mul_UUt_mul_B_rc(A, U, B, r, c)
  end
  out
end

#####
#
# functions to operate with sparse lower triangular
#

function ind2subLowerTriangular(p::T, ind::T) where {T<:Integer}
  rvLinear = div(p*(p+1), 2) - ind
  k = trunc(T, (sqrt(1+8*rvLinear)-1.)/2. )
  j = rvLinear - div(k*(k+1), 2)
  (p-j, p-k)
end

sub2indLowerTriangular(p::T, r::T, c::T) where {T<:Integer} = p*(c-1)-div(c*(c-1),2)+r

function vec2tril(x::SparseVector, p::Int64)
  nx = nnz(x)
  nzval = SparseArrays.nonzeros(x)
  nzind = SparseArrays.nonzeroinds(x)

  I = zeros(Int64, nx)
  J = zeros(Int64, nx)
  for i=1:nx
    I[i], J[i] = ind2subLowerTriangular(p, nzind[i])
  end

  sparse(I,J,nzval, p, p)
end

function vec2tril(x::SparseIterate, p::Int64)
  nx = nnz(x)

  I = zeros(Int64, nx)
  J = zeros(Int64, nx)
  for i=1:nx
    I[i], J[i] = ind2subLowerTriangular(p, x.nzval2full[i])
  end

  sparse(I,J,x.nzval[1:nx], p, p)
end

function tril2symmetric(Δ::SparseMatrixCSC)
  lDelta = tril(Δ, -1)
  dDelta = spdiagm(diag(Δ))
  (lDelta + lDelta') + dDelta
end

#######

function norm_diff(A::AbstractArray{T}, B::AbstractArray{T}, p::Real=2.) where {T<:AbstractFloat}
  size(A) == size(B) || throw(DimensionMismatch())

  v = zero(T)
  if p == 2.
    @inbounds @simd for i in eachindex(A)
      t =
      v += abs2( A[i] - B[i] )
    end
    return sqrt(v)
 elseif p == Inf
   @inbounds @simd for i in eachindex(A)
     t = abs( A[i] - B[i] )
     if t > v
       v = t
     end
   end
   return v
 else
   throw(ArgumentError("p should be 2 or Inf"))
 end
end


text:
# # Specifying and fine-tuning prior distributions
#
# This example shows how to specify priors in `ParameterEstimocean`
# and illustrates some of their properties.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add ParameterEstimocean, Oceananigans, Distributions, CairoMakie"
# ```
#
# Then we import packages

using ParameterEstimocean
using CairoMakie
using Distributions

# # What's a prior?
#
# A "prior" is a probability distribution for an known parameter that
# represents "prior knowledge" of the problem at hand, before we generate
# new knowledge with a new model-data comparison via ensemble Kalman inversion (EKI).
#
# In the (EKI) context, a prior distribution determines (i) the initial distribution
# of parameter values in the ensemble and (ii) the mapping between physical space
# where constraints on parameter values are applied, and the "unconstrained" space
# in which the EKI dynamics evolve (and where parameters are assumed normally-distributed).
#
# # Prior flavors
#
# `ParameterEstimocean` supports three types of prior
# distributions. The normal prior distribution is provided by `Distributions.jl`,

normal_prior = Normal(0.1, 0.1)

# For `LogNormal` priors, we provide a special constructor that takes the `mean` and
# standard deviation `std` of the prior distribution in _constrained_ space:

lognormal_prior = lognormal(mean=0.1, std=0.1)

# `lognormal` calculates `μ` and `σ` for `Distributions.Lognormal(μ, σ)`.
# Parameters with a `Lognormal` prior distribution are _strictly_ positive.
# Finally, the "scaled" logit-normal prior,

logitnormal_prior = ScaledLogitNormal(bounds=(0, 0.2))

# can be used to provide strict upper and lower bounds for parameter values.
#
# ## Samples
#
# Sampling the distributions reveals their properties. We define a function to use
# later, too.

function visualize_priors(priors; samples = 10^6)
    fig = Figure()
    ax = Axis(fig[1, 1:9], xlabel="Random samples of priors", ylabel="Density")

    for (prior, name) in zip(values(priors), keys(priors))
        density!(ax, rand(prior, samples), label=replace(string(name), "_" => " "))
    end

    axislegend(ax)

    return fig, ax
end

fig, ax = visualize_priors(Dict("Normal" => normal_prior,
                                "Log normal" => lognormal_prior,
                                "Scaled logit-normal" => logitnormal_prior))

xlims!(ax, -0.5, 1.5)
save("prior_flavors.svg", fig); nothing # hide

# ![](prior_flavors.svg)

# We note three important features:
#
# 1. `Normal` samples can be negative.
# 2. `Lognormal` samples cannot be negative, but can have large positive values.
# 3. `ScaledLogitNormal` samples have compact support, and thus strict upper and lower bounds.
#
# Boundedness is a very useful property of `ScaledLogitNormal`, so we
# explore specifying `ScaledLogitNormal` priors in more detail.
#
# # Fine-tuning a `ScaledLogitNormal` prior
#
# We can tune the "width" of a `ScaledLogitNormal` by supplying the
# standard-deviation `σ`:

logit_priors = Dict("Narrow (σ = 0.1)"                 => ScaledLogitNormal(σ = 0.1),
                    "Default (σ = 1)"                  => ScaledLogitNormal(σ = 1),
                    "Weird (σ = 4 with bounds=(0, 1))" => ScaledLogitNormal(σ = 4))

fig, ax = visualize_priors(logit_priors)
xlims!(ax, -0.1, 1.1)
save("logit_normal_widths.svg", fig); nothing #hide

# ![](logit_normal_widths.svg)

# Smaller standard deviations produce narrower prior distributions.
# However, weird shapes can result when a standard deviation is specified
# that's larger than the bounds interval (eg `σ > diff(bounds)`).
# In the above example we're using the default `bounds = (0, 1)`.
#
# Another way to build `ScaledLogitNormal` prior is to specify
# a central `interval` and `mass`. This allows us to shift the center of
# mass relative to the bounds (here we use the default `bounds=(0, 1)`),
# and provides a slightly safer way to tune a distribution.

shifted_priors = Dict("Shifted left"  => ScaledLogitNormal(interval=(0.2, 0.3), mass=0.9),
                      "Shifted right" => ScaledLogitNormal(interval=(0.7, 0.8), mass=0.9),
                      "Centered"      => ScaledLogitNormal(interval=(0.2, 0.8), mass=0.6))

fig, ax = visualize_priors(shifted_priors)
xlims!(ax, -0.1, 1.1)
save("logit_normal_shifting.svg", fig); nothing #hide

# ![](logit_normal_shifting.svg)
#
# The "Centered" case uses a `mass` that's proportional to the relative width spanned
# by `interval` (eg `mass = diff(interval) / diff(bounds)`) to produce a nearly
# uniform distribution within the bounds.


text:
@recipe function f(p::Passive{T}) where {T}
    seriestype --> :csg
    linecolor --> :silver
    seriescolor --> :silver
    fillalpha --> 0.2
    l = p.name != "" ? p.name : "Passive $(p.id)"
    @series begin
        #add empty line so that label is shown with alpha = 1
        seriestype := :path
        label --> l
        linewidth := 1
        T[], T[]
    end
    label := ""
    p.geometry
end

@recipe function f(sc::Semiconductor{T}) where {T}
    seriestype --> :csg
    linecolor --> :grey
    seriescolor --> :grey
    @series begin
        #add empty line so that label is shown with alpha = 1
        seriestype := :path
        label --> "Semiconductor"
        linewidth := 1
        T[], T[]
    end
    label := ""
    sc.geometry
end

@recipe function f(contact::Contact{T}) where {T}
    seriestype --> :csg
    seriescolor --> contact.id
    linecolor --> contact.id
    fillalpha --> 0.2
    l = contact.name != "" ? "$(contact.name) (id: $(contact.id))" : "Contact - id: $(contact.id)"
    @series begin
        #add empty line so that label is shown with alpha = 1
        seriestype := :path
        label --> l
        linewidth := 1
        T[], T[]
    end
    label := ""
    contact.geometry
end

@recipe function f(det::SolidStateDetector; show_semiconductor = false, show_passives = true, n_samples = 40)

    show_normal --> false

    plot_objects = []
    append!(plot_objects, det.contacts)
    if show_semiconductor
        push!(plot_objects, det.semiconductor)
    end
    if show_passives && !ismissing(det.passives)
        append!(plot_objects, det.passives)
    end

    if haskey(plotattributes, :seriestype) && plotattributes[:seriestype] == :samplesurface
        scales = [get_scale(o.geometry) for o in plot_objects]
        max_scale = maximum(scales)
    end

    for (i,o) in enumerate(plot_objects)
        @series begin
            if haskey(plotattributes, :seriestype) && plotattributes[:seriestype] == :samplesurface
                n_samples --> Int(ceil(n_samples * scales[i]/max_scale))
                CSG_scale := scales[i]
            end
            o
        end
    end
end


text:
module Mex

using MATLAB, Libdl

export jl_mex, input, call_matlab, is_interrupt_pending, check_for_interrupt

# --- Initialization --- #

# MATLAB mex interface
const libmex = Ref{Ptr{Cvoid}}()
const mex_call_matlab_with_trap = Ref{Ptr{Cvoid}}()

# MATLAB interrupt handling (from undocumented MATLAB)
const libut = Ref{Ptr{Cvoid}}()
const ut_is_interrupt_pending = Ref{Ptr{Cvoid}}()

function __init__()
    libmex[] = Libdl.dlopen(joinpath(MATLAB.matlab_libpath(), "libmex"), Libdl.RTLD_GLOBAL)
    mex_call_matlab_with_trap[] = Libdl.dlsym(libmex[], :mexCallMATLABWithTrap)

    libut[]  = Libdl.dlopen(joinpath(MATLAB.matlab_libpath(), "libut"), Libdl.RTLD_GLOBAL)
    ut_is_interrupt_pending[] = Libdl.dlsym(libut[], :utIsInterruptPending)
end

include("exceptions.jl")
include("juliatomatlab.jl")
include("matlabtojulia.jl")

end # module


text:
module RGBDSegmentation
using ColorTypes
using Statistics
using LinearAlgebra
using StaticArrays
using JLD
export clusterize
abstract type RGBDSegmentationAlgorithm end
import ImageSegmentation:
   SegmentedImage

include("Datasets/RGBDDataset.jl")
include("Processing.jl")

struct ColorDistanceNormalElement{T <: AbstractFloat }<: Colorant{T, 9}
    color::SVector{3, T}
    distance::SVector{3,T}
    normal::SVector{3, T}
end
struct ColorDistanceNormalSerializer{T <: AbstractFloat}
    data::Vector{T}
end

import Base.zero
import Base.show
import Base.+
import Base.-
function zero(c::Type{ColorDistanceNormalElement})
    return ColorDistanceNormalElement(
        zeros(SVector{3}),
        zeros(SVector{3}),
        zeros(SVector{3}))
end
function show(io::IO, elem::ColorDistanceNormalElement)
    show(io, elem.color)
    show(io, elem.distance)
    show(io, elem.normal)
end
function -(v1::ColorDistanceNormalElement, v2::ColorDistanceNormalElement)
    n = v1.normal - v2.normal
    n = n / norm(n)
    return ColorDistanceNormalElement(
        v1.color - v2.color,
        v1.distance - v2.distance,
        n)
end
function +(v1::ColorDistanceNormalElement, v2::ColorDistanceNormalElement)
    n = v1.normal + v2.normal
    n = n / norm(n)
    return ColorDistanceNormalElement(
        v1.color + v2.color,
        v1.distance + v2.distance,
        n)
end
function /(v1::ColorDistanceNormalElement, v::Integer)
    n = v1.normal / v
    n = n / norm(n)
    return ColorDistanceNormalElement(
        v1.color / v,
        v1.distance / v,
        n)
end
import Images.channelview
import Images.colorview
function channelview(img::Array{ColorDistanceNormalElement{T}, 2}) where T
    return reshape(reinterpret(T,  img), (9, size(img)...))
end
function colorview(::Type{RGB}, img::Array{ColorDistanceNormalElement{T}, 2}) where T
    return colorview(RGB, @view channelview(img)[1:3, :, :])
end
function colors(img::Array{ColorDistanceNormalElement, 2})
    colors = zeros(3, size(img)...)
    for I in CartesianIndices(size(img))
        colors[:,I] .= img[I].color
    end
    return colors
end

function JLD.readas(serdata::ColorDistanceNormalSerializer)
    return ColorDistanceNormalElement(
        SVector{3}(serdata.data[1:3]),
        SVector{3}(serdata.data[4:6]),
        SVector{3}(serdata.data[7:9]))
end
function JLD.writeas(data::ColorDistanceNormalElement)
    return ColorDistanceNormalSerializer(
        Vector(vcat(data.color, data.distance, data.normal))
    )
end

function SegmentedImage(img::CDNImage, labels::Matrix)
    label_list = convert.(Integer, sort(unique(labels)))
    region_means = Dict{Int, ColorDistanceNormalElement}()
    region_pix_count = Dict{Int, Int}()
    for label in label_list
        segment = findall(labels .== label)
        region_pix_count[label] = length(segment)
        mean_color = vec(mean(img[1:3, segment], dims=2))
        mean_dist = vec(mean(img[4:6, segment], dims=2))
        mean_normal = sum(img[7:9, segment], dims=2)
        mean_normal = vec(mean_normal / norm(mean_normal))
        region_means[label] = ColorDistanceNormalElement(
            SVector{3,Float64}(mean_color...),
            SVector{3,Float64}(mean_dist...),
            SVector{3,Float64}(mean_normal...))
    end
    return SegmentedImage(labels,
                          label_list,
                          region_means,
                          region_pix_count)
end

# Datasets
include("Datasets/NYUDepthV2.jl")

function clusterize(cfg::T, img::CDNImage, k::Integer) where T
    throw("Not implemented")
end

# Algorithms
include("Algorithms/DepthWMM/DepthWMM.jl")
include("Algorithms/JCSA/JCSA.jl")
include("Algorithms/CDNGraph/CDNGraph.jl")
include("Algorithms/GCF/GCF.jl")
end


text:
using LinearAlgebra
using CUDA

"""
	sim_activity_gpu(J::Matrix{Float32}, g::Float32, h0::Vector{Float32};
		s=0f0, b=nothing, dt=0.1f0, t_bound=1000.0f0, savestep=1,
		minimum_norm=0)

Integrates a tanh rate network with coupling matrix `J` and strength scaling `g`
from initial condition `h0`, returning a `N x NT` matrix of network unit activity.

When `minimum_norm > 0` the integration will be cut off when the network state
falls too close to the origin to prevent floating point underflow. The keyword
arguments `s` and `b` add self-coupling to `J` and add a bias term, respectively.
"""
function sim_activity_gpu(J::Matrix{Float32}, g::Float32, h0::Vector{Float32};
		s=0f0, b=nothing, dt=0.1f0, t_bound=1000.0f0, savestep=1,
		minimum_norm=0)

	N = size(J)[1]
	# Check network size and type consistency
	@assert N == size(J)[2] == length(h0)
	@assert Float32 == typeof(dt) == typeof(t_bound) == typeof(s)

	if b isa Array{Float32} && (size(b) == (N,) || size(b) == (N, 1))
		b = CuArray{Float32}(b)
		use_bias = true
	else
		use_bias = false
	end

	# Discretized simulation length
	maxstep = Int64(t_bound/dt)

	# Allocate memory for output
	out = CuArray{Float32,2}(undef, N, maxstep ÷ savestep)

	h = CuArray{Float32,1}(h0) # Set initial condition
	r = CuArray{Float32,1}(undef, N)

	J = CuArray{Float32,2}(g .* J + s * I)

	step = 0
	while step < maxstep
		step += 1
		map!(tanh, r, h)
		mul!(h, J, r, dt, 1 - dt)
		if use_bias
			h .+= dt .* b
		end
		if step % savestep == 0 # Save every specified time interval
			savept = step ÷ savestep
			out[:,savept] .= h
		end
		if norm(h) < minimum_norm
			D = -I + J .* ((sech.(h)).^2)'
			if all(real(eigvals(D)) .< 0)
				break
			end
		end
	end

	lastindex = step ÷ savestep
	return Matrix{Float32}(view(out, :, 1:lastindex))
end

"""
	sim_lyapunov_gpu(J::Matrix{Float32}, g::Float32, h0::Vector{Float32},
		dh::Vector{Float32}, Q0::Matrix{Float32}; s=0f0, b=nothing, dt=0.1f0,
		t_bound=1000.0f0, t_ONS=1.0f0, wait=0.0f0, return_max_LE=false,
		printstuff=false, scale_limit=100, p_factor=1f-1)

Integrates a tanh network with a given weight matrix and gain parameter from a
specified initial condition, and returns a Lyapunov spectrum calculated along
the trajectory.

Optionally returns the maximum Lyapunov exponent calculated using a direct
numerical simulation of a perturbed trajectory.
"""
function sim_lyapunov_gpu(J::Matrix{Float32}, g::Float32, h0::Vector{Float32},
		dh::Vector{Float32}, Q0::Matrix{Float32}; s=0f0, b=nothing, dt=0.1f0,
		t_bound=1000.0f0, t_ONS=1.0f0, wait=0.0f0, return_max_LE=false,
		printstuff=false, scale_limit=100, p_factor=1f-1)

	N = size(J, 1)
	@assert(size(J) == (length(h0), length(dh)) == size(Q0) == (N,N),
			"Dimensions of the coupling matrix, IC, initial ONS, and " *
			"perturbation must be consistent"
	)
	@assert(Float32 == typeof(s) == typeof(dt) == typeof(t_bound) == typeof(t_ONS)
			== typeof(wait), "Make sure all scalar parameters are the same type"
	)

	maxstep = round(Int64, t_bound / dt)
	startstep = round(Int64, wait / dt)
	tONSstep = round(t_ONS / dt)
	@assert(maxstep % tONSstep == startstep % tONSstep == 0,
		"t_ONS must exactly divide t_bound and wait"
	)

	if b isa Array{Float32} && (size(b) == (N,) || size(b) == (N, 1))
		b = CuArray{Float32,1}(reshape(b, (N,)))
		use_bias = true
	else
		use_bias = false
	end

	# Allocate memory for output
	spectrum = CUDA.zeros(N)

	# Allocate vector for activity and firing rate of current state
	h = CuArray{Float32,1}(h0)
	r = CuArray{Float32,1}(undef, N)

	# Load network structure and ICs to GPU
	J = CuArray{Float32,2}(g .* J + s * I)
	Q = CuArray{Float32,2}(Q0)

	# Allocate memory for Jacobian calculation and orthonormal basis
	sech2(x::Float32)::Float32 = 1 / (CUDA.cosh(x))^2
	sech2h = CuArray{Float32,2}(undef, 1, N)
	D = CuArray{Float32,2}(undef, N, N)
	dQ = CuArray{Float32,2}(undef, N, N)

	if return_max_LE
		max_LE = 0
		p_factor = convert(Float32, p_factor)
		dh = CuArray{Float32,1}(dh ./ norm(dh))
		h2 = h .+ p_factor .* dh
		r2 = CuArray{Float32,1}(undef, N)
	end

	step = 0
	while step < maxstep
		step += 1

		#Evolve the trajectory and the orthonormal basis in the tangent space
		map!(tanh, r, h)
		mul!(h, J, r, dt, 1 - dt)
		if use_bias
			h .+= dt .* b
		end

		map!(sech2, sech2h, h)
		copyto!(D, J)
		D .*= sech2h
		D[diagind(D)] .-= 1

		mul!(dQ, D, Q, dt, 0.0f0)
        Q .+= dQ

		#Evolve the perturbed trajectory and calculate how much it has grown
		if return_max_LE
			map!(tanh, r2, h2)
			mul!(h2, J, r2, dt, 1 - dt)
			if use_bias
				h2 .+= dt .* b
			end

			dh .= (h2 .- h) ./ p_factor
			scale = norm(dh)

			# If the perturbation scale has grown or shrunk too much, then
			# reset it to 1 and track the cumulative sum of log scales
			if scale > scale_limit || scale < 1 / scale_limit || step == startstep
				h2 .= h .+ dh ./ scale .* p_factor
				if step > startstep
					max_LE += log(scale)
				end
			end
		end

		# At each re-orthonormalization step, the log diagonal entries of R are
		# saved in a cumulative sum
		if step % tONSstep == 0
			printstuff && println("Time step: ", step)
			copyto!(dQ, Q)
			F = qr!(dQ)
			diags = view(F.factors, diagind(F.factors))
			if step > startstep  # forcing positive R-diagonal elements
				spectrum .+= log.(abs.(diags))
			end
			fill!(Q, 0)
			map!(sign, view(Q, diagind(Q)), diags)
            lmul!(F.Q, Q)
		end
	end

	# Update the maximum Lyapunov exponent after the simulation has ended so
	# that the full simulation time is used for the calculation
	true_time = dt * (maxstep - startstep)
	if return_max_LE
		max_LE += log(norm((h2 .- h) ./ p_factor))
		return Vector{Float32}(spectrum ./ true_time), max_LE / true_time
	else
		return Vector{Float32}(spectrum ./ true_time)
	end
end


text:
using Test
using LogicCircuits

include("helper/little_circuits.jl")

@testset "Plain logic nodes" begin
    n1 = little_2var()
    n0 = little_3var()
    n0c = little_3var_constants()
    r1 = fully_factorized_circuit(PlainLogicCircuit,10)

    @test length(and_nodes(n0c)) == 6
    @test length(or_nodes(n0c)) == 5

    @test variable(left_most_descendent(n1)) == Var(1)
    @test ispositive(left_most_descendent(n1))
    @test !isnegative(left_most_descendent(n1))

    @test istrue(compile(PlainLogicCircuit,true))
    @test istrue(PlainLogicCircuit(true))
    @test istrue(LogicCircuit(true))
    @test isfalse(compile(PlainLogicCircuit,false))
    @test isfalse(PlainLogicCircuit(false))
    @test isfalse(LogicCircuit(false))
    @test !istrue(compile(PlainLogicCircuit,false))
    @test !isfalse(compile(PlainLogicCircuit,true))
    @test !istrue(compile(PlainLogicCircuit,Lit(2)))
    @test !isfalse(compile(PlainLogicCircuit,Lit(2)))
    @test !istrue(n1)
    @test !isfalse(n1)

    @test literal(compile(PlainLogicCircuit,Lit(2))) == Lit(2)
    @test literal(PlainLogicCircuit(Lit(2))) == Lit(2)
    @test literal(LogicCircuit(Lit(2))) == Lit(2)

    @test_throws Exception conjoin(PlainLogicCircuit[])
    @test_throws Exception disjoin(PlainLogicCircuit[])
    a1 = conjoin(n1, n0, n0c)
    o1 = disjoin(n1, n0, n0c)

    @test conjoin([n1, n0, n0c]; reuse=a1) == a1
    @test disjoin([n1, n0, n0c]; reuse=o1) == o1
    @test disjoin([n1, n0, n0c]; reuse=a1) != a1
    @test conjoin([n1, n0, n0c]; reuse=o1) != o1

    @test num_nodes(fully_factorized_circuit(PlainLogicCircuit,10)) == 32
    @test num_edges(fully_factorized_circuit(PlainLogicCircuit,10)) == 31

    @test tree_formula_string(n1) == "((1 ⋀ (2 ⋁ -2)) ⋁ (-1 ⋀ (2 ⋁ -2)))"
    @test tree_formula_string(n0) == "((3 ⋀ ((1 ⋀ (2 ⋁ -2)) ⋁ (-1 ⋀ (2 ⋁ -2)))) ⋁ (-3 ⋀ ((1 ⋀ (2 ⋁ -2)) ⋁ (-1 ⋀ (2 ⋁ -2)))))"
    @test tree_formula_string(n0c) == "(((3 ⋀ true) ⋀ ((1 ⋀ (2 ⋁ -2)) ⋁ (-1 ⋀ (2 ⋁ -2)))) ⋁ ((-3 ⋀ false) ⋀ ((1 ⋀ (2 ⋁ -2)) ⋁ (-1 ⋀ (2 ⋁ -2)))))"

    con_map = canonical_constants(n0c)
    @test isfalse(con_map[1])
    @test istrue(con_map[2])
    @test length(con_map) == 2

    lit_map = canonical_literals(r1)
    @test literal(lit_map[Lit(1)]) == Lit(1)
    @test literal(lit_map[Lit(-5)]) == Lit(-5)
    @test length(lit_map) == 20

    @test canonical_constants(r1) == (nothing, nothing)

    @test node_stats(n0c)[PlainConstantNode] == 2

    io = IOBuffer()
    show(io,n0c)
    @test length(String(take!(io))) > 0

    n0c2 = LogicCircuit(n0c)
    @test n0c2 isa PlainLogicCircuit
    @test n0c2 !== n0c
    @test num_edges(n0c2) == num_edges(n0c)
    @test num_nodes(n0c2) == num_nodes(n0c)
    @test all(isleaf, intersect(linearize(n0c2),linearize(n0c)))

    n0c2 = PlainLogicCircuit(n0c)
    @test n0c2 isa PlainLogicCircuit
    @test n0c2 !== n0c
    @test num_edges(n0c2) == num_edges(n0c)
    @test num_nodes(n0c2) == num_nodes(n0c)
    @test all(isleaf, intersect(linearize(n0c2),linearize(n0c)))

end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
using Documenter
using DocStringExtensions
using Hypatia

makedocs(
    sitename = "Hypatia",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [Hypatia, Hypatia.Cones, Hypatia.Models, Hypatia.Solvers, Hypatia.PolyUtils],
    pages = [
        "Home" => "index.md",
        "Cones" => "cones.md",
        "Modeling" => "modeling.md",
        "Solving" => "solving.md",
        "Examples" => "examples.md",
        "Benchmarks" => "benchmarks.md",
        "API reference" => [
            "api/hypatia.md",
            "api/cones.md",
            "api/models.md",
            "api/solvers.md",
            "api/polyutils.md",
            ],
        ],
    # linkcheck = true, # test URLs
    )

deploydocs(
    repo = "github.com/chriscoey/Hypatia.jl.git",
    target = "build",
    push_preview = true,
    )


text:
println("started")

# ]add BenchmarkTools
# ]add MAT
using BenchmarkTools # for benchmark
const BenchmarkTools.DEFAULT_PARAMETERS.samples = 700;
const BenchmarkTools.DEFAULT_PARAMETERS.evals = 1;
using DelimitedFiles # for readdlm, writedlm
using Statistics # for median
using LinearAlgebra # for eigen, pinv, svd, cholesky
using Random # for randperm
using MAT

const operationMode = 2; # 0 for test only # 1 for partial benchmark # 2 for full benchmark
const saveData = true; # true for saving Runtimes to csv and mat
const JuliaFlag = true; # run Julia benchmark
const JuliaSIMDFlag = true; # run Julia SIMD benchmark
const debugging = false;

if BLAS.vendor() == :mkl
    # MKL

    ### Julia
    if JuliaFlag
        include("JuliaBench.jl")
        tRunTime, mRunTime = JuliaBench(operationMode)

        if saveData
            # Main RunTime Table write
            writedlm(
                joinpath(
                    "RunTimeData",
                    "RunTimeJulia$(BLAS.vendor())Table.csv",
                ),
                tRunTime,
                ',',
            )
            # RunTime save
            file = matopen(
                joinpath("RunTimeData", "RunTimeJulia$(BLAS.vendor()).mat"),
                "w",
            )
            write(file, "mRunTime", mRunTime)
            close(file)
        end
    end

    ### Julia SIMD
    if JuliaSIMDFlag
        include("JuliaBenchSIMD.jl")
        tRunTime, mRunTime = JuliaBenchSIMD(operationMode)

        if saveData
            # Main RunTime Table write
            writedlm(
                joinpath(
                    "RunTimeData",
                    "RunTimeJulia$(BLAS.vendor())SIMDTable.csv",
                ),
                tRunTime,
                ',',
            )
            # RunTime save
            file = matopen(
                joinpath("RunTimeData", "RunTimeJulia$(BLAS.vendor())SIMD.mat"),
                "w",
            )
            write(file, "mRunTime", mRunTime)
            close(file)
        end
    end

else
    # openBLAS

    ## BLAS threads = 1
    BLAS.set_num_threads(1)

    ### Julia
    if JuliaFlag
        include("JuliaBench.jl")
        tRunTime, mRunTime = JuliaBench(operationMode)

        if saveData
        # Main RunTime Table write
            writedlm(
                joinpath(
                    "RunTimeData",
                    "RunTimeJulia$(BLAS.vendor())Table.csv",
                ),
                tRunTime,
                ',',
            )
        # RunTime save
            file = matopen(
                joinpath("RunTimeData", "RunTimeJulia$(BLAS.vendor()).mat"),
                "w",
            )
            write(file, "mRunTime", mRunTime)
            close(file)
        end
    end

    ### Julia SIMD
    if JuliaSIMDFlag
        include("JuliaBenchSIMD.jl")
        tRunTime, mRunTime = JuliaBenchSIMD(operationMode)

        if saveData
        # Main RunTime Table write
            writedlm(
                joinpath(
                    "RunTimeData",
                    "RunTimeJulia$(BLAS.vendor())SIMDTable.csv",
                ),
                tRunTime,
                ',',
            )
        # RunTime save
            file = matopen(
                joinpath("RunTimeData", "RunTimeJulia$(BLAS.vendor())SIMD.mat"),
                "w",
            )
            write(file, "mRunTime", mRunTime)
            close(file)
        end
    end

    ## BLAS threads = 4
    BLAS.set_num_threads(4)

    ### Julia
    if JuliaFlag
        include("JuliaBench.jl")
        tRunTime, mRunTime = JuliaBench(operationMode)

        if saveData
        # Main RunTime Table write
            writedlm(
                joinpath(
                    "RunTimeData",
                    "RunTimeJulia$(BLAS.vendor())Table_4Thread.csv",
                ),
                tRunTime,
                ',',
            )
        # RunTime save
            file = matopen(
                joinpath(
                    "RunTimeData",
                    "RunTimeJulia$(BLAS.vendor())_4Thread.mat",
                ),
                "w",
            )
            write(file, "mRunTime", mRunTime)
            close(file)
        end
    end

    ### Julia SIMD
    if JuliaSIMDFlag
        include("JuliaBenchSIMD.jl")
        tRunTime, mRunTime = JuliaBenchSIMD(operationMode)

        if saveData
        # Main RunTime Table write
            writedlm(
                joinpath(
                    "RunTimeData",
                    "RunTimeJulia$(BLAS.vendor())SIMDTable_4Thread.csv",
                ),
                tRunTime,
                ',',
            )
        # RunTime save
            file = matopen(
                joinpath(
                    "RunTimeData",
                    "RunTimeJulia$(BLAS.vendor())SIMD_4Thread.mat",
                ),
                "w",
            )
            write(file, "mRunTime", mRunTime)
            close(file)
        end
    end

end

if debugging
    # Debug and performance trace
    include("JuliaBench.jl")
    matrixSize = 20
    mX = randn(matrixSize, matrixSize)
    mY = randn(matrixSize, matrixSize)
    allFunctions = [
        MatrixGeneration,
        MatrixAddition,
        MatrixMultiplication,
        MatrixQuadraticForm,
        MatrixReductions,
        ElementWiseOperations,
        MatrixExp,
        MatrixSqrt,
        Svd,
        Eig,
        CholDec,
        MatInv,
        LinearSystem,
        LeastSquares,
        CalcDistanceMatrix,
        KMeans,
    ]
    # debug the code
    # Juno.@enter allFunctions[2](matrixSize, mX, mY) # copy paste to REPL
    # timing and memory
    @time allFunctions[2](matrixSize, mX, mY)
    # bad practices used in coding
    # Juno.@trace allFunctions[1](matrixSize, mX, mY)
    # profiles the code
    # using Profile
    # Profile.clear();
    # @profile allFunctions[2](matrixSize, mX, mY)
    # Profile.print(combine = true)
    # profiletree=Juno.profiletree()
    # Juno.profiler()
    # advanced tool for diagnosing type-related problems
    # @code_warntype allFunctions[2](matrixSize, mX, mY)
end

max_stars_repo_name: aminya/Julia-Matlab-Benchmark

max_issues_repo_name: aminya/Julia-Matlab-Benchmark

max_forks_repo_name: aminya/Julia-Matlab-Benchmark


text:
abstract type _HTTP_METHOD end
abstract type _HTTP_METHOD_GET <: _HTTP_METHOD end
abstract type _HTTP_METHOD_DELETE <: _HTTP_METHOD end
abstract type _HTTP_METHOD_POST <: _HTTP_METHOD end
print(io::IO, ::Type{_HTTP_METHOD_GET}) = print(io, "GET")
print(io::IO, ::Type{_HTTP_METHOD_DELETE}) = print(io, "DELETE")
print(io::IO, ::Type{_HTTP_METHOD_POST}) = print(io, "POST")

function _sign(signature::String, data::String)
    return base64encode(hmac_sha256(Vector{UInt8}(signature), data))
end

function _generate_headers(
    api_data::ApiData,
    http_method_type::Type{<:_HTTP_METHOD},
    endpoint_path::String,
    request_data::String,
)
    ts_str = string(round(Int64, time() * 1000))
    str_to_sign = "$ts_str$http_method_type$endpoint_path$request_data"
    return (
        "KC-API-SIGN" => _sign(api_data.secret, str_to_sign),
        "KC-API-TIMESTAMP" => ts_str,
        "KC-API-KEY" => api_data.key,
        "KC-API-PASSPHRASE" => api_data.encrypted_passphrase,
        "KC-API-KEY-VERSION" => 2,
        "Content-Type" => "application/json",
    )
end

function _create_api_uri(endpoint_path::String, query::String)
    return URI(; scheme="https", host="api.kucoin.com", path=endpoint_path, query=query)
end

function _create_api_uri(endpoint_path::String)
    return URI(; scheme="https", host="api.kucoin.com", path=endpoint_path)
end

function _handle(response::Response)::Union{JSON3.Array{JSON3.Object},JSON3.Object}
    iserror(response) && error(response.body)
    data = JSON3.read(response.body)
    if haskey(data, "code")
        data["code"] == "200000" || error(JSON3.write(data))
        return data["data"]
    end
    return data
end

function _kucoin_request(
    http_method_type::Type{<:_HTTP_METHOD}, endpoint_path::String; kw...
)
    query_str = join("$argument=$value" for (argument, value) ∈ kw if value ≢ nothing)
    uri = _create_api_uri(endpoint_path, query_str)
    return request(http_method_type, uri; query=query_str)
end

function _kucoin_request(
    api_data::ApiData,
    http_method_type::Union{Type{_HTTP_METHOD_GET},Type{_HTTP_METHOD_DELETE}},
    endpoint_path::String;
    kw...,
)
    query_str = join(
        ("$argument=$value" for (argument, value) ∈ kw if value ≢ nothing), "&"
    )
    uri = _create_api_uri(endpoint_path, query_str)
    request_data = isempty(query_str) ? "" : "?$query_str"
    headers = _generate_headers(api_data, http_method_type, endpoint_path, request_data)
    return request(http_method_type, uri; query=query_str, headers=headers)
end

function _kucoin_request(
    api_data::ApiData,
    http_method_type::Type{_HTTP_METHOD_POST},
    endpoint_path::String;
    kw...,
)
    uri = _create_api_uri(endpoint_path)
    json_str = JSON3.write(kw)
    headers = _generate_headers(api_data, http_method_type, endpoint_path, json_str)
    return request(http_method_type, uri; body=json_str, headers=headers)
end


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

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export DropoutLayer
# We use @eval begin ... end blocks here to delay the expansion of the @neuron
# macro which causes side effects, that is, it should only happen once for the
# active definition of dropout neuron
if LATTE_BATCH_DROPOUT
@eval begin
    @neuron type DropoutNeuron
        ratio   :: Float32
        randval :: Float32
        scale   :: Float32
    end
    DropoutNeuron(ratio::Float32) = DropoutNeuron(ratio, 0.0f0, 1.0f0 / (1.0f0 - ratio))
end
else
@eval begin
    @neuron type DropoutNeuron
        ratio   :: Float32
        randval :: Batch{Float32}
        scale   :: Float32
    end
    DropoutNeuron(ratio::Float32) = DropoutNeuron(ratio, Batch(0.0f0), 1.0f0 / (1.0f0 - ratio))
end
end

# FIXME: CGen does not support rand()
# ccall((:srand48, "libc"), Void, (Clong,), ccall((:time, "libc"), Clong, (Ptr{Void},), C_NULL))
# function float_rand()
#     val = ccall((:drand48, "libc"), Cdouble, ())
# end

# if LATTE_BATCH_DROPOUT
    @neuron forward(neuron::DropoutNeuron) do
        if neuron.randval > neuron.ratio
            neuron.value = neuron.inputs[1] * neuron.scale
        else
            neuron.value = 0.0
        end
    end
# else
#     @neuron forward(neuron::DropoutNeuron) do
#         neuron.randval = float_rand()
#         if neuron.randval > neuron.ratio
#             neuron.value = neuron.inputs[1] / (1.0f0 - neuron.ratio)
#         else
#             neuron.value = 0.0
#         end
#     end
# end

@neuron backward(neuron::DropoutNeuron) do
    if neuron.randval > neuron.ratio
        neuron.∇inputs[1] = neuron.∇ * neuron.scale
    else
        neuron.∇inputs[1] = 0.0
    end
end

function DropoutLayer(name::Symbol, net::Net, input_ensemble::AbstractEnsemble, ratio=0.5f0)
    neurons = Array(DropoutNeuron, size(input_ensemble)...)
    for i in 1:length(neurons)
        neurons[i] = DropoutNeuron(ratio)
    end
    ActivationEnsemble(net, name, neurons, input_ensemble; phase=Train)
end


text:
using MPSToolkit
using Test

@testset "TransverseIsing1D" begin
    model = TransverseIsing1D(3)
    @test length(getGatesEven(model, 0.1)) == 1
    @test length(getGatesOdd(model, 0.1)) == 2

    model = TransverseIsing1D(4)
    @test length(getGatesEven(model, 0.1)) == 2
    @test length(getGatesOdd(model, 0.1)) == 2

    # test ground state in ordered phase
    dmrg_options = DMRGOptions()
    for L in [4, 5, 8, 16]
        model = TransverseIsing1D(L)
        model.parameters["hx"] = 0.0
        model.parameters["hz"] = 1.0
        results = runDMRG(model, dmrg_options)
        @test results.observables[1]["H"][1] ≈ -2 * L + 1
        @test abs(results.observables[1]["Sx"][1]) < 1e-12
        @test results.observables[1]["Sz"][1] ≈ L
        @test results.localOperators[1]["sz"] ≈ fill(1.0, (L,))
        # FIXME: test sx and sz_sz
    end
end

@testset "Gates" begin
    model = TransverseIsing1D(2)
    model.parameters["J"] = 1.0
    model.parameters["hx"] = 0.0
    model.parameters["hz"] = 0.1
    even = getGatesEven(model, 0.1)
    odd = getGatesOdd(model, 0.3)
    @test length(even) == 1
    @test length(odd) == 1

    psi = productMPS(model.sites, n -> "Dn")

    @test MPSToolkit.ITensors.expect(psi, "Sz") ≈ [-0.5, -0.5]
    for i = 1:10
        psi = MPSToolkit.ITensors.apply(odd, psi)
    end
    @test MPSToolkit.ITensors.expect(psi, "Sz") ≈ [-0.5, -0.5]
    for i = 1:10
        psi = MPSToolkit.ITensors.apply(even, psi)
    end
    @test MPSToolkit.ITensors.expect(psi, "Sz") ≈ [-0.5, -0.5]
end


text:
Triangle -> Triangle
Incremental -> Incremental
Accumulated -> Accumulated
IncrementalCompleted -> IncrementalCompleted
AccumulatedCompleted -> AccumulatedCompleted
taccum -> taccum
tdeacum -> tdeacum
tfactors -> tfactors
fillcl -> fillcl
ibnrcl -> ibnrcl
tquotas -> tquotas
expositionbf -> expositionbf
fillbf -> fillbf
ibnrbf -> ibnrbf
randtri
lossesbf -> lossesbf
incurred -> incurred
lossratiosbf -> lossratiosbf


text:
module Soundcard
using SharedArrays
using Libaudio

# int record(float * pcm_record, int64_t record_channels, int64_t record_frames, int64_t samplerate);
# int play(const float * pcm_play, int64_t play_channels, int64_t play_frames, int64_t samplerate);
# int playrecord(const float * pcm_play, int64_t play_channels, float * pcm_record, int64_t record_channels, int64_t common_frames, int64_t samplerate);

modulepath(name) = realpath(joinpath(dirname(pathof(name)),".."))
folder() = "C:/DRIVERS/Julia/"

"""
    init()

install binary dependencies to "C:\\Drivers\\Julia\\"
"""
function __init__()
    root = folder()
    isdir(root) || mkpath(root)
    mp = modulepath(Soundcard)

    t1 = joinpath(root, "portaudio_x64.dll")
    t2 = joinpath(root, "soundcard_api.dll")

    isfile(t1) || cp(joinpath(mp, "deps/usr/lib/portaudio_x64.dll"), t1, force=true)
    isfile(t2) || cp(joinpath(mp, "deps/usr/lib/soundcard_api.dll"), t2, force=true)
end

"""
    device()

populate soundcard devices attached to the workstation
"""
function device()
    buffer = zeros(Int8, 8192)
    numdev = ccall((:list_devices, "C:\\Drivers\\Julia\\soundcard_api"), Int32, (Ptr{Int8},), buffer)
    digest = ""
    for i in buffer
        digest = digest * string(Char(i))
    end
    report = split(digest,'\n')
    numdev, report[1:end-1]
end

"""
    record(dim, fs)

low level api for recording
"""
function record(dim::Tuple{Int64, Int64}, fs::Int64)    # -> Matrix{Float32}
    pcm = zeros(Float32, dim[2] * dim[1])
    ccall((:record, "C:\\Drivers\\Julia\\soundcard_api"), Int32, (Ptr{Float32}, Int64, Int64, Int64), pcm, dim[2], dim[1], fs)
    return transpose(reshape(pcm, dim[2], dim[1]))
end

"""
    play(dat, fs)

low level api for playback
"""
function play(dat::Matrix{Float32}, fs::Int64)
    pcm = interleave(dat)
    ccall((:play, "C:\\Drivers\\Julia\\soundcard_api"), Int32, (Ptr{Float32}, Int64, Int64, Int64), pcm, size(dat)[2], size(dat)[1], fs)
    return nothing
end

"""
    playrecord(dat, ch, fs)

low level api for simultaneous playback and recording
usually we have ch = size(mixmic,1)
"""
function playrecord(dat::Matrix{Float32}, ch::Int64, fs::Int64)    # -> Matrix{Float32}
    pcmo = interleave(dat)
    pcmi = zeros(Float32, size(dat)[1] * ch)
    ccall((:playrecord, "C:\\Drivers\\Julia\\soundcard_api"), Int32, (Ptr{Float32}, Int64, Ptr{Float32}, Int64, Int64, Int64), pcmo, size(dat)[2], pcmi, ch, size(dat)[1], fs)
    return transpose(reshape(pcmi, ch, size(dat)[1]))
end

function interleave(x::Matrix{T}) where T <: Number
    fr,ch = size(x)
    y = zeros(T, ch * fr)
    k::Int64 = 0
    for i = 1:fr
        y[k+1:k+ch] = x[i,:]
        k += ch
    end
    return y
end

"""
    mixer(x, mix)

apply mixing matrix 'mix' to input signal matrix 'x' for logical signal routing,
mixing matrix will be converted to the same type as the input matrix, note that
this function is the last gate to the physical world of soundcard so use clamp
"""
function mixer(x::AbstractMatrix{T}, mix::AbstractMatrix) where T <: AbstractFloat
    mm = convert(Matrix{T}, mix)
    y = x * mm
    maximum(abs.(y)) >= one(T) && Libaudio.printl(joinpath(Libaudio.folder(), Libaudio.logfile()), :light_red, Libaudio.nows() * " | Soundcard.mixer: sample amplitude clipping")
    return clamp.(y, -one(T), one(T))
end

                                    # --------------------------
                                    #    API for ease of use
                                    # --------------------------

function playrecord(playing::Matrix, ms::Matrix, mm::Matrix, fs)::Matrix{Float32}
    out = convert(Matrix{Float32}, playing)
    recording = mixer(playrecord(mixer(out,ms), size(mm,1), convert(Int64,fs)), mm)
end

function play(playing::Matrix, ms::Matrix, fs)
    out = convert(Matrix{Float32}, playing)
    play(mixer(out, ms), Int64(fs))
end

function record(samples::Integer, mm::Matrix, fs)::Matrix{Float32}
    recording = mixer(record((convert(Int64,samples),size(mm,1)), Int64(fs)), mm)
end

                                    # --------------------------
                                    #  API for ultralow latency
                                    # --------------------------
# x = Soundcard.mixer(convert(Matrix{Float32},out), ms)
# y = SharedArray{Float32,1}(Soundcard.interleave(x))
play(size_x::Tuple{Int64,Int64}, y::SharedArray{Float32,1}, fs) =
    ccall((:play,"C:\\Drivers\\Julia\\soundcard_api"), Int32, (Ptr{Float32}, Int64, Int64, Int64), y, size_x[2], size_x[1], fs)  # remotecall
# async tasks...
# fetch(done)

# y = SharedArray{Float32,1}(zeros(Float32, size(mm,1)*samples))
record(y::SharedArray{Float32,1}, size_mm::Tuple{Int64,Int64}, samples, fs) =
    ccall((:record,"C:\\Drivers\\Julia\\soundcard_api"), Int32, (Ptr{Float32}, Int64, Int64, Int64), y, size_mm[1], samples, fs)  # remotecall
# async tasks...
# fetch(done)
# recording = Soundcard.mixer(transpose(reshape(y,size(mm,1),samples)), mm)

# samples = size(out,1)
# x = Soundcard.mixer(convert(Matrix{Float32},out), ms)
# y = SharedArray{Float32,1}(Soundcard.interleave(x))
# z = SharedArray{Float32,1}(zeros(Float32,size(mm,1)*samples))
playrecord(size_x::Tuple{Int64,Int64}, y::SharedArray{Float32,1}, z::SharedArray{Float32,1}, size_mm::Tuple{Int64,Int64}, fs) =
    ccall((:playrecord, "C:\\Drivers\\Julia\\soundcard_api"), Int32, (Ptr{Float32}, Int64, Ptr{Float32}, Int64, Int64, Int64), y, size_x[2], z, size_mm[1], size_x[1], fs)
# async tasks...
# fetch(done)
# recording = Soundcard.mixer(transpose(reshape(z,size(mm,1),samples)), mm)

end # module


text:
using SimpleDirectMediaLayer.LibSDL2
using Test

SDL2_pkg_dir = joinpath(@__DIR__, "..", "..")

@testset "Init+Quit" begin
    # Test that you can init and quit SDL_ttf multiple times correctly.
    @test SDL_Init(SDL_INIT_EVERYTHING) == 0
    @test TTF_Init() == 0
    TTF_Quit()
    SDL_Quit()

    @test SDL_Init(SDL_INIT_EVERYTHING) == 0
    @test TTF_Init() == 0
    TTF_Quit()
    SDL_Quit()
end

@testset "Text" begin
    @test SDL_Init(SDL_INIT_EVERYTHING) == 0
    @test TTF_Init() == 0

    @testset "Simple text" begin
        font = TTF_OpenFont(joinpath(SDL2_pkg_dir, "assets/fonts/FiraCode/ttf/FiraCode-Regular.ttf"), 14)
        @test font != C_NULL
        txt = "Hello World!"
        text = TTF_RenderText_Blended(font, txt, SDL_Color(20, 20, 20, 255))
        @test text != C_NULL

        fx, fy = Ref{Cint}(0), Ref{Cint}(0)
        @test TTF_SizeText(font, txt, fx, fy) == 0
        @test fx[] > 0
        @test fy[] > 0
    end

    TTF_Quit()
    SDL_Quit()
end

@testset "Rendering" begin
    @test SDL_Init(SDL_INIT_EVERYTHING) == 0
    @test TTF_Init() == 0
    win = SDL_CreateWindow("Hello World!", 100, 100, 300, 400, 0)
    renderer = SDL_CreateRenderer(win, -1, 0)

    # try rendering with complicated text
    @testset "Special Characters text" begin
        font = TTF_OpenFont(joinpath(SDL2_pkg_dir, "assets/fonts/FiraCode/ttf/FiraCode-Regular.ttf"), 14)
        @test font != C_NULL
        txt = "@BinDeps.install Dict([ (:glib, :libglib) ])"
        text = TTF_RenderText_Blended(font, txt, SDL_Color(20, 20, 20, 255))
        @test text != C_NULL
        tex = SDL_CreateTextureFromSurface(renderer, text)
        @test tex != C_NULL

        fx, fy = Ref{Cint}(0), Ref{Cint}(0)
        @test TTF_SizeText(font, txt, fx, fy) == 0
        @test fx[] > 0
        @test fy[] > 0
        @test SDL_RenderCopy(renderer, tex, C_NULL, Ref(SDL_Rect(0, 0, fx[], fy[]))) == 0

        SDL_RenderPresent(renderer)
    end

    TTF_Quit()
    SDL_Quit()
end


text:
# including files in the path
include("04a-prismMod.jl")
include("04b-circleMod.jl")

# importing modules into code
using .Prism_Mod
import .Circle_Mod

function main()

    # accessing Prism module
    # note how these functions use its namespace to access its use
    # namespace as in Prism_Mod.function()
    p = Prism_Mod.Prism()
    p1 = Prism_Mod.Prism(1, 2, 3)

    # Prism module exported function
    # while this function doesn't use the namespace
    println(volume(p))

    # accessing Circle module
    # note how even though this does have export in its module
    # it still needs its namespace i.e. Circle_Mod.function()
    c = Circle_Mod.Circle_const(5)
    println(c.radius)

end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
# This part is used for loading the data for testing purpose
const test_file_url_M6707HH = "https://www.cfa.harvard.edu/~dfabricant/huchra/ay16/fits/M6707HH.fits"
const test_file_path_M6707HH = joinpath(@__DIR__, "data", "M6707HH.fits")
if !isfile(test_file_path_M6707HH)
    mkpath(dirname(test_file_path_M6707HH))
    download(test_file_url_M6707HH, test_file_path_M6707HH)
end
const M6707HH =  FITS(test_file_path_M6707HH)

const test_file_url_M35070V = "https://www.cfa.harvard.edu/~dfabricant/huchra/ay16/fits/M35070V.fits"
const test_file_path_M35070V = joinpath(@__DIR__, "data", "M35070V.fits")
if !isfile(test_file_path_M35070V)
    mkpath(dirname(test_file_path_M35070V))
    download(test_file_url_M35070V, test_file_path_M35070V)
end
const M35070V =  FITS(test_file_path_M35070V)


text:
""" `zpk(gain[, Ts]), zpk(num, den, k[, Ts]), zpk(sys)`

Create transfer function on zero pole gain form. The numerator and denominator are represented by their poles and zeros.

`sys::TransferFunction{SisoZpk{T,TR}} = k*numerator/denominator`
where `T` is the type of `k` and `TR` the type of the zeros/poles, usually Float64 and Complex{Float64}.

`num`: the roots of the numerator polynomial. Either scalar or vector to create SISO systems
or an array of vectors to create MIMO system.

`den`: the roots of the denominator polynomial. Either vector to create SISO systems
or an array of vectors to create MIMO system.

`k`: The gain of the system. Obs, this is not the same as `dcgain`.

`Ts`: Sample time if discrete time system.

Other uses:

`zpk(sys)`: Convert `sys` to `zpk` form.

`zpk("s")`: Create the transferfunction `s`.

"""
function zpk(z::VecOrMat{<:AbstractVector{TZ}}, p::VecOrMat{<:AbstractVector{TP}}, k::VecOrMat{T0}, Ts::TE) where {TE<:TimeEvolution, T0<:Number, TZ<:Number, TP<:Number}
    # Validate input and output dimensions match
    if !(size(z) == size(p) == size(k))
        error("Dimensions for z, p, and k must match")
    end

    TR = promote_type(T0,TZ,TP)
    T = promote_type(T0, real(TR)) # Ensuring this avoids many problems, e.g., with SisoZpk{Int64,ComplexF64}

    matrix = Array{SisoZpk{T, TR}}(undef, size(z,1), size(z,2)) # TODO: make this nicer
    for o=1:size(z,1)
        for i=1:size(z,2)
            matrix[o, i] = SisoZpk{T,TR}(z[o, i], p[o, i], k[o, i])
        end
    end
    return TransferFunction{TE,SisoZpk{T,TR}}(matrix, Ts)
end
function zpk(z::AbstractVector{TZ}, p::AbstractVector{TP}, k::T, Ts::TE) where {TE<:TimeEvolution, T<:Number, TZ<:Number, TP<:Number}
    return zpk(fill(z, 1, 1), fill(p, 1, 1), fill(k, 1, 1), Ts)
end

function zpk(z::AbstractVector, p::AbstractVector, k::T, Ts::TE) where {TE<:TimeEvolution, T<:Number} # To be able to send in empty vectors [] of type Any
    if eltype(z) == Any
        z == [] || throw(ArgumentError("non numeric vectors must be empty."))
        z = T[]
    end
    if eltype(p) == Any
        p == [] || throw(ArgumentError("non numeric vectors must be empty."))
        p = T[]
    end
    zpk(z, p, k, Ts)
end

function zpk(gain::Matrix{T}, Ts::TE; kwargs...) where {TE<:TimeEvolution, T <: Number}
    TR = promote_type(ComplexF64,T)
    ny, nu = size(gain)
    matrix = [SisoZpk{T, TR}(TR[],TR[], gain[o, i]) for o=1:ny, i=1:nu]
    return TransferFunction{TE,SisoZpk{T,TR}}(matrix, Ts)
end

zpk(k::Real, Ts::TimeEvolution) = zpk(eltype(k)[], eltype(k)[], k, Ts)

zpk(sys::StateSpace) = zpk(zpkdata(sys)...)

function zpk(G::TransferFunction{TE,S}) where {TE<:TimeEvolution,T0, S<:SisoTf{T0}}
    T = typeof(one(T0)/one(T0))
    convert(TransferFunction{TE,SisoZpk{T, complex(T)}}, G)
end

zpk(var::AbstractString) = zpk(tf(var))
zpk(var::AbstractString, Ts::Real) = zpk(tf(var, Ts))

# Catch all 3(4) argument versions
zpk(z, p, k, Ts::Number) = zpk(z, p, k, Discrete(Ts))
zpk(z, p, k) = zpk(z, p, k, Continuous())
# Catch all 1(2) argument versions
zpk(gain, Ts::Number; kwargs...) where {T <: Number} = zpk(gain, Discrete(Ts))
zpk(gain; kwargs...) where {T <: Number} = zpk(gain, Continuous())

# This method is required since the Base method uses promote_op(matprod, eltype(A), eltype(B)) which fails to infer correctly
function Base.:(*)(A::AbstractMatrix{<:Union{<:SisoZpk, <:SisoRational}}, B::AbstractMatrix{<:Union{<:SisoZpk, <:SisoRational}})
    TS = promote_op(*, eltype(A), eltype(B))
    mul!(similar(B, TS, (size(A,1), size(B,2))), A, B)
end