text:
export DMABufferVector, PageAlignedArray, PageAlignedVector

"""
    struct DMABufferVector{S,T} <: AbstractVector{S}

AlazarTech digitizers use direct memory access (DMA) to transfer data from digitizers to
the computer's main memory. This struct abstracts memory buffers on the host. The elements
of a `DMABufferVector` are pointers to the individual buffers, which are each page-aligned,
provided a page-aligned backing array is used (e.g. `Base.SharedVector` or
`Alazar.PageAlignedVector`).

`DMABufferVector` may be constructed as, for example,
`DMABufferVector(SharedVector{UInt8}, bytes_buf, n_buf)` or
`DMABufferVector(Alazar.PageAlignedArray{UInt8}, bytes_buf, n_buf)`.

The fields of a `DMABufferVector{S,T}` are:
- `bytes_buf::Int`: The number of bytes per buffer. If there is more than one buffer it
  should be a multiple of Base.Mmap.PAGESIZE. This is enforced in the inner constructor.
- `n_buf::Int`: The number of buffers to allocate.
- `backing::T`: The page-aligned backing array. `S` is `Ptr{eltype(T)}`.

This code may not support 32-bit systems.
"""
struct DMABufferVector{S,T} <: AbstractVector{S}
    bytes_buf::Int
    n_buf::Int
    backing::T

    function DMABufferVector{S,T}(bytes_buf::Integer, n_buf::Integer) where {S,T<:AbstractVector}
        @assert n_buf >= 0
        @assert bytes_buf * n_buf % sizeof(eltype(T)) == 0
        n_buf > 1 && bytes_buf % Base.Mmap.PAGESIZE != 0 &&
            error("bytes per buffer must be a multiple of Base.Mmap.PAGESIZE when ",
                  "there is more than one buffer.")

        backing = T(div(bytes_buf * n_buf, sizeof(eltype(T))))
        Integer(pointer(backing)) % Base.Mmap.PAGESIZE != 0 &&
            error("array type $T is not page-aligned.")

        return new{S,T}(bytes_buf, n_buf, backing)
    end
end
DMABufferVector(::Type{T}, bytes_buf::Integer, n_buf::Integer) where {S, T<:AbstractVector{S}} =
    DMABufferVector{Ptr{S}, T}(bytes_buf, n_buf)

### DMABufferVector interface

bytespersample(buf_array::DMABufferVector{Ptr{T}}) where {T} = sizeof(T)
sampletype(buf_array::DMABufferVector{Ptr{T}}) where {T} = T

### AbstractArray methods

Base.size(dma::DMABufferVector) = (dma.n_buf,)
Base.IndexStyle(::Type{<:DMABufferVector}) = Base.IndexLinear()
Base.getindex(dma::DMABufferVector, i::Int) =
    pointer(dma.backing) + (i-1) * dma.bytes_buf
Base.length(dma::DMABufferVector) = dma.n_buf

## PageAlignedArray definitions

"""
    mutable struct PageAlignedArray{T,N} <: AbstractArray{T,N}

An `N`-dimensional array of eltype `T` which is guaranteed to have its memory be page-aligned.

This has to be a mutable struct because finalizers are used to clean up the memory allocated
by C calls when there remain no references to the PageAlignedArray object in Julia.
"""
mutable struct PageAlignedArray{T,N} <: AbstractArray{T,N}
    backing::Array{T,N}
    addr::Ptr{T}

    PageAlignedArray{T,0}() where {T} = PageAlignedArray{T,0}(())
    PageAlignedArray{T,1}(a::Integer) where {T} = PageAlignedArray{T,1}((a,))
    PageAlignedArray{T,2}(a::Integer, b::Integer) where {T} = PageAlignedArray{T,0}((a,b))
    PageAlignedArray{T,3}(a::Integer, b::Integer, c::Integer) where {T} =
        PageAlignedArray{T,3}((a,b,c))

    function PageAlignedArray{T,N}(dims::NTuple{N,Integer}) where {T,N}
        n = N == 0 ? 1 : reduce(*, dims)
        addr = virtualalloc(sizeof(T) * n, T)
        backing = unsafe_wrap(Array, addr, dims, false)
        array = new{T,N}(backing, addr)
        finalizer(array, x->virtualfree(x.addr))
        return array
    end
end
PageAlignedArray{T}(dims::Integer...) where {T} = PageAlignedArray{T,length(dims)}(dims)
const PageAlignedVector{T} = PageAlignedArray{T,1}

Base.size(A::PageAlignedArray) = size(A.backing)
Base.IndexStyle(::Type{<:PageAlignedArray}) = Base.IndexLinear()
Base.getindex(A::PageAlignedArray, idx...) = A.backing[idx...]
Base.setindex!(A::PageAlignedArray, v, idx...) = setindex!(A.backing, v, idx...)
Base.length(A::PageAlignedArray) = length(A.backing)
Base.unsafe_convert(::Type{Ptr{T}}, A::PageAlignedArray{T}) where {T} =
    Base.unsafe_convert(Ptr{T}, A.backing)

"""
    virtualalloc(size_bytes::Integer, ::Type{T}) where {T}

Allocate page-aligned memory and return a `Ptr{T}` to the allocation. The caller is
responsible for de-allocating the memory using `virtualfree`, otherwise it will leak.
"""
function virtualalloc(size_bytes::Integer, ::Type{T}) where {T}
    @static is_windows() ? begin
        MEM_COMMIT = 0x1000
        PAGE_READWRITE = 0x4
        addr = ccall((:VirtualAlloc, "Kernel32"), Ptr{T},
            (Ptr{Void}, Culonglong, Culong, Culong),
            C_NULL, size_bytes, MEM_COMMIT, PAGE_READWRITE)
    end : @static is_linux() ? begin
        addr = ccall((:valloc, linux_libc), Ptr{T}, (Culonglong,), size_bytes)
    end : @static is_apple() ? begin
        addr = ccall((:valloc, "libSystem.dylib"), Ptr{T}, (Culonglong,), size_bytes)
    end : throw(SystemError())

    addr == C_NULL && throw(OutOfMemoryError())
    return addr::Ptr{T}
end

"""
    virtualfree(addr::Ptr{T}) where {T}

Free memory that has been allocated using `virtualalloc`. Undefined, likely very bad
behavior if called on a pointer coming from elsewhere.
"""
function virtualfree(addr::Ptr{T}) where {T}
    @static is_windows() ? begin
        MEM_RELEASE = 0x8000
        return ccall((:VirtualFree, "Kernel32"), Cint, (Ptr{Void}, Culonglong, Culong),
            addr, 0, MEM_RELEASE)
    end : @static is_linux() ? begin
        return ccall((:free, linux_libc), Void, (Ptr{Void},), addr)
    end : @static is_apple() ? begin
        return ccall((:free, "libSystem.dylib"), Void, (Ptr{Void},), addr)
    end : error("OS not supported")
end


text:
const COST_ITEMS = [
    :income_tax,
    :national_insurance,
    :employers_ni,
    :scottish_income_tax,

    :total_benefits,

    :means_tested_bens,
    :legacy_mtbs,
    :universal_credit,
    :non_means_tested_bens,
    :sickness_illness,
    :scottish_benefits]

const COST_LABELS = [
    "Total Income Tax",
    "Employee's National Insurance",
    "Employer's National Insurance",
    "Scottish Income Tax",

    "Total Benefit Spending",

    "All Means Tested Benefits",
    "Legacy Means-Tested Benefits",
    "Universal Credit",
    "Non Means Tested Benefits",
    "Disability, Sickness-Related Benefits",
    "Scottish Benefits" ]

const MR_LABELS =
    ["Negative or Zero",
     "Under 10%",
     "10-20%",
     "20-30%",
     "30-50%",
     "50-80%",
     "80-100%",
     "Above 100%"]

function extract_incs( d :: DataFrame, targets :: Vector{Symbol}, row = 1 ) :: Vector
    n = length( targets )[1]
    out = zeros(n)
    for i in 1:n
        out[i] = d[row, targets[i]]
    end
    return out
end

function f0( n :: Number ) :: String
	format(n, commas=true, precision=0 )
	# , autoscale=:finance
end

function fp( n :: Number ) :: String
	format(n, precision=2, signed=true )
	# , autoscale=:finance
end

function f2( n :: Number ) :: String
	format(n, commas=true, precision=2 )
	# , autoscale=:finance
end

function costs_dataframe(  incs1 :: DataFrame, incs2 :: DataFrame ) :: DataFrame
    pre = extract_incs( incs1, COST_ITEMS ) ./ 1_000_000
    post = extract_incs( incs2, COST_ITEMS ) ./ 1_000_000
    diff = post-pre
    return DataFrame( Item=COST_LABELS, Before=pre, After=post, Change=diff )
end

function mr_dataframe( mr1::Histogram, mr2::Histogram, mean1::Real, mean2 :: Real ) :: DataFrame
    println( "mr1.weights=$(mr1.weights) mean1=$mean1")
    change = mr2.weights - mr1.weights
    df = DataFrame( Item=MR_LABELS, Before=mr1.weights, After=mr2.weights, Change=change)
    mchange = mean2 - mean1
    push!( df, (Item="Average METR", Before=mean1, After=mean2, Change=mchange ))
    return df
end

function ineq_dataframe( ineq1 :: InequalityMeasures, ineq2 :: InequalityMeasures )
    names = ["Gini", "Palma", "Atkinson(ϵ=0.5)", "Atkinson(ϵ=1)", "Atkinson(ϵ=2)", "Hoover"]
    v1 = [ineq1.gini, ineq1.palma, ineq1.atkinson[2], ineq1.atkinson[4], ineq1.atkinson[8], ineq1.hoover] .* 100
    v2 = [ineq2.gini, ineq2.palma, ineq2.atkinson[2], ineq2.atkinson[4], ineq2.atkinson[8], ineq2.hoover] .* 100
    diff = v2 -v1
    return DataFrame( Item=names, Before=v1, After=v2, Change=diff)
end

function pov_dataframe(
    pov1 :: PovertyMeasures,
    pov2 :: PovertyMeasures,
    ch1 :: GroupPoverty,
    ch2 :: GroupPoverty )
    println( "got child poverty[1] as $(ch1)")
    names = ["Headcount (All)", "Child Poverty", "Gap", "FGT(α=2)", "Watts", "Sen", "Shorrocks"]
    # child povs already %s ..
    v1 = [pov1.headcount, ch1.prop, pov1.gap, pov1.foster_greer_thorndyke[5], pov1.watts, pov1.sen, pov1.shorrocks]  .* 100
    v2 = [pov2.headcount, ch2.prop, pov2.gap, pov2.foster_greer_thorndyke[5], pov2.watts, pov2.sen, pov2.shorrocks]  .* 100
    diff = v2 - v1
    return DataFrame( Item=names, Before=v1, After=v2, Change=diff)
end

function detailed_cost_dataframe( inc1 :: DataFrame, inc2 :: DataFrame )
    i1 = collect(values(inc1[1,1:98])) ./ 1_000_000
    c1 = collect(values(inc1[2,1:98]))
    i2 = collect(values(inc2[1,1:98])) ./ 1_000_000
    c2 = collect(values(inc2[2,1:98])) # fixme parameterise 98
    dc = c2 - c1
    di = i2 - i1
    names = pretty.(collect((keys(inc1[end,1:98]))))
    return  DataFrame(
        :Item=>names,
        :value1=>i1,
        :count1=>c1,
        :value2=>i2,
        :count2=>c2,
        :dval => di,
        :dcount => dc)
end

ext: jl


text:
using ThreeBodyDecay
using Parameters
using Plots
using PartialWaveFunctions
theme(:wong)

urand(x1,x2) = x1+(x2-x1)*rand()
const mΞcc = 3.6212
const mΞc  = 2.46794
const mΞc′  = 2.5784
const mπ   = 0.14
const mπ0  = 0.135
const mγ   = 0.0
#
function four_body_decay(vars; ms)
    @unpack s,σ1,cosθ1,cosθ23,ϕ23 = vars
    @unpack m0,m1,m2,m3,m4 = ms
    p0 = [0,0,0,m0]
    p123, p4 = four_vectors_in_binary_decay(p0,1,0; m1sq=s, m2sq=m4^2)
    p23, p1 = four_vectors_in_binary_decay(p123,cosθ1,0; m1sq=σ1, m2sq=m1^2)
    p2,p3 = four_vectors_in_binary_decay(p23,cosθ23,ϕ23; m1sq=m2^2, m2sq=m3^2)
    return p1,p2,p3,p4
end
#
phsp(s,σ;ms) = sqrt(λ(ms.m0^2,s,ms.m4^2)*λ(s,σ,ms.m1^2)*λ(σ,ms.m2^2,ms.m3^3))/(σ*s)
#
function rand_Ξcpipi_phsp(; m_miss)
    #
    ms = (m0 = mΞcc, m1 = mπ, m2 = mπ, m3 = mΞc, m4 = m_miss)
    #
    s = urand((ms.m1+ms.m2+ms.m3)^2,(ms.m0-ms.m4)^2)
    σ1 = urand((ms.m2+ms.m3)^2,(ms.m0-(ms.m4+ms.m1))^2)
    σ1 > (√s-ms.m1)^2 && return rand_Ξcpipi_phsp(; m_miss=m_miss)
    cosθ1,cosθ23,ϕ23 = urand(-1,1),urand(-1,1),urand(-π,π)
    vars = (s=s,σ1=σ1,cosθ1=cosθ1,cosθ23=cosθ23,ϕ23=ϕ23)
    #
    p1,p2,p3,p4 = four_body_decay(vars; ms=ms)
    !(invmasssq(p1+p2+p3+p4) ≈ ms.m0^2) && error("something is wrong with s")
    !(invmasssq(p1+p2+p3) ≈ s) && error("something is wrong")
    !(invmasssq(p2+p3) ≈ σ1) && error("something is wrong with σ")
    (m=sqrt(invmasssq(p1+p2+p3)), w=phsp(s,σ1;ms=ms))
end
#
function rand_Ξcpipi_Ξc′(;mΞc′=mΞc′)
    ms = (m0 = mΞcc, m2 = mΞc, m3 = mγ, m1 = mπ, m4 = mπ)
    #
    s = urand((ms.m1+ms.m2+ms.m3)^2,(ms.m0-ms.m4)^2)
    σ1 = mΞc′^2
    √s-ms.m1 < √σ1 && return rand_Ξcpipi_Ξc′(;mΞc′=mΞc′)
    cosθ1,cosθ23,ϕ23 = urand(-1,1),urand(-1,1),urand(-π,π)
    vars = (s=s,σ1=σ1,cosθ1=cosθ1,cosθ23=cosθ23,ϕ23=ϕ23)
    #
    p1,p2,p3,p4 = four_body_decay(vars; ms=ms)
    (m=sqrt(invmasssq(p1+p2+p4)), w=phsp(s,σ1;ms=ms))
end
#
function rand_Ξcpipi_ω(;mω=0.78)
    ms = (m0 = mΞcc, m2 = mπ, m3 = mπ, m1 = mπ0, m4 = mΞc)
    #
    s = mω^2
    σ1 = urand((ms.m2+ms.m3)^2,(√s-ms.m1)^2)
    cosθ1,cosθ23,ϕ23 = urand(-1,1),urand(-1,1),urand(-π,π)
    vars = (s=s,σ1=σ1,cosθ1=cosθ1,cosθ23=cosθ23,ϕ23=ϕ23)
    #
    p1,p2,p3,p4 = four_body_decay(vars; ms=ms)
    (m=sqrt(invmasssq(p2+p3+p4)), w=phsp(s,σ1;ms=ms))
end

@time let Nev = 1e6, bins = range(2.8,3.7,length=300)
    plot(leg=:left, title="feed-down Ξcc")
    calv = [rand_Ξcpipi_Ξc′() for _ in 1:Nev]; w = getindex.(calv,:w)
    stephist!(getindex.(calv,:m), weights=w/sum(w)*Nev , bins=bins, α=0.7, lab="Ξc′ (Ξcγ)")
    calv = [rand_Ξcpipi_ω() for _ in 1:Nev]; w = getindex.(calv,:w)
    stephist!(getindex.(calv,:m), weights=w/sum(w)*Nev, bins=bins, α=0.7, lab="ω (π⁻π⁺π⁰)")
    calv = [rand_Ξcpipi_phsp(m_miss=mπ0) for _ in 1:Nev]; w = getindex.(calv,:w)
    stephist!(getindex.(calv,:m), weights=w/sum(w)*Nev, bins=bins, α=0.7, lab="ph. sp.(miss π0)")
    calv = [rand_Ξcpipi_phsp(m_miss=mγ) for _ in 1:Nev]; w = getindex.(calv,:w)
    stephist!(getindex.(calv,:m), weights=w/sum(w)*Nev, bins=bins, α=0.7, lab="ph. sp.(miss γ)")
    vline!([mΞcc], c=:black, lab="mΞcc", xlab="m(Ξcππ) (GeV)")
end
savefig(joinpath("feeddown_Xicc.pdf"))

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
#####
##### Convert struct to dictionary
#####
using Base.Iterators

export @struct_2_dict

"""
    add_leaves!(D_full::Dict,
                     s::T,
                     is_leaf::F,
                     stop_recursion_type::SRT,
                     _basename::S) where {S,T,F<:Function,SRT}

Add leaves of data structure `s` to dictionary `D_full`
"""
function add_leaves!(D_full::Dict,
                     s::T,
                     is_leaf::F,
                     _basename::S) where {S,T,F<:Function,SRT}
  fns = fieldnames(T)
  for fn in fns
    prop = getproperty(s, fn)
    if triggers(is_leaf(s, prop))
      D_full[_basename*"."*string(fn)] = prop
    else
      add_leaves!(D_full, prop, is_leaf, _basename*"."*string(fn))
    end
  end
end

function struct_2_dict(s, obj_name::AbstractString; ft::AbstractFilterTypes=EntireStruct())
  D_full = Dict()
  is_leaf(s, prop) = __is_leaf(s, prop, ft)
  add_leaves!(D_full, s, is_leaf, obj_name)
  return D_full
end

"""
    @struct_2_dict(s::T) where {T}

Convert struct to a `Dict` with
a hierarchical naming scheme.
"""
macro struct_2_dict(s)
  obj_name = string(s)
  return esc(quote
               FreeParameters.struct_2_dict($s, $obj_name)
             end)
end


text:
function sigmoid(x::Float64)::Float64
    return (MathConstants.e^x)/(1 + MathConstants.e^x)
end

function sigmoid_p(x::Float64)::Float64
    return (MathConstants.e^x)/(1 + 2MathConstants.e^x + MathConstants.e^(2x))
end

function MSE(a::Vector{Float64},l::Vector{Float64})::Float64
    return 0.5*(sum((a.-l).^2))
end

function MSE_p(a::Vector{Float64},l::Vector{Float64})::Float64
    return sum(a.-l)
end

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_name: ian-double-u/NNfromScratch


text:
# test that the SdErgm interface works for a concrete model<:SdErgm

using Pkg
Pkg.activate(".")
Pkg.instantiate()

using Test
using ScoreDrivenERGM
import ScoreDrivenERGM:StaticNets,DynNets, Utilities
using LinearAlgebra
using SharedArrays
using PyPlot
using Statistics
using StatsBase
using DataStructures

# model = DynNets.SdErgmDirBin0Rec0_pmle("HESS_D")
model = DynNets.SdErgmDirBin0Rec0_pmle(scoreScalingType ="FISH_D")
model_mle = DynNets.SdErgmDirBin0Rec0_mle(scoreScalingType="FISH_D")
model_2 = model_mle
# model_3 =  DynNets.SdErgmDirBin0Rec0_pmle(scoreScalingType="FISH_D", options = SortedDict("Firth" => true))

# sample dgp
@time begin
N = 50
T =300
quantileVals = [[0.975, 0.025]]
listDgpSettigns = DynNets.list_example_dgp_settings(model_mle)
# dgpSet = listDgpSettigns.dgpSetSD
# dgpSet.opt.A[1] = 1
dgpSet = listDgpSettigns.dgpSetARlow

parDgpT = DynNets.sample_time_var_par_from_dgp(model_mle, dgpSet.type, N, T;  dgpSet.opt...)

A_T = DynNets.sample_ergm_sequence(model_mle, N, parDgpT, 1)[:,:,:,1]

    obsT, vEstSdResPar, fVecT_filt, target_fun_val_T, sVecT_filt, conv_flag, ftot_0, confBandsFiltPar, confBandsPar, errFlag = DynNets.estimate_filter_and_conf_bands(model, A_T, parDgpT=parDgpT, quantileVals, plotFlag=true);
    res = (;  obsT, vEstSdResPar, fVecT_filt, target_fun_val_T, sVecT_filt, conv_flag, ftot_0, confBandsFiltPar, confBandsPar, errFlag)
end

@time begin
    obsT, vEstSdResPar, fVecT_filt, target_fun_val_T, sVecT_filt, conv_flag, ftot_0, confBandsFiltPar, confBandsPar, errFlag = DynNets.estimate_filter_and_conf_bands(model_2, A_T, parDgpT=parDgpT, quantileVals, plotFlag=true);
    res_2 =
    (;  obsT, vEstSdResPar, fVecT_filt, target_fun_val_T, sVecT_filt, conv_flag, ftot_0, confBandsFiltPar, confBandsPar, errFlag)
end

@time begin
    #     obsT, vEstSdResPar, fVecT_filt, target_fun_val_T, sVecT_filt, conv_flag, ftot_0, confBandsFiltPar, confBandsPar, errFlag = DynNets.estimate_filter_and_conf_bands(model_3, A_T, parDgpT=parDgpT, quantileVals, plotFlag=true);
    #     res_3 =
    #     (;  obsT, vEstSdResPar, fVecT_filt, target_fun_val_T, sVecT_filt, conv_flag, ftot_0, confBandsFiltPar, confBandsPar, errFlag)
end

# Single snapshot estimators
begin
fVecT_filt_SS =  DynNets.estimate_single_snap_sequence(model, A_T)

fVecT_filt_SS_2 =  DynNets.estimate_single_snap_sequence(model_2, A_T)

fig, ax =  DynNets.plot_filtered(model, N, fVecT_filt_SS; lineType = ".", lineColor="r", parDgpTIn=parDgpT)
DynNets.plot_filtered(model, N, fVecT_filt_SS_2; lineType = ".", lineColor="b", parDgpTIn=parDgpT, fig=fig, ax=ax)
DynNets.plot_filtered(model, N, fVecT_filt; lineType = "-", parDgpTIn=parDgpT, fig=fig, ax=ax)
DynNets.plot_filtered(model, N, mapslices(rolling_mean,fVecT_filt_SS, dims=2); lineType = "-", parDgpTIn=parDgpT, fig=fig, ax=ax)

end

@show var(parDgpT, dims=2)
@show var(fVecT_filt, dims=2)
@show var(fVecT_filt_SS, dims=2)

# Compare estimators of variances of static parameter's estimators
mvSDUnParEstCovWhite, errFlagEstCov = DynNets.white_estimate_cov_mat_static_sd_par(model_2, N, res_2.obsT, model_2.indTvPar, res_2.ftot_0, res_2.vEstSdResPar)

mvSDUnParEstCovWhite[I(6)]

distribFilteredSDParBoot, filtCovHatSample, errFlagVec , vEstSdUnParBootDist = DynNets.par_bootstrap_distrib_filtered_par(model_2, N, res_2.obsT, model_2.indTvPar, res_2.ftot_0, res_2.vEstSdResPar)

vEstSdUnParNPBootDist = DynNets.non_par_bootstrap_distrib_filtered_par(model_2, N, obsT, model_2.indTvPar, ftot_0; nBootStrap = 100)

figure()
hist(filtCovHatSample')
mapslices( x -> mean(winsor(x, prop=0.005)), filtCovHatSample, dims=2)

using ScikitLearn
@sk_import covariance: MinCovDet

DynNets.conf_bands_given_SD_estimates(model_2, N, res_2.obsT, vEstSdUnPar, res_2.ftot_0, [[0.975, 0.025]]; parDgpT=parDgpT, plotFlag=true, parUncMethod = "WHITE-MLE", offset = 0, mvSDUnParEstCov = 0.3 .*mvSDUnParEstCovWhite, winsorProp=0.05)

DynNets.conf_bands_given_SD_estimates(model_2, N, res_2.obsT, vEstSdUnPar, res_2.ftot_0, [[0.975, 0.025]]; parDgpT=parDgpT, plotFlag=true, parUncMethod = "NPB-MVN", offset = 0, mvSDUnParEstCov = MinCovDet().fit(Utilities.drop_bad_un_estimates(vEstSdUnParNPBootDist)').covariance_, winsorProp=0.05)

DynNets.conf_bands_given_SD_estimates(model_2, N, res_2.obsT, vEstSdUnPar, res_2.ftot_0, [[0.975, 0.025]]; parDgpT=parDgpT, plotFlag=true, parUncMethod = "PB-MVN", offset = 0, mvSDUnParEstCov = MinCovDet().fit(Utilities.drop_bad_un_estimates(vEstSdUnParBootDist)').covariance_, winsorProp=0.05)

DynNets.conf_bands_given_SD_estimates(model_2, N, res_2.obsT, vEstSdUnPar, res_2.ftot_0, [[0.975, 0.025]]; parDgpT=parDgpT, plotFlag=true, parUncMethod = "PB-SAMPLE", offset = 0, sampleStaticUnPar=Utilities.drop_bad_un_estimates(vEstSdUnParBootDist), winsorProp=0 )

if true
parNames = ["w_theta", "B_theta", "A_theta", "w_eta", "B_eta", "A_eta"]
sample1 = vEstSdUnParBootDist
sample1 = Utilities.drop_bad_un_estimates(sample1)

fig, ax = subplots(3,2)
redLines = DynNets.unrestrict_all_par(model_2, res_2.vEstSdResPar)
    if true
        sample1 = mapslices(x -> DynNets.restrict_all_par(model_2, x), sample1, dims=1)
        redLines = res_2.vEstSdResPar
    end

for i = 1:6
    ax[i].hist(sample1[i,:], range=quantile(sample1[i,:], [0.01, 0.99]), 20, alpha = 0.4, density=true)
    ax[i].set_title(parNames[i])
    ylim = ax[i].get_ylim()
    ax[i].vlines(redLines[i], ylim[1], ylim[2], color = "r" )

end
end

plt.figure()
plt.hist(mvSDUnParEstCovParBoot[3, :])

distribFilteredSD, filtCovHatSample,  _, errFlagSample = DynNets.distrib_filtered_par_from_mv_normal(model_2, N, res_2.obsT, model_2.indTvPar, res_2.ftot_0, res_2.vEstSdUnPar, mvSDUnParEstCov)

# nBootStrap = 150
# vEstSdUnParBootDist = SharedArray(zeros(3*sum(model.indTvPar), nBootStrap))

# @time Threads.@threads for k=1:nBootStrap
#     vEstSdUnParBootDist[:, k] = rand(1:T, T) |> (inds->(inds |> x-> DynNets.estimate(model, N, res.obsT; indTvPar=model.indTvPar, ftot_0 = res.ftot_0, shuffleObsInds = x) |> x-> getindex(x, 1) |> x -> DynNets.array_2_vec_all_par(model, x, model.indTvPar))) |> x -> DynNets.unrestrict_all_par(model, x)
# end

# @time begin
#     obsT, vEstSdResPar, ftot_0, fVecT_filt, confBandsFiltPar, confBandsPar, errFlag = DynNets.estimate_filter_and_conf_bands(DynNets.SdErgmDirBin0Rec0_pmle("FISH_D"), A_T, parDgpT=parDgpT, quantileVals, plotFlag=true);
#     resFish = (;obsT, vEstSdResPar, fVecT_filt, target_fun_val_T, sVecT_filt, conv_flag, ftot_0)
# end

# IL PROBLEMA € NELLA MATRICE DI VARIANZA COVARIANZA STIMATA CON WHITE. QUELLA NON HA LA GIUSTA MAGNITUDE, quando i parametri time varying variano molto

# @time fVecT_filt, confBandsFiltPar, confBandsPar, errFlag, mvSDUnParEstCov, distribFilteredSD = DynNets.conf_bands_given_SD_estimates(model, N, res.obsT, res.vEstSdResPar, res.ftot_0, quantileVals; indTvPar = model.indTvPar, parDgpT=parDgpT, plotFlag=true, parUncMethod = "NPB-MVN")

# correctly specified filter

# estimate

# store the sufficient statistics and change statistics in R


text:
# Estimate the best-case resolution (in nanoseconds) of benchmark timings based
# on the system clock.
#
# Arguments:
#
#     n_samples::Integer: The number of times we call the system clock to
#         estimate its resolution. Defaults to 10,000 calls.
#
# Returns:
#
#     min_Δt::UInt: The estimated clock resolution in nanoseconds.
#
# TODO:
#
#     This function is known to not work on Windows because of the behavior
#     of `time_ns()` on that platform.

function estimate_clock_resolution(n_samples::Integer = 10_000)
    min_Δt = typemax(UInt)

    for _ in 1:n_samples
        t1 = Base.time_ns()
        t2 = Base.time_ns()
        # On Linux AArch64 it seems that t1 and t2 could be the same sometimes
        while t2 <= t1
            t2 = Base.time_ns()
        end
        Δt = t2 - t1
        min_Δt = min(min_Δt, Δt)
    end

    min_Δt
end


text:
## GeomUtils.jl --- constructs for geometry

export distance, anglespan, dihedral
export centroid, RMSD

"""
    distance(a, b) -> Real

Calculate Euclidean distance from `a` to `b`.
"""
distance(a, b) = norm(a-b)

"""
    anglespan(a, b, c) ∈ [0, π]

Calculate smaller angle between two vectors ``AB`` and ``BC`` meeting at ``B``.
"""
function anglespan(a, b, c)
    Δ1 = a - b; n1 = norm(Δ1)
    Δ2 = c - b; n2 = norm(Δ2)
    (n1 == 0 || n2 == 0) && error("two positions overlap")
    cos_ = Δ1⋅Δ2 / norm(Δ1) / norm(Δ2)
    cos_ >= 1 && return 0        ## handle "overflow?" cases where cos_ ∉ [-1,1]
    cos_ <= -1 && return π
    acos(cos_)
end

"""
    dihedral(a, b, c, d) ∈ [-π, π]

Calculate right handed dihedral (torsion) angle between the plane ``ABC`` and
``BCD`` intersecting along ``BC``. ``ABC`` is determined by the first three
arguments as points, while ``BCD`` is determined by the last three.
"""
function dihedral(a, b, c, d)
    Δ1 = a - b
    Δ2 = b - c; n2 = norm(Δ2)
    Δ3 = d - c
    n2 == 0 && error("torque axis is undefined; `b`, `c` overlap")

    Δ23 = Δ2×Δ3
    Δ21 = Δ2×Δ1
    (norm(Δ23) == 0 || norm(Δ21) == 0) && error("some points overlap/colinear")

    atan(Δ23⋅(Δ21×Δ2) / n2, Δ23⋅Δ21)
end

"""
    centroid(collection) -> eltype(collection)

Get the centroid of a collection of points. The preferred point representation
is a subtype of `StaticVector`. This centroid is an _unweighted_ average of
positions.
"""
function centroid(collection::Vector{StaticVector})
    sm_ = sum(collection)
    sm_ / length(collection)
end

"""
    centered(collection) -> centered positions

Get a copy of the centered positions in collection. The preferred point
representation is a subtype of `StaticVector`.

Equivalent to `collection .- centroid(collection)`.
"""
centered(collection) = collection .- centroid(collection)

"""
    RMSD(cA, cB) -> Real

Calculate root-mean-square deviation of points in `cA` from those in `cB`.
Preferred point representation is a subtype of `StaticVector`.
"""
function RMSD(cA::Vector{M}, cB::Vector{M}) where M <: StaticVector
    len = length(cA)
    @assert length(cB) == len "arguments don't have the same length"
    msd = sum(cA, cB; init=zero(eltype(M))) do (a, b)
        distance(a, b)
    end / len

    sqrt(msd)
end


text:
module SFrames

export head, tail, materialize, ismaterialized, sample, dropna, unpack, save, load, SArray, SArrayTyped, nummissing, clip, clip_lower, clip_upper, fillna, count_ngrams, count_words, FlexibleType, SFrame, column_names, read_csv, unpack, pack_columns, stack, unstack, addrownumber, topk, Agg, groupby

using Cxx

"""
`head(s::T, n) -> T`

Returns the first `n` elements of `s`.

`s` can be either an `SFrame` or an `SArray`.
"""
function head end
function tail end
function materialize end
function ismaterialized end
function sample end
function dropna end
function unpack end
function save end
function load end

include("Util.jl")
include("FlexibleType.jl")
include("SArray.jl")
include("SFrame.jl")
include("Agg.jl")

# import .SFrameMod: SFrame
# import .SArrayMod: SArray
# import .FlexibleTypeMod: FlexibleType

using .FlexibleTypeMod
using .SArrayMod
using .SFrameMod

export SArray, FlexibleType

function __init__()
    const SFRAME_PATH = get(ENV, "SFRAME_PATH",
        joinpath(homedir(), "SFrame"))
    println("Assuming SFrame is in $SFRAME_PATH")
    Libdl.dlopen(
        joinpath(SFRAME_PATH, "debug/oss_src/unity/libunity_shared.dylib"), Libdl.RTLD_GLOBAL)

    map(path->addHeaderDir(joinpath(SFRAME_PATH,path), kind=C_System), [
        "oss_src",
        "deps/local/include",
    ])

    map(cxxinclude, [
        "unity/lib/gl_sframe.hpp",
        "unity/lib/gl_sgraph.hpp",
        "unity/lib/gl_sarray.hpp",
        "flexible_type/flexible_type.hpp",
        "flexible_type/flexible_type_base_types.hpp"
    ])

    cxx"""
    #include <iostream>
    #include <ostream>
    #include <sstream>
    #include <string>
    #include <vector>
    #include <map>
    #include <utility>

    using namespace graphlab;
    using namespace std;
    """
end

end # module


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ f6d978ea-b349-11eb-309e-f5a517efc907
begin
	using LinearAlgebra
	using Plots
	using QuadGK
end

# ╔═╡ 99268461-9e3f-4685-b530-e24febfb539e
#=
IDEA OF MULTIGRID :

	G is a mesh of a 1D domain.
	G1 is another, but coarser.
	G1's operator has higher frequency eigenvalues as opposed to  G.
	Use G1 to weed out these values.
	Use smoothers alternatively.
=#

# ╔═╡ b422b21d-b9d3-4aad-b29b-d0669331abc8

begin
	e = 5; #  Tune e for different basis functions
	N = 10;
	Ni = 100;
	xmin = 0;
	xmax = 1;
	x = LinRange(xmin,xmax,N);
	xi = LinRange(x[e-1],x[e+1],Ni);
	h = x[2] - x[1];
end

# ╔═╡ b11596d4-ee5d-41ec-af8f-96734ba2f579
# This shows a heirarchy of grids
grid1 = plot(x,zeros(length(x)),marker="*")

# ╔═╡ c1463b68-12e6-4e8c-9902-c2cc7e6caa91
begin
	# Plot the basis functions
	y = zeros(length(xi));
	for i = 1:length(xi);
		if(xi[i]>=xi[1] && xi[i]<=x[e])
			y[i] = 1/h * (xi[i]-x[e-1]);
		elseif (xi[i]>x[e] && xi[i]<x[e+1])
			y[i] = 1/h * (x[e+1]-xi[i]);
		end
	end
	plot!(xi,y,title="Linear Basis Functions");
end

# ╔═╡ d5da54e6-0a4f-4f43-a7ab-8791574b32c6
# Now we use the linear basis functions and obtain an approximation of the Laplacian operator.

begin
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	function θ(x,i,N)
		n = length(x);
		y = zeros(n);
		if i==1
			for k=1:N
				if x[k]<= x[i+1]
					y[k] = (N-1)*(x[i+1]-x[k]);
				end
			end
		elseif i==N
			for l = 1:N
				if x[l] >= x[i-1]
					y[l] = (N-1) * (x[l]-x[i]);
				end
			end
		else
			for j=2:n-1
				if(x[j]>=x[i-1] && x[j]<=x[i])
					y[j] = (N-1) * (x[j]-x[i-1]);
				elseif(x[j]>=x[i] && x[j]<=x[i+1])
					y[j] = (N-1) * (x[i+1]-x[j]);
				end
			end
		end
		return y;
	end
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	function θ_scalar(x,i,xi,N)
			if(x>=xi[i-1] && x<=xi[i])
				y = (N-1) * (x-xi[i-1]);
			elseif(x>=xi[i] && x<=xi[i+1])
				y = (N-1) * (xi[i+1]-x);
			else
				y=0;
			end
		return y;
	end
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	function stiffness_matrix(n)
		md = 2.0*ones(n-2);
		ld = -1.0*ones(n-3);
		ud = -1.0*ones(n-3);
		return (n-1)*Tridiagonal(ld,md,ud);
	end
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	function intergrate(x,i,N)
		domain = LinRange(0,1,N+2);
		return sin(pi*x) * θ_scalar(x,i,domain,N);
	end

	function compute_rhs(n)
		n = n-2;
		integral = zeros(n);
		for i=1:length(integral)
		integral[i] = 1.5*quadgk(q->intergrate(q,i+1,n),0,1,rtol=1e-8)[1];
		end
		return integral;
	end
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
end

# ╔═╡ a5f1a911-3f71-4630-a872-99514d52e5f3
begin
	nk = 12;
	a = LinRange(xmin,xmax,nk);
	z = θ(a,3,nk);
	plot(a,z)
end

# ╔═╡ feba23c9-1b66-400e-9ffb-0dc5fc65896f
begin
	# Get the stiffness matrix and the RHS : The RHS is evaluated and plotted.
	n = 7;
	test_A = stiffness_matrix(n)
	test_b = compute_rhs(n);
	plot(test_b)
	#test_A
	#test_b
end

# ╔═╡ b8b7b670-6039-4bd9-9f95-8652f6e03528
md"""

	#NESTED MULTIGRID

	Start on a coarse grid h_0

	Iterate using a smoother and remove one range of high frequencies.

	Interpolate solution to a finer grid - h_1

	Start step 1 with a smaller time step h_1.
"""

# ╔═╡ a217bef6-3c9f-4d12-8b3f-eb937987c709
# Description

md"""

To implement a nested multigrid method.

We would need access to a smoother. The kind of smoother to choose depends on the nature of the matrix A that we solve.

For the finite element approximation of the poisson equation we have a matrix A that is not diagonally dominant.

So solve the problem we need a method for which convergence could be proven.

In comes the Gauss Seidel method.

"""

# ╔═╡ f2a7eb59-3617-47fe-be19-f6873d0f3c56
function gauss_seidel_relaxation(A,b,x0,max_iterations)
	x = x0;
	res = b - A*x;
	for iter = 1:max_iterations
		for i = 1:length(x)
			x[i] = (1/A[i,i])*(b[i] - dot(A[i,:],x) + A[i,i]*x[i]);
		end
		res = b - A*x;
	end
	return x
end

# ╔═╡ abf72f94-406b-41a6-bbdd-fa3ac2919531
# To perform a nested iteration
begin
	n_nodes = 7;
	A = stiffness_matrix(n_nodes);
	b = compute_rhs(n_nodes);
	x0 = zeros(n_nodes-2);

	y_appx = gauss_seidel_relaxation(A,b,x0,10);
	y_appx = vcat(0,y_appx,0);
	println(y_appx)

	x_nodes = LinRange(0,1,100);
	y_actual = zeros(length(x_nodes));
	for i=1:length(x_nodes)
		y_actual[i] = 1/pi^2 * sin(pi*x_nodes[i]);
	end
	plot(x_nodes,y_actual,label="Actual Solution");
	plot!(LinRange(0,1,length(y_appx)),y_appx,label="Approximate solution");
end

# ╔═╡ Cell order:
# ╠═f6d978ea-b349-11eb-309e-f5a517efc907
# ╠═99268461-9e3f-4685-b530-e24febfb539e
# ╠═b422b21d-b9d3-4aad-b29b-d0669331abc8
# ╠═b11596d4-ee5d-41ec-af8f-96734ba2f579
# ╠═c1463b68-12e6-4e8c-9902-c2cc7e6caa91
# ╠═d5da54e6-0a4f-4f43-a7ab-8791574b32c6
# ╠═a5f1a911-3f71-4630-a872-99514d52e5f3
# ╠═feba23c9-1b66-400e-9ffb-0dc5fc65896f
# ╟─b8b7b670-6039-4bd9-9f95-8652f6e03528
# ╟─a217bef6-3c9f-4d12-8b3f-eb937987c709
# ╠═f2a7eb59-3617-47fe-be19-f6873d0f3c56
# ╠═abf72f94-406b-41a6-bbdd-fa3ac2919531


text:
module MathematicalFunctions

export finv, grad

function finv end
function grad end

end

max_stars_repo_name: ettersi/MathematicalFunctions

max_issues_repo_name: ettersi/MathematicalFunctions

max_forks_repo_name: ettersi/MathematicalFunctions


text:
export
    ZoomNode

mutable struct ZoomNode <: Ranger.AbstractNode
    base::Nodes.NodeData
    zoom::Math.ZoomTransform

    # Collection of nodes.
    children::Array{Ranger.AbstractNode,1}

    # State management
    mouse_window_position::Geometry.Point{Float64}
    zoom_point::Geometry.Point{Float64}
    wheel_direction::Int32   # 0 = not active, 1 = zoom-in, -1 = zoom-out

    function ZoomNode(world::Ranger.World, name::String, parent::Ranger.AbstractNode)
        o = new()

        o.base = Nodes.NodeData(name, parent, world)
        o.zoom = Math.ZoomTransform{Float64}()
        o.children = Array{Ranger.AbstractNode,1}[]

        o.mouse_window_position = Geometry.Point{Float64}()
        o.zoom_point = Geometry.Point{Float64}()
        o.wheel_direction = Int32(0)

        o
    end
end

function Nodes.draw(node::ZoomNode, context::Rendering.RenderContext)
    # Zoom node doesn't need to render but it does require the context
    # for space mapping.

    handle_movement(node)

    handle_zoom(node)

end

function handle_movement(node::ZoomNode)
    # Map window/device-space to zoom-space
    Nodes.map_device_to_node!(node.base.world,
        Int32(node.mouse_window_position.x), Int32(node.mouse_window_position.y),
        node,
        node.zoom_point)

    Custom.set_focal_point!(node, node.zoom_point.x, node.zoom_point.y)
end

function handle_zoom(node::ZoomNode)
    if node.wheel_direction == 1
        Custom.zoom_in!(node)
    elseif node.wheel_direction == -1
        Custom.zoom_out!(node)
    end

    node.wheel_direction = 0
end
# --------------------------------------------------------
# Grouping
# --------------------------------------------------------
function Nodes.get_children(node::ZoomNode)
    node.children
end

# --------------------------------------------------------
# Transforms
# --------------------------------------------------------
# Zoom node manages its own transform differently.
function Nodes.calc_transform!(node::ZoomNode)
    if Nodes.is_dirty(node)
        Math.update!(node.zoom)
        Nodes.set_dirty!(node, false)
    end

    node.zoom.transform
end

# --------------------------------------------------------
# Zooming
# --------------------------------------------------------
function set_position!(node::ZoomNode, x::Float64, y::Float64)
    Math.set_position!(node.zoom, x, y)
    Nodes.ripple_dirty!(node, true);
end

function set_focal_point!(node::ZoomNode, x::Float64, y::Float64)
    Math.set_zoom_at!(node.zoom, x, y)
    Nodes.ripple_dirty!(node, true);
end

function zoom_by!(node::ZoomNode, delta_x::Float64, delta_y::Float64)
    Math.zoom_by!(node.zoom, delta_x, delta_y)
    Nodes.ripple_dirty!(node, true);
end

function translate_by!(node::ZoomNode, delta_x::Float64, delta_y::Float64)
    Math.translate_by!(node.zoom, delta_x, delta_y)
    Nodes.ripple_dirty!(node, true);
end

function zoom_in!(node::ZoomNode)
    Math.zoom_by!(node.zoom, 0.1, 0.1)
    Nodes.ripple_dirty!(node, true);
end

function zoom_out!(node::ZoomNode)
    Math.zoom_by!(node.zoom, -0.1, -0.1)
    Nodes.ripple_dirty!(node, true);
end

# --------------------------------------------------------
# Events
# --------------------------------------------------------
function Nodes.io_event(node::ZoomNode, event::Events.MouseMotionEvent)
    Geometry.set!(node.mouse_window_position, Float64(event.x), Float64(event.y))
end

function Nodes.io_event(node::ZoomNode, event::Events.MouseWheelEvent)
    # Mouse.y:
    # Forward = 1 = zoom in
    # Toward = -1 = zoom out
    node.wheel_direction = event.y
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
function loosen(problem; percent=.1, id_increment=1000, label::Char='?')
    new_lower_bounds = []
    for lb in problem.lower_bounds
        new_lb = convert(Int, floor(lb[2] * percent))
        push!(new_lower_bounds, (lb[1], new_lb))
    end

    MDMKP.MDMKP_Prob(problem.objective, problem.upper_bounds, new_lower_bounds,
        MDMKP.Problem_ID(
            problem.id.id + id_increment,
            problem.id.dataset,
            problem.id.instance,
            problem.id.case,
            problem.id.tightness,
            problem.id.n_vars,
            problem.id.n_demands,
            problem.id.n_dimensions,
            problem.id.mixed_obj,
            true,
            '?'
        ))
end

function decide(problem)
    id = problem.id
    if id.loosened
        if id.n_dimensions < 20
            'A'
        else
            if id.n_vars < 175
                'A'
            else
                if id.tightness >= .625
                    'A'
                else
                    if id.mixed_obj
                        'A'
                    else
                        'C'
                    end
                end
            end
        end
    else
        if id.dataset == 7 && in(id.case, [3, 6])
            'D'
        elseif id.n_demands >= 12.5
            'C'
        else
            if id.n_dimensions < 7.5
                if !id.mixed_obj
                    'A'
                else
                    if id.n_vars < 175
                        'A'
                    else
                        'B'
                    end
                end
            else
                if id.n_vars < 175
                    'A'
                else
                    'B'
                end
            end
        end
    end
end

function set_decision(problem)
    decision = decide(problem)

    MDMKP.MDMKP_Prob(
        problem.objective,
        problem.upper_bounds,
        problem.lower_bounds,

        MDMKP.Problem_ID(
            problem.id.id,
            problem.id.dataset,
            problem.id.instance,
            problem.id.case,
            problem.id.tightness,
            problem.id.n_vars,
            problem.id.n_demands,
            problem.id.n_dimensions,
            problem.id.mixed_obj,
            problem.id.loosened,
            decision
        ))
end

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
export linop, linop_oriig, fallback_demo
using BenchmarkTools
using Random

"""
    fallback_demo()

```julia-repl
julia> n = 10000
10000
julia> Random.seed!(0); Asym = Symmetric(sprandn(n, n, 0.01)); b = randn(n); nnz(Asym.data)
1000598

julia> @benchmark linop_orig(Asym) # abstractarray fallback
BenchmarkTools.Trial:
  memory estimate:  16 bytes
  allocs estimate:  1
  --------------
  minimum time:     2.406 s (0.00% GC)
  median time:      2.409 s (0.00% GC)
  mean time:        2.439 s (0.00% GC)
  maximum time:     2.503 s (0.00% GC)
  --------------
  samples:          3
  evals/sample:     1

julia> @benchmark linop(Asym) # sparsearray fallback in action
BenchmarkTools.Trial:
  memory estimate:  15.32 MiB
  allocs estimate:  8
  --------------
  minimum time:     10.980 ms (0.00% GC)
  median time:      11.669 ms (3.40% GC)
  mean time:        11.807 ms (3.59% GC)
  maximum time:     64.878 ms (82.49% GC)
  --------------
  samples:          423
  evals/sample:     1

julia> @benchmark linop(\$(sparse(Asym))) # only sparse operation
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     957.077 μs (0.00% GC)
  median time:      1.001 ms (0.00% GC)
  mean time:        1.008 ms (0.00% GC)
  maximum time:     1.091 ms (0.00% GC)
  --------------
  samples:          4923
  evals/sample:     1
```
"""
function fallback_demo(n::Integer=10000)
    Random.seed!(0)
    Asym = Symmetric(sprandn(n, n, 0.01))
    b1 = @benchmark linop_orig($Asym) # abstractarray fallback
    b2 = @benchmark linop($Asym) # sparsearray fallback in action
    b3 = @benchmark linop($(sparse(Asym))) # only sparse operation
    b1, b2, b3
end

# This file demonstrates how to change the abstractarray fallback to a
# sparse fallback.
#
#
# Sample function with abstractarray fallback
function linop_orig(A::AbstractMatrix)
    n, m = size(A)
    sum = zero(eltype(A))
    for j = 1:n
        for i = 1:m
            sum += A[i,j]
        end
    end
    sum
end

# An implementation for sparse matrices exists
function linop(A::SparseMatrixCSC{T}) where T
    m, n = size(A)
    sum = zero(T)
    nz = nonzeros(A)
    for k = 1:nnz(A)
        sum += nz[k]
    end
    # error("provoke error")
    return sum
end

# An implementation for sparse matrices exists
function linop(A::Symmetric{T,<:SparseMatrixCSC}) where T
    A = A.data
    m, n = size(A)
    sum = zero(T)
    nz = nonzeros(A)
    for k = 1:nnz(A)
        sum += nz[k]
    end
    # error("symmetric error")
    return sum
end

# The original function is modified
function linop(@nospecialize A::AbstractMatrix{T}) where T
    # inflate reduces depth by 1 - try to catch other specific methods
    if iswrsparse(A) && depth(A) > 0
        linop(inflate(A))
    else
        linop_aaf(A)
    end
end

# The original method is renamed to break potential calling loop
function linop_aaf(A::AbstractMatrix)
    n, m = size(A)
    sum = zero(eltype(A))
    for j = 1:n
        for i = 1:m
            sum += A[i,j]
        end
    end
    sum
end

# find all methods with at least one argument of type `AbstractMatrix`
function allmethods()
    meths = Method[]
    for mod = Base.loaded_modules_array()
        for nm in names(mod)
            if isdefined(mod, nm)
                f = getfield(mod, nm)
                if isa(f, Function)
                    append!(meths, methods(f))
                end
            end
        end
    end
    meths
end

parameters(m::Method) = parameters(m.sig)
parameters(u::UnionAll) = parameters(u.body)
parameters(d::DataType) = d.parameters

name(u::UnionAll) = name(u.body)
name(d::DataType) = d.name
name(u::Union) = u

# for one method check if there is a related method with same name and `AbstractMatrix`
# replaced by `AbstractSparseMatrix`
function relatedmethods(m::Method)
    p = parameters(m)
    f = p[1].instance
    an = (AbstractMatrix)
    repl(x) = (x isa Type && an <: x) ? (SparseMatrixCSC{T} where T) : x
    args = Tuple(repl(x) for x in p[2:end])
    # println("which($f, $args)")
    rel = which(f, args)
    rel == m ? nothing : (m, rel)
end


text:
# CUDA related types

import Base: length, size, eltype

export
    CuIn, CuOut, CuInOut

#
# Device pointer
#

# TODO: a custom DevicePtr class would be nice, but Julia doesn't reason through
#       nested immutable classes (CuDeviceArray -> DevicePtr -> Ptr),
#       resulting in boxed accesses in CuDeviceArray.get/setindex
# immutable DevicePtr{T}
#     ptr::Ptr{T}
# end

typealias DevicePtr{T} Ptr{T}

isnull{T}(p::DevicePtr{T}) = (p == 0)

#
# Managed data containers
#

abstract CuManaged{T}

length(i::CuManaged) = length(i.data)
size(i::CuManaged) = size(i.data)
eltype{T}(::CuManaged{T}) = T

type CuIn{T} <: CuManaged{T}
    data::T
end

type CuOut{T} <: CuManaged{T}
    data::T
end

type CuInOut{T} <: CuManaged{T}
    data::T
end

# FIXME: define on CuManaged instead of on each instance?
#eltype{T}(::Type{CuManaged{T}}) = T
eltype{T}(::Type{CuIn{T}}) = T
eltype{T}(::Type{CuOut{T}}) = T
eltype{T}(::Type{CuInOut{T}}) = T
eltype{T<:CuManaged}(::Type{T}) =
    error("missing eltype definition for this managed type")


text:
module Server

using DotEnv
DotEnv.config()

include("utils.jl")
include("controller/requests.jl")
include("types.jl")
include("controller/network.jl")
include("controller/commands.jl")

using .Types, .Requests, Sockets, Dates

function is_valid_command(command::String)
    if (!startswith(command, "/"))
        return (false)
    end
    true
end

function parse_line!(command::String, conn::IO)
    if (is_valid_command(command) == false)
        return ([])
    end
    command = replace(command, "/" => "")
    command = split(command, " ")
    command
end

function wait_client(conn::IO, storage)
    line = readline(conn)
    parsed_line = parse_line!(line, conn)

    if (size(parsed_line)[1] > 0)
        Commands.exec_command(parsed_line, storage, conn)
    end
    Network.fancy_write(storage, conn)
    if (length(line) <= 0)
        return
    end
    Network.broadcast_channel(storage, line, conn)

    # Store logs in files
    Utils.log(conn, line)
end

function add_default_client(storage, ip_addr::String, conn::IO)
    storage.active_clients[ip_addr] = Client(ip_addr, conn)
end

function serve(port::Int)
    storage = Storage(listen(IPv4(0), port), SQL())
    Network.init_lobby!(storage)

    print("Server listening on port $port\n")
    while true
        conn = accept(storage.listener)
        ip_addr = string(first(getpeername(conn)))

        if (Network.is_connected(storage, ip_addr))
            write(conn, "Someone is already using this ip address ($ip_addr)\n")
            continue
        end

        write(conn, Types.welcome_msg)
        add_default_client(storage, ip_addr, conn)
        Network.fancy_write(storage, conn)

        @async begin
            try
                while isopen(conn)
                    wait_client(conn, storage)
                end
            catch err
                Network.disconnect!(storage, conn)
            end
        end
    end
end

export serve

end # Server


text:
"""
Create an experiment
"""
function create_experiment(api::BaseAPI=BASE_API)::Int
	response = send_request(api, "experiments"; verb="POST")
	# JSON RESPONSE:
	# {"result":"success","id":147}
	to_json(response)["id"]
end

"""
Create an item where `id` is the `items_types` id.
"""
function create_item(id::T, api::BaseAPI=BASE_API) where T <: Integer
	response = send_request(api, "items/$(id)"; verb="POST")
	to_json(response)
end

"""
Create a template
"""
function create_template(api::BaseAPI=BASE_API)
	response = send_request(api, "templates"; verb="POST")
end

"""
Returns an array with all experiments.

## Keyword Arguments
"""
function get_all_experiments(api::BaseAPI=BASE_API; limit=15, offset=0)
	params = Dict(
		"limit" => limit,
		"offset" => offset,
	) |> HTTP.Form

	response = send_request(api, "experiments/"; body=params)
	parse_response(response, Vector{Experiment})
end

"""
Return the `Experiment` with id `id`.
"""
function get_experiment(id::T, api::BaseAPI=BASE_API) where T <: Integer
	response = send_request(api, "experiments/$(id)")
	parse_response(response, Experiment)
end

"""
Return all items.

## Keyword Arguments
"""
function get_all_items(api::BaseAPI=BASE_API; limit=15, offset=0)
	params = Dict(
		:limit => limit,
		:offset => offset,
	)

	response = send_request(api, "items/"; query=params)
	parse_response(response, Vector{Item})
end

"""
Return the `Item` with id `id`.
"""
function get_item(id::T, api::BaseAPI=BASE_API) where T <: Integer
	response = send_request(api, "items/$(id)")
	parse_response(response, Experiment)
end

function get_item_types(api::BaseAPI=BASE_API)
	send_request(api, "item_types/")
end

function get_tags(api::BaseAPI=BASE_API)
	response = send_request(api, "tags/")
end

function get_upload(id::T, api::BaseAPI=BASE_API) where T <: Integer
	send_request(api, "uploads/$(id)")
end

function get_status(api::BaseAPI=BASE_API)
	send_request(api, "status/")
end

function get_all_templates(api::BaseAPI=BASE_API)
	response = send_request(api, "templates/")
end

function get_template(id::T, api::BaseAPI=BASE_API) where T <: Integer
	response = send_request(api, "templates/$(id)")
end

function post_experiment(x::Experiment, api::BaseAPI=BASE_API)
	body = struct_to_form_data(x)
	response = send_request(api, "experiments/$(x.id)"; verb="POST", body=body)
end

function post_item(x::Item, api::BaseAPI=BASE_API)
	body = struct_to_form_data(x)
	response = send_request(api, "items/$(x.id)"; verb="POST", body=body)
end

function post_template(x::Template, api::BaseAPI=BASE_API)
	body = struct_to_form_data(x)
	response = send_request(api, "templates/$(x.id)"; verb="POST", body=body)
end

function add_link_to_experiment(experiment_id::T, item_id::T, api::BaseAPI=BASE_API) where T <: Integer
	@show body = Dict("link" => item_id) |> HTTP.Form
	response = send_request(api, "templates/$(experiment_id)";
		verb="POST", body=body
	)
end

function add_link_to_experiment(x::Experiment, y::Item, api::BaseAPI=BASE_API)
	add_link_to_experiment(x.id, y.id, api)
end

function add_link_to_item(item_id::T, link_id::T, api::BaseAPI=BASE_API) where T <: Integer
	body = Dict("link" => link_id) |> HTTP.Form
	response = send_request(api, "items/$(item_id)"; verb="POST", body=body)
end

function add_link_to_item(x::Item, link::Item, api::BaseAPI=BASE_API) where T <: Integer
	add_link_to_item(x.id, link.id, api)
end

function add_link_to_experiment(x::Experiment, y::Item, api::BaseAPI=BASE_API)
	add_link_to_experiment(x.id, y.id, api)
end

function upload_to_experiment(x::Experiment, data::IO, api::BaseAPI=BASE_API)
	upload_to_experiment(x.id, data, api)
end

function upload_to_experiment(experiment_id::T, data::IO, api::BaseAPI=BASE_API) where T <: Integer
	response = send_request(api, "experiments/$(experiment_id)", data)
end

function upload_to_item(x::Item, data::IO, api::BaseAPI=BASE_API)
	upload_to_item(x.id, data, api)
end

function upload_to_item(item_id::T, data::IO, api::BaseAPI=BASE_API) where T <: Integer
	response = send_request(api, "items/$(item_id)", data)
end

function add_tag_to_experiment(x::Experiment, tag::Union{Symbol,String}, api::BaseAPI=BASE_API) where T
	add_tag_to_experiment(x.id, tag, api)
end

function add_tag_to_experiment(id::T, tag::Union{Symbol,String}, api::BaseAPI=BASE_API) where T
	body = Dict("tag" => tag) |> HTTP.Form
	send_request(api, "experiments/$(id)"; verb="POST", body=body)
end

function add_tag_to_item(x::Item, tag::Union{Symbol,String}, api::BaseAPI=BASE_API) where T
	add_tag_to_item(x.id, tag, api)
end

function add_tag_to_item(id::T, tag::Union{Symbol,String}, api::BaseAPI=BASE_API) where T
	body = Dict("tag" => tag) |> HTTP.Form
	send_request(api, "items/$(id)"; verb="POST", body=body)
end


text:
function _child_modules(m::Module)
    ns = names(m, imported = false, all = true)
    ms = []
    for n in ns
        if n != nameof(m)
            try
                x = Core.eval(m, n)
                if x isa Module
                    push!(ms, x)
                end
            catch
            end
        end
    end
    return ms
end

"""
    glasscatalogs()

Prints the complete list of glass catalogs available from GlassCat.

## Example
```julia-repl
julia> glasscatalogs()
41 glass catalogs:
GlassCat.AMTIR
GlassCat.ANGSTROMLINK
GlassCat.APEL
GlassCat.ARCHER
GlassCat.ARTON
GlassCat.AUER_LIGHTING
GlassCat.BIREFRINGENT
...
```
"""
function glasscatalogs()
    println("$(length(_child_modules(GlassCat))) glass catalogs:")
    for c in _child_modules(GlassCat)
        println("$c")
    end
end
export glasscatalogs

"""
    glasses(catalog::Union{Module,Nothing} = nothing)

Prints the glasses available from a given catalog, or if none if specified from all catalogs.

# Example
```julia-repl
julia> glasses(GlassCat.CARGILLE)
GlassCat.CARGILLE - 3 glasses:
	OG0607
	OG0608
	OG081160
```
"""
function glasses(catalog::Union{Module,Nothing} = nothing)
    function printforcat(cat::Module)
        glass_names = names(cat, all = true, imported = false)
        glasses = []
        for glass_name in glass_names
            glass_name_str = string(glass_name)
            if !occursin("#", glass_name_str) && glass_name_str != "eval" && glass_name_str != "include" && glass_name != nameof(cat)
                push!(glasses, "$glass_name_str")
            end
        end
        println("$cat - $(length(glasses)) glasses:")
        for g in glasses
            println("\t$g")
        end
    end
    if catalog === nothing
        for c in _child_modules(GlassCat)
            printforcat(c)
        end
    else
        printforcat(catalog)
    end
end
export glasses

"""
    findglass(condition::Function) -> Vector{Glass}

Returns the list of glasses which satisfy `condition` where `condition::(Glass -> Bool)`.

TODO - make the condition easier to specify (accessor functions for fields?)

# Example
```julia-repl
julia> findglass(x -> (x.Nd > 2.3 && x.λmin < 0.5 && x.λmax > 0.9))
8-element Array{GlassCat.Glass,1}:
 BIREFRINGENT.TEO2_E
 BIREFRINGENT.PBMOO4
 BIREFRINGENT.LINBO3
 INFRARED.CLEARTRAN_OLD
 INFRARED.CLEARTRAN
 INFRARED.SRTIO3
 INFRARED.ZNS_BROAD
 INFRARED.ZNS_VIS
```
"""
function findglass(condition::Function)
    out = Vector{Glass}(undef, 0)
    for g in AGF_GLASSES
        if condition(g)
            push!(out, g)
        end
    end
    for g in OTHER_GLASSES
        if condition(g)
            push!(out, g)
        end
    end
    return out
end


text:
using GoogleMaps
using Base.Test

# write your own tests here
@testset "Decode Polyline Algorithm" begin
    lonlats = [(-120.2, 38.5), (-120.95, 40.7), (-126.453, 43.252)]
    for (i,lonlat) in enumerate(GoogleMaps.decodepolyline("_p~iF~ps|U_ulLnnqC_mqNvxq`@"))
        @test all(lonlat .≈ lonlats[i])
    end
end

@testset "Geocoding" begin
    result = GoogleMaps.geocode("1600 Amphitheatre Parkway, Mountain View, CA")
    @test !isempty(result)
end

@testset "Routing" begin
    sleep(10)
    result1 = GoogleMaps.directions(
        "MBTA Kendall Square",
        "Logan International Airport"
    )
    sleep(10)
    result2 = GoogleMaps.directions(
        (-71.0862274, 42.3622273), # Kendall Square
        (-71.0095602, 42.3656132) # Logan Airport
    )
    @test !isempty(result1)
    @test !isempty(result2)
    @test length(result1) == length(result2)
end

# results3 = GoogleMaps.details("AIzaSyDuGVOl_C0G8vGT8buABSZcYNl8X7Dhho4", "ChIJN1t_tDeuEmsRUsoyG83frY4")

# result = GoogleMaps.matchroads("AIzaSyDuGVOl_C0G8vGT8buABSZcYNl8X7Dhho4", [(149.12958,-35.27801), (149.12907,-35.28032), (149.12929,-35.28099), (149.12984,-35.28144), (149.13003,-35.28194), (149.12956,-35.28282), (149.12881,-35.28302), (149.12836,-35.28473)], interpolate=true, snap=true)
# result1 = GoogleMaps.matchroads("AIzaSyDuGVOl_C0G8vGT8buABSZcYNl8X7Dhho4", [(24.942795,60.170880), (24.942796,60.170879), (24.942796, 60.170877)])

# results2 = GoogleMaps.nearby("AIzaSyDuGVOl_C0G8vGT8buABSZcYNl8X7Dhho4", (151.1957362,-33.8670522),500, keyword="cruise")

# results4 = GoogleMaps.searchfor("AIzaSyDuGVOl_C0G8vGT8buABSZcYNl8X7Dhho4", "restaurants in Sydney")

# results5 = GoogleMaps.photo("AIzaSyDuGVOl_C0G8vGT8buABSZcYNl8X7Dhho4", "CnRtAAAATLZNl354RwP_9UKbQ_5Psy40texXePv4oAlgP4qNEkdIrkyse7rPXYGd9D_Uj1rVsQdWT4oRz4QrYAJNpFX7rzqqMlZw2h2E2y5IKMUZ7ouD_SlcHxYq1yL4KbKUv3qtWgTK0A6QbGh87GB3sscrHRIQiG2RrmU_jF4tENr9wGS_YxoUSSDrYjWmrNfeEHSGSc3FyhNLlBU")


text:
__precompile__()

module Atom

using Juno, Lazy, JSON, MacroTools, Reexport, Media, Base.StackTraces

import Requires
import Media: @dynamic

function __init__()
  Juno.activate()

  atreplinit(i -> fixdisplayorder())

  Atom.handle("connected") do
    if isREPL()
      reset_repl_history()
      fixdisplayorder()
    end
    nothing
  end

  Requires.@require WebIO="0f1e0344-ec1d-5b48-a673-e5cf874b6c29" begin
    include("display/webio.jl")
  end
end

include("comm.jl")
include("display/display.jl")
include("progress.jl")
include("eval.jl")
include("workspace.jl")
include("repl.jl")
include("docs.jl")
include("completions.jl")
include("misc.jl")
include("frontend.jl")
include("utils.jl")

include("debugger/load.jl")

include("profiler/profiler.jl")

end # module


text:
const BOOL_VALS = Union{Val{:(&&)}, Val{:(||)}, Val{:(⊻)}}

bool_val_to_set(::Val{:(&&)}) = AndSet
bool_val_to_set(::Val{:(||)}) = OrSet
bool_val_to_set(::Val{:(⊻)}) = XorSet

function _build_bool_constraint(
    _error::Function,
    lhs,
    rhs,
    set_type
)
    lhs_set = JuMP.moi_set(lhs)
    rhs_set = JuMP.moi_set(rhs)

    lhs_jump_func = JuMP.jump_function(lhs)
    rhs_jump_func = JuMP.jump_function(rhs)

    lhs_func = JuMP.moi_function(lhs)
    rhs_func = JuMP.moi_function(rhs)

    func = [lhs_jump_func..., rhs_jump_func...]
    return JuMP.VectorConstraint(
        func, set_type{typeof(lhs_func), typeof(rhs_func)}(lhs_set, rhs_set)
    )
end

"""
    transform_binary_expr(sym::Symbol)

Transform a symbol to a constraint of the form Symbol == 1
"""
function transform_binary_expr(sym::Symbol)
    return :($sym == 1)
end

"""
    transform_binary_expr(expr::Expr)

Transform a ! (symbol) to a constraint of the form Symbol == 0
or x[...] to x[...] == 1
"""
function transform_binary_expr(expr::Expr)
    if expr.head == :ref
        expr = :($expr == 1)
    elseif expr.head == :call && expr.args[1] == :! && (expr.args[2] isa Symbol || expr.args[2].head == :ref)
        expr = :($(expr.args[2]) == 0)
    end
    return expr
end

function parse_bool_constraint(_error, bool_val::BOOL_VALS, lhs, rhs)
    _error1 = deepcopy(_error)
    # allow a || b instead of a == 1 || b == 1
    lhs = transform_binary_expr(lhs)

    lhs_vectorized, lhs_parsecode, lhs_buildcall =
        JuMP.parse_constraint_expr(_error, lhs)

    if lhs_vectorized
        _error("`$(lhs)` should be non vectorized. There is currently no vectorized support for `and` constraints. Please open an issue at ConstraintSolver.jl")
    end

    rhs = transform_binary_expr(rhs)
    rhs_vectorized, rhs_parsecode, rhs_buildcall =
        JuMP.parse_constraint_expr(_error1, rhs)

    if rhs_vectorized
        _error("`$(rhs)` should be non vectorized. There is currently no vectorized support for `and` constraints. Please open an issue at ConstraintSolver.jl")
    end

    # TODO implement vectorized version
    vectorized = false
    complete_parsecode = quote
        $lhs_parsecode
        $rhs_parsecode
    end

    bool_set = bool_val_to_set(bool_val)

    buildcall = :($(esc(:(CS._build_bool_constraint)))(
        $_error,
        $lhs_buildcall,
        $rhs_buildcall,
        $bool_set
    ))
    return vectorized, complete_parsecode, buildcall
end

function JuMP.parse_constraint_head(_error::Function, bool_val::BOOL_VALS, lhs, rhs)
    return parse_bool_constraint(_error, bool_val, lhs, rhs)
end

function JuMP.parse_one_operator_constraint(_error::Function, vectorized::Bool, bool_val::BOOL_VALS, lhs, rhs)
    @assert !vectorized
    v,c,b = parse_bool_constraint(_error, bool_val, lhs, rhs)
    return c,b
end


text:
include("silicon_runners.jl")
using DoubleFloats

@testset "SCF of silicon without exchange-correlation (small)" begin
    run_silicon_redHF(Float64, Ecut=5, test_tol=0.05, n_ignored=0, grid_size=15)
end

if !isdefined(Main, :FAST_TESTS) || !FAST_TESTS
    @testset "SCF of silicon without exchange-correlation (large)" begin
        run_silicon_redHF(Float64, Ecut=25, test_tol=5e-6, n_ignored=2, grid_size=35,
                          scf_tol=1e-7)
    end
end

@testset "Silicon LDA (small, Double32)" begin
    run_silicon_redHF(Double32, Ecut=5, test_tol=0.05, n_ignored=0, grid_size=15)
end


text:
# sets can only contain unique elements - adding an additional
# 1 will appear to do 'nothing'.
s1 = Set([1,2,3,"Richard", 3])
println(typeof(s1))
push!(s1, 1,55,98)
println(s1)
println("******")

# if I create a set of integers I cannot
# add a string
s2 = Set([1,2,3])
println(typeof(s2))

# to explicitly create a Set of a particular type
s3 = Set{Float64}([1,2,3,4,5,66])
println(typeof(s3))
println(s3)

# I can check if a Set contains a value
println(1 in s3)
println(34 in s3)

# I can iterate over a Set - but order is not guarenteed
# and I can't  use an index to retrieve a element from the Set
for e in s3
    println("$(e * 5)")
end

# with Sets it is possible to calculate the interset (common elements) of two Sets
s4 = Set([1,4,6,7])
s5 = Set([4,7,8,9])
println(intersect(s4, s5))

# to output the 'union' (all unique elements from both sets) of two sets
println(union(s4, s5))

# to calculate the 'elements' that are in set1 but not in set2
println(setdiff(s5, s4))

max_stars_repo_name: rpillar/JuliaExamples

max_issues_repo_name: rpillar/JuliaExamples

max_forks_repo_name: rpillar/JuliaExamples


text:
# this script produces Figure 4

using CSV, DataFrames, DelimitedFiles, Glob, GMT, Statistics
stats = expanduser("~/SeisNoise_paper/DATA/gmap-stations.txt")
df = DataFrame!(CSV.File(stats,delim="|"))
lat,lon = df[:,:Latitude], df[:,:Longitude]
latmin, latmax = minimum(lat), maximum(lat)
lonmin, lonmax = minimum(lon), maximum(lon)
latpad = 0.025
lonpad = 0.05

# create figure
coast(
    region=(lonmin-lonpad,lonmax+lonpad,latmin-latpad,latmax+latpad),
    proj=:Mercator,
    res=:full,
    land=:white,
    water=:green,
    # scale=0.03,
    par=(:MAP_FRAME_TYPE,"fancy+"),
    map_scale="jBR+c$(mean(lat))+w10k+f+o2.5/2.5+at+u",
    figsize=15,
)
scatter!(lon,lat,fill=:gold,marker=:triangle,markeredgecolor=:black,markersize=0.2)
basemap!(inset=(anchor=:TL, width=0.25, offset=(1, 5), save="xx000"))
t = readdlm("xx000")
coast!(region=[-107 -91 26 41], proj=:merc,
           land=:lightgray, area=5000, shore=:faint,
           x_off=t[1], y_off=t[2],
           N="a",
           figsize=4,
           frame=:bare,
)
scatter!(
     [mean(lon)],
     [mean(lat)],
     marker=:star,
     markeredgecolor=:black,
     fill=:gold,
     markersize=0.5,
     savefig=expanduser("~/SeisNoise_paper/LASSOmap.png"),
)


text:
using ArrayAllocators
using SafeByteCalculators
using NumaAllocators
using Documenter

DocMeta.setdocmeta!(ArrayAllocators, :DocTestSetup, :(using ArrayAllocators); recursive=true)

makedocs(;
    modules=[ArrayAllocators, SafeByteCalculators, NumaAllocators],
    authors="Mark Kittisopikul <kittisopikulm@janelia.hhmi.org> and contributors",
    repo="https://github.com/mkitti/ArrayAllocators.jl/blob/{commit}{path}#{line}",
    sitename="ArrayAllocators.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://mkitti.github.io/ArrayAllocators.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Detailed Overview" => "detailed_overview.md",
        "Byte Calculators" => "bytecalculators.md",
        "NUMA Allocators" => "numa.md",
        "Compositions" => "compositions.md",
        "API" => [
            "ArrayAllocators.jl" => "api.md",
            "ArrayAllocators.ByteCalculators" => "bytecalculators_api.md",
            "NumaAllocators.jl" => "numa_api.md",
        ]
    ],
)

deploydocs(;
    repo="github.com/mkitti/ArrayAllocators.jl",
    devbranch="main",
    push_preview = true
)


text:
module MPSKitModels
    using TensorKit,MPSKit
    using LinearAlgebra:Diagonal,diag
    using MPSKit:@plansor;

    export spinmatrices,nonsym_spintensors,nonsym_bosonictensors
    include("utility.jl")

    export nonsym_ising_ham,nonsym_ising_mpo,z2_ising_mpo
    include("ising.jl");

    export nonsym_xxz_ham,su2_xxx_ham,u1_xxz_ham
    include("xxz.jl");

    export su2u1_grossneveu,su2u1_orderpars
    include("grossneveu.jl");

    export nonsym_qstateclock_mpo
    include("qstateclock.jl");

    export nonsym_qed_qlm_ham,qed_qlm_G2,u1_qed_ham
    include("qed_qlm.jl")

    export nonsym_sixvertex_mpo,u1_sixvertex_mpo,cu1_sixvertex_mpo
    include("sixvertex.jl")

    export U1_strip_harper_hofstadter
    include("hofstadter.jl")
end

max_stars_repo_name: maartenvd/MPSKitModels

max_issues_repo_name: maartenvd/MPSKitModels

max_forks_repo_name: maartenvd/MPSKitModels


text:
function LeastSquaresAdapter(arg0::LeastSquaresProblem)
    return LeastSquaresAdapter((LeastSquaresProblem,), arg0)
end

function evaluate(obj::LeastSquaresAdapter, arg0::RealVector)
    return jcall(obj, "evaluate", LeastSquaresProblem_Evaluation, (RealVector,), arg0)
end

function get_convergence_checker(obj::LeastSquaresAdapter)
    return jcall(obj, "getConvergenceChecker", ConvergenceChecker, ())
end

function get_evaluation_counter(obj::LeastSquaresAdapter)
    return jcall(obj, "getEvaluationCounter", Incrementor, ())
end

function get_iteration_counter(obj::LeastSquaresAdapter)
    return jcall(obj, "getIterationCounter", Incrementor, ())
end

function get_observation_size(obj::LeastSquaresAdapter)
    return jcall(obj, "getObservationSize", jint, ())
end

function get_parameter_size(obj::LeastSquaresAdapter)
    return jcall(obj, "getParameterSize", jint, ())
end

function get_start(obj::LeastSquaresAdapter)
    return jcall(obj, "getStart", RealVector, ())
end


text:
#=networking.jl - Auth: Kino Roy for Shogi, an assignment for CMPT276:
 provides a server on port 8080 where up to two client
may connect to play Shogi.
Current state: accepts starting payload to set game settings, allows users to make
moves with wincode 2, sends move to opponent =#
module networking

function getAt(arr, index,backup)
  local res
  try
    res = arr[index]
  catch err
    res = backup
  end
  return res
end

function serverStart()
  global server = listen(8080)
  global settingsSet = false
  global output = ""
  global player1
  global player2
  global timeLeftPlayer1
  global timeLeftPlayer2
  global whosTurn
  global timeElapsed = Dict()
  global limitAdd
  global timeLimit = 1
#  global wincode,sL
  global connections = Dict() #Array(TCPSocket,0)
  listenForClientsTask = Task(listenForClients)
  schedule(listenForClientsTask)
end

function serverClose()
  global server
  for i in values(connections)
    close(i)
  end
  close(server)
end

function listenForClients()
  while true
    local client
    try
      client = accept(server)
    catch
      println("goodbye!")
      close(client)
      break;
    end
    (host,port)=getsockname(client)
    connections[port] = client
    timeElapsed[port] = 0
    clientCommTask = @task(handleClientComm(client))
    schedule(clientCommTask)
  end
end

function handleClientComm(client)
  (host,port) = getsockname(client)
  message = Array(UInt8,4096)
   while (true)
  global timeLimit
      if timeElapsed[port]>timeLimit #player went over time limit, player resigns
        opponent = client == player1 ? player2 : player1
        write(opponent,"""9:0:3:0:0:0:0:0:0:0:0\n""") #send resign move to opponent
        write(client,"""e:you ran out of time and were forced to resign\n""")
        prinln("resigned due to timeout")
      end
       bytesRead = 1;
       local wincode,sL
       try
           #blocks until a client sends a message
           if settingsSet
             tic()
           end
           message = (readline(client))
           println(message)
           if (isempty(message))
             println("disconnect")
             close(client)
             delete!(connections,port)
             for i in values(connections)
               try
                 write(i,"\n3\n")
               end
             end
             break;
           end
           message = chomp(message)
           sL = split(message,r":",keep=false)
           wincode = getAt(sL,1,-1)
        #=  (wincode,authString,movenum,movetype,sourcex,sourcey,targetx,targety,option,cheating,
targetx2,targety2) = (getAt(sL,1,-1),getAt(sL,2,-1),getAt(sL,3,-1),getAt(sL,1,-1),getAt(sL,1,-1),getAt(sL,1,-1),
getAt(sL,1,-1),getAt(sL,1,-1),getAt(sL,1,-1),getAt(sL,1,-1))=# #possibly iterate through these size of sL and set or null??
          # write(client,message)
       catch err
         #println("socket error: $err")
           #a socket error has occured
        #   write(client,"e:Readline failed with $err")
           break;
       end
       if (isempty(message)) #Doesnt detect

           #the client has disconnected from the server
           #println("Player on port $port disconected :(")
           #write(client,"e:sorry, player disconected")
          # break;
        end
        # message recieved
        global settingsSet
        if wincode == "0" && !settingsSet

          #"<wincode>: <gametype>: <legality>: <timeLimit>: <limitAdd>"
          global timeLimit
          global limitAdd
          (gametype,legality,timeLimit,limitAdd) = (sL[2],sL[3],sL[4],sL[5])
          timeLimit = parse(timeLimit)
          limitAdd = parse(limitAdd)
          global output = "$gametype:$legality:$timeLimit:$limitAdd"
          #prints to player 1
          #=try
            write(connections[player1],"0:$(player1*10^3):$gametype:$legality:$timeLimit:$limitAdd")
          catch err
            println("$err")
          end =#
          #prints to player 2

        elseif wincode == "1" #Quit the game
          #Draw

        elseif wincode == "2" #Play a move
          timeElapsedCurrentMove = toc()
          (wincode,authString,movenum,movetype,sourcex,sourcey,targetx,targety,option,cheating,
  targetx2,targety2) = (getAt(sL,1,-1),getAt(sL,2,-1),getAt(sL,3,-1),getAt(sL,4,-1),getAt(sL,5,-1),getAt(sL,6,-1),
  getAt(sL,7,-1),getAt(sL,8,-1),getAt(sL,9,-1),getAt(sL,10,-1),getAt(sL,11,-1),getAt(sL,12,-1))

          timeElapsed[port]+=(timeElapsedCurrentMove-limitAdd)
          if whosTurn == port #check if its the clients turn
            opponent = client == player1 ? player2 : player1
            write(opponent,"\n9:$movenum:$movetype:$sourcex:$sourcey:$targetx:$targety:$option:$cheating
$targetx2:$targety2\n") #send move to opponent
          else
          #  println("$whosTurn:$port")
            write(client,"\n8\n")
          end
          (host1,port1) = getsockname(player1)
          (host2,port2) = getsockname(player2)
          global whosTurn = port1 == port ? port2 : port1  #After move is made, change the turn
        elseif wincode == "3" #Accuse opponent of cheating

        elseif wincode == "10" #Bad payload

        end
        if !settingsSet
          global settingsSet = true
            println("settings set")
          bothPlayersConnected = false
          while !(size(collect(keys(connections)))[1]==2) #waits until both players connected
            sleep(.5)
          end
          bothPlayersConnected = size(collect(keys(connections)))[1]==2

          if bothPlayersConnected
            println("bothPlayersConnected")
            #Roll to determine player 1 and player 2
            player1index = rand([1,2])
            player2index = player1index == 1 ? 2 : 1
            global player1 = collect(values(connections))[player1index]
            global player2 = collect(values(connections))[player2index]
            (host1,port1) = getsockname(player1)
            global whosTurn = port1
          #  println("$port1:first turn")
            try
              (host1,port1) = getsockname(player1)
              (host2,port2) = getsockname(player2)
                write(player1,"""$(player1index-1):$(port1*10^3):$gametype:$legality:$timeLimit:$limitAdd\n""")
                write(player2,"""$(player2index-1):$(port2*10^3):$gametype:$legality:$timeLimit:$limitAdd\n""")
            catch err
              println(err)
            end
          end

        end
        #Prints output to all connections

    #= for i in values(connections)
          try
            global output
            write(i,output)
          catch err #Someone disconected
            disconects+=1
            println(err)
            continue
          end
        end=# #need to check if socket is closed, then print to opponent wincode 3
      end
    end

end #end module

max_stars_repo_name: kinoroy/Shogi

max_issues_repo_name: kinoroy/Shulia

max_forks_repo_name: kinoroy/Shulia


text:
module BayesianQuadrature

using AbstractGPs
using AbstractMCMC
using Distributions
using KernelFunctions
using LinearAlgebra
using PDMats
using Random

export BayesQuad, LogBayesQuad
export PriorSampling
export BayesModel
export prior, integrand, logintegrand, logprior, logjoint
export BQ # Short version for calling BayesianQuadrature
export KernelEmbedding
export kernel_mean, kernel_variance

const BQ = BayesianQuadrature

"""
    AbstractBQ

General class of models to perform quadrature.
It needs to implement `quadrature(::AbstractBQ, model::AbstractBQModel, samples)`
"""
abstract type AbstractBQ{TK} end

"""
    AbstractBQSampler

General class for sampling.
Should implement the `AbstracMCMC` API
"""
abstract type AbstractBQSampler <: AbstractMCMC.AbstractSampler end

"""
    AbstractBQModel

General model class.
Should implement `p_0(m)` and `logintegrand(m)`
"""
abstract type AbstractBQModel{Tp,Ti} <: AbstractMCMC.AbstractModel end

include("bayesquads/abstractbq.jl")
include("samplers/abstractbqsampler.jl")
include("kernelembeddings/kernelembedding.jl")
include("interface.jl")
include("models.jl")
include("utils.jl")

end


text:
"""
    reorder(A::AbstractDimArray, order::Pair) => AbstractDimArray
    reorder(A::Dimension, order::Order) => AbstractDimArray

Reorder every dims index/array to `order`, or reorder index for
the the given dimension(s) to the `Order` they wrap.

`order` can be an [`Order`](@ref), or `Dimeension => Order` pairs.

If no axis reversal is required the same objects will be returned, without allocation.
"""
function reorder end

reorder(x, p::Pair, ps::Vararg{<:Pair}) = reorder(x, (p, ps...))
reorder(x, ps::Tuple{Vararg{<:Pair}}) = reorder(x, Dimensions.pairdims(ps...))
# Reorder specific dims.
reorder(x, dimwrappers::Tuple) = _reorder(x, dimwrappers)
# Reorder all dims.
reorder(x, ot::Order) = reorder(x, typeof(ot))
reorder(x, ot::Type{<:Order}) = _reorder(x, map(d -> rebuild(d, ot), dims(x)))
reorder(dim::Dimension, ot::Type{<:Order}) =
    ot <: basetypeof(order(dim)) ? dim : reverse(dim)

# Recursive reordering. x may be reversed here
_reorder(x, orderdims::DimTuple) = _reorder(reorder(x, orderdims[1]), tail(orderdims))
_reorder(x, orderdims::Tuple{}) = x

reorder(x, orderdim::Dimension) = _reorder(val(orderdim), x, dims(x, orderdim))

_reorder(neworder::Order, x, dim::DimOrDimType) = _reorder(basetypeof(neworder), x, dim)
# Reverse the dimension index
_reorder(::Type{O}, x, dim::DimOrDimType) where O<:Ordered =
    order(dim) isa O ? x : reverse(x; dims=dim)
_reorder(ot::Type{Unordered}, x, dim::DimOrDimType) = x

"""
    modify(f, A::AbstractDimArray) => AbstractDimArray
    modify(f, s::AbstractDimStack) => AbstractDimStack
    modify(f, dim::Dimension) => Dimension
    modify(f, x, lookupdim::Dimension) => typeof(x)

Modify the parent data, rebuilding the object wrapper without
change. `f` must return a `AbstractArray` of the same size as the original.

This method is mostly useful as a way of swapping the parent array type of
an object.

## Example

If we have a previously-defined `DimArray`, we can copy it to an Nvidia GPU with:

```julia
A = DimArray(rand(100, 100), (X, Y))
modify(CuArray, A)
```

This also works for all the data layers in a `DimStack`.
"""
function modify end
modify(f, s::AbstractDimStack) = map(a -> modify(f, a), s)
function modify(f, A::AbstractDimArray)
    newdata = f(parent(A))
    size(newdata) == size(A) || error("$f returns an array with a different size")
    rebuild(A, newdata)
end
modify(f, x, dim::DimOrDimType) = set(x, modify(f, dims(x, dim)))
modify(f, dim::Dimension) = rebuild(dim, modify(f, val(dim)))
function modify(f, lookup::LookupArray)
    newindex = modify(f, parent(lookup))
    rebuild(lookup; data=newindex)
end
function modify(f, index::AbstractArray)
    newindex = f(index)
    size(newindex) == size(index) || error("$f returns a vector with a different size")
    newindex
end

"""
    brodadcast_dims(f, sources::AbstractDimArray...) => AbstractDimArray

Broadcast function `f` over the `AbstractDimArray`s in `sources`, permuting and reshaping
dimensions to match where required. The result will contain all the dimensions in
all passed in arrays in the order in which they are found.

## Arguments

- `sources`: `AbstractDimArrays` to broadcast over with `f`.

This is like broadcasting over every slice of `A` if it is
sliced by the dimensions of `B`.
"""
function broadcast_dims(f, As::AbstractDimArray...)
    dims = combinedims(As...)
    T = Base.Broadcast.combine_eltypes(f, As)
    broadcast_dims!(f, similar(first(As), T, dims), As...)
end

"""
    brodadcast_dims!(f, dest::AbstractDimArray, sources::AbstractDimArray...) => dest

Broadcast function `f` over the `AbstractDimArray`s in `sources`, writing to `dest`.
`sources` are permuting and reshaping dimensions to match where required.

The result will contain all the dimensions in all passed in arrays, in the order in
which they are found.

## Arguments

- `dest`: `AbstractDimArray` to update.
- `sources`: `AbstractDimArrays` to broadcast over with `f`.
"""
function broadcast_dims!(f, dest::AbstractDimArray{<:Any,N}, As::AbstractDimArray...) where {N}
    As = map(As) do A
        isempty(otherdims(A, dims(dest))) || throw(DimensionMismatch("Cannot broadcast over dimensions not in the dest array"))
        # comparedims(dest, dims(A, dims(dest)))
        # Lazily permute B dims to match the order in A, if required
        if !dimsmatch(commondims(A, dest), commondims(dest, A))
            PermutedDimsArray(A, commondims(dest, A))
        else
            A
        end
    end
    od = map(A -> otherdims(dest, dims(A)), As)

    # Broadcast over b for each combination of dimensional indices D
    if all(map(isempty, od))
        dest .= f.(As...)
    else
        not_shared_dims = combinedims(od...)
        reshaped = map(As) do A
            all(hasdim(A, dims(dest))) ? parent(A) : _insert_length_one_dims(A, dims(dest))
        end
        dest .= f.(reshaped...)
    end
    return dest
end

function _insert_length_one_dims(A, alldims)
    lengths = map(alldims) do d
        hasdim(A, d) ? size(A, d) : 1
    end
    return reshape(parent(A), lengths)
end

@deprecate dimwise broadcast_dims
@deprecate dimwise! broadcast_dims!

# Get a tuple of unique keys for DimArrays. If they have the same
# name we call them layerI.
function uniquekeys(das::Tuple{AbstractDimArray,Vararg{<:AbstractDimArray}})
    uniquekeys(Symbol.(map(name, das)))
end
function uniquekeys(keys::Tuple{Symbol,Vararg{<:Symbol}})
    ids = ntuple(x -> x, length(keys))
    map(keys, ids) do k, id
        count(k1 -> k == k1, keys) > 1 ? Symbol(:layer, id) : k
    end
end


text:
Base.@kwdef mutable struct Resources <: AMGXObject
    handle::API.AMGX_resources_handle = API.AMGX_resources_handle(C_NULL)
    ref_count::Threads.Atomic{Int} = Threads.Atomic{Int}(0)
    cfg::Union{Config, Nothing} = nothing
    function Resources(handle::API.AMGX_resources_handle, ref_count::Threads.Atomic{Int},
                       cfg::Union{Config, Nothing})
        resources = new(handle, ref_count, cfg)
        finalizer(warn_not_destroyed_on_finalize, resources)
        return resources
    end
end
get_api_destroy_call(::Type{Resources}) = API.AMGX_resources_destroy
function dec_refcount_parents(resources::Resources)
    dec_refcount!(resources.cfg)
    resources.cfg = nothing
    nothing
end

function create!(resources::Resources, cfg::Config)
    res_handle_ptr = Ref{API.AMGX_resources_handle}()
    @checked API.AMGX_resources_create_simple(res_handle_ptr, cfg.handle)
    resources.handle = res_handle_ptr[]
    resources.cfg = cfg
    inc_refcount!(resources.cfg)
    return resources
end

Resources(cfg::Config) = create!(Resources(), cfg)


text:
export HydraulicPositioningSystem

"""
Control of a carriage along a track using a hydraulic positioning system.

The plant model for this problem is defined by the following parameters:
* `B` - Compressibility modulus
* `m` - Mass
* `R` - Friction coefficient
* `Rc` - Coulomb friction coefficient
* `KL` - Bypass oil coefficient
* `σ` - Hydraulic constant
* `V` - Volume
* `A` - Piston area

From E. J. Davison, ‘Benchmark Problems for Control System Design: Report of the IFAC Theory Committee’, IFAC, 1990.
"""
mutable struct HydraulicPositioningSystem <: ControlBenchmarkProblem
    B::Float64      # Compressibility modulus
    m::Float64      # Mass
    R::Float64      # Friction coefficient
    Rc::Float64     # Coulomb friction coefficient
    KL::Float64     # Bypass oil coefficient
    σ::Float64      # Hydraulic constant
    V::Float64      # Volume
    A::Float64      # Piston area
end

function HydraulicPositioningSystem()
    HydraulicPositioningSystem(
        14000.0,   # Nominal compressibility modulus
        0.1287,    # Nominal mass
        0.150,     # Nominal friction coefficient
        0.01,      # Nominal Coulomb friction coefficient
        0.002,     # Nominal bypass oil coefficient
        0.24,      # Nominal hydraulic constant
        874.0,     # Nominal volume
        10.75 )    # Nominal piston area
end

@addbenchmark HydraulicPositioningSystem "ifac07" "Hydraulic positioning system" "Hydraulically controlled carriage on a track" Set([:linear, :stable, :continuous, :mechanical, :ifac])

"""
    parameterbounds(::HydraulicPositioningSystem)
"""
function parameterbounds( ::HydraulicPositioningSystem )
    # The intervals the parameters should be inside
    return Dict( "B"  => Bound( 9000, 16000 ),
                 "m"  => Bound( 0.05, 0.3 ),
                 "R"  => Bound( 0.05, 5.0 ),
                 "Rc" => Bound( 0.0, 0.05 ),
                 "KL" => Bound( 0.000103, 0.0035 ),
                 "σ"  => Bound( 0.001, 15.0 ),
                 "V"  => Bound( 0.0, Inf ),
                 "A"  => Bound( 10.5, 11.10 ) )
end

"""
    controlbenchmark(p::HydraulicPositioningSystem; ignoreBounds::Bool = false)

```@example
hps = HydraulicPositioningSystem()

# Change the mass from its default
hps.m = 0.25

sys, inputs, disturbance = controlbenchmark( hps )
```
"""
function controlbenchmark( p::HydraulicPositioningSystem; ignoreBounds::Bool = false )

    if( !ignoreBounds )
        @checkbounds HydraulicPositioningSystem p
    end

    # The state transition matrix
    A = fill( 0.0, (3, 3) )
    A[1, 2] = 1.0
    A[2, 2] = -( p.R + 4.0 / π * p.Rc ) / p.m
    A[3, 2] = -( 4.0 * p.B ) / p.V * p.A
    A[2, 3] = p.A / p.m
    A[3, 3] = -( 4.0 * p.B ) / p.V * ( p.σ + p.KL )

    # The input matrix
    B = fill( 0.0, (3, 1) )
    B[3, 1] = -( 4.0 * p.B ) / p.V;

    # The output matrix
    C = fill( 0.0, (1, 3) )
    C[1, 1] = 1.0

    # The feed-forward matrix
    D = fill( 0.0, (1, 1) )

    # The disturbance matrix
    # TODO: Extend ss to take this disturbance matrix
    E = fill( 0.0, (3, 1) )
    E[2, 1] = -1.0 / p.m

    # Input constraints
    inputBound = Bound( -866.0, 866.0 )

    # Disturbance bound for FL
    disturbanceBound = Bound( 0.0, 800.0 )

    # Create the constrained system
    return (; :sys => StateSpace( A, B, C, D ),
              :inputBounds => inputBound,
              :disturbanceBounds => disturbanceBound )
end


text:
using Luna

a = 15e-6
gas = :Ar
pres = 25

τfwhm = 600e-15
λ0 = 800e-9
flength = 80e-2
energy = 10e-6

grid = Grid.EnvGrid(flength, λ0, (220e-9, 3000e-9), 4e-12)

m = Capillary.MarcatiliMode(a, gas, pres, loss=false)
aeff = let m=m
    z -> Modes.Aeff(m, z=z)
end

energyfun, energyfunω = Fields.energyfuncs(grid)

densityfun = let dens0=PhysData.density(gas, pres)
    z -> dens0
end

linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0)

ionpot = PhysData.ionisation_potential(gas)
ionrate = Ionisation.ionrate_fun!_ADK(ionpot)

responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),)
            # Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot))

inputs = (Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy), Fields.ShotNoise())

Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)

statsfun = Stats.default(grid, Eω, m, linop, transform; gas=gas, windows=((150e-9, 300e-9),))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
Luna.run(Eω, grid, linop, transform, FT, output)

##
Plotting.pygui(true)
Plotting.stats(output)
Plotting.prop_2D(output; trange=(-1e-12, 1e-12), λrange=(220e-9, 2000e-9))
Plotting.time_1D(output; trange=(-1e-12, 1e-12))
Plotting.spec_1D(output)


text:
const FORM_TEXT_MIMES = (MIME"text/plain", MIME"text/latex")

form_chars_for_mime(::MIME"text/plain") = ("{", "|", "}")
form_chars_for_mime(::MIME"text/latex") = (raw"\left\{ ", raw" \middle| ", raw" \right\}")

"""
Show the canonical form of the surreal number x.
"""
function show_form(io::IO, mime::Union{FORM_TEXT_MIMES...}, x)
	l, m, r = form_chars_for_mime(mime)
    print(io, l)
    show_set(io, mime, left(x))
    print(io, m)
    show_set(io, mime, right(x))
    print(io, r)
end

show_form_outer(io::IO, mime::MIME, x) = show_form(io, mime, x)

function show_form_outer(io::IO, mime::MIME"text/latex", x)
	print(io, "\$")
	show_form(io, mime, x)
	print(io, "\$")
end

"""
Show a set of surreal numbers within the representation of a form.
"""
function show_set(io::IO, mime::Union{FORM_TEXT_MIMES...}, s)
    first = true
    for x in s
        first || print(io, ", ")
        show_inner(io, mime, x)
        first = false
    end
end

"""
Show the surreal number `x` inside the left or right subset of a form.
"""
show_inner(io::IO, mime::MIME, x) = show(io, mime, x)


text:
function prob(beta::Real, h::Real)
	return 1.0/(1.0+exp(-beta*h))    # beta = 1/T
end

function stochastic_sign{T<:Real}(p::T)
    return p > rand() ? one(T) : -one(T)
end

function activatefunc(beta::Float64, h::Float64)
    return stochastic_sign(prob(beta,h))
end

function activatefunc(beta::Float64, h::Complex{Float64})
    r_val = stochastic_sign(prob(beta,real(h)))
	i_val = stochastic_sign(prob(beta,imag(h)))
	return complex(r_val, i_val)
end

function activatefunc(beta::Float64, h::Quaternion{Float64})
    r_val = stochastic_sign(prob(beta,real(h)))
	i_val = stochastic_sign(prob(beta,imagi(h)))
	j_val = stochastic_sign(prob(beta,imagj(h)))
	k_val = stochastic_sign(prob(beta,imagk(h)))
	return quaternion(r_val, i_val, j_val, k_val)
end

function internal_potential(W::AbstractArray{Quaternion{Float64}}, x::AbstractVector{Quaternion{Float64}},cue_pattern::AbstractVector{Float64},pos::Integer,s::Float64,L::Integer)
	e1 = dot(W[:,pos],x)

	if L==2
		val = x[pos]
		idx = pos*2
		e2 = quaternion(cue_pattern[idx-1] * imagi(val),
						cue_pattern[idx-1] * real(val),
						cue_pattern[idx  ] * imagk(val),
						cue_pattern[idx  ] * imagj(val))
	elseif L==4
		val = x[pos]
		idx = pos*3
		e2 = quaternion(cue_pattern[idx-2] * imagi(val),
						cue_pattern[idx-2] * real(val),
						cue_pattern[idx-1] * real(val),
						cue_pattern[idx  ] * real(val))
	else
		@assert(false,"unimplemented")
	end

	return (e1 + s*e2)
end

function internal_potential(W::AbstractArray{Complex{Float64}}, x::AbstractVector{Complex{Float64}},cue_pattern::AbstractVector{Float64},pos::Integer,s::Float64,L::Integer)
	e1 = dot(W[:,pos],x)

	if L==2
		e2 = cue_pattern[pos] * complex(imag(x[pos]), real(x[pos]))
	else
		@assert(false,"unimplemented")
	end

	return (e1 + s*e2)
end

function internal_potential(W::AbstractArray{Float64}, x::AbstractVector{Float64},cue_pattern::AbstractVector{Float64},pos::Integer,s::Float64,L::Integer)
	e1 = dot(W[:,pos],x)

	if L==2
		n = length(cue_pattern)
		pos_cue = (pos -1)%n+1
		if pos > n
			e2 = cue_pattern[pos_cue] * x[pos_cue]
		else
			e2 = cue_pattern[pos_cue] * x[pos+n]
		end
	elseif L==4
		n = round(Int,length(x)/L)
		pos_cue = (pos -1)%n+1
		block = ifloor((pos-1)/n)
		@switch block begin
			0; e2 = cue_pattern[pos_cue] * x[pos_cue+1*n]
			1; e2 = cue_pattern[pos_cue] * x[pos_cue]
			2; e2 = cue_pattern[pos_cue+n] * x[pos_cue]
			3; e2 = cue_pattern[pos_cue+2*n] * x[pos_cue]
		end
	else
		@assert(false,"unimplemented")
	end

	return (e1 + s*e2)
end

function update_cue_async!(net::HopfieldNetwork, cue_pattern::AbstractVector,L::Integer, s::Real, beta::Real)
	ord = randperm(length(net.state))
	for i in ord
		h = internal_potential(net.weight,net.state,cue_pattern,i,s,L)
		net.state[i] = activatefunc(beta, h)
	end
end

function update_cue_sync!(net::HopfieldNetwork, cue_pattern::AbstractVector, L::Integer, s::Real, beta::Real)
	state = copy(net.state)
	for i = 1:length(net.state)
		h = internal_potential(net.weight,state,cue_pattern,i,s,L)
		net.state[i] = activatefunc(beta, h)
	end
end

function associate!{T1 <: Number, T2 <: Real}(net::HopfieldNetwork, init_pattern::AbstractVector{T1}, cue_pattern::AbstractVector{T2},L::Integer, s::Real, beta::Real, gamma::Real, iterations::Integer = 1_000, sync::Bool = false)
	copy!(net.state, init_pattern)
	associate!(net,cue_pattern,L,s,beta,gamma, iterations,sync)
end

function associate!{T <: Number}(net::HopfieldNetwork, cue_pattern::AbstractVector{T},L::Integer, s::Real, beta::Real, gamma::Real, iterations::Integer = 1_000, sync::Bool = false)
	if sync
		update! = update_cue_sync!
	else
		update! = update_cue_async!
	end
	beta_step = beta
	for i in 1:iterations
		update!(net, cue_pattern,L, s, beta_step)
		beta_step = beta_timestep(beta_step, gamma)
	end
	return copy(net.state)
end

function associate_step_async!{T <: Number}(net::HopfieldNetwork, cue_pattern::AbstractVector{T},L::Integer, s::Real, beta::Real, gamma::Real)
	update_cue_async!(net, cue_pattern, L, s, beta)
	return beta_timestep(beta, gamma)
end

function associate_step_sync!{T <: Number}(net::HopfieldNetwork, cue_pattern::AbstractVector{T},L::Integer, s::Real, beta::Real, gamma::Real)
	update_cue_sync!(net, cue_pattern, L, s, beta)
	return beta_timestep(beta, gamma)
end

function beta_timestep(beta::Real, gamma::Real)
	return gamma * beta;
end