text:
struct CutMesh
    mesh::Mesh
    cellsign::Vector{Int}
    activecells::Matrix{Bool}
    cutmeshnodeids::Matrix{Int}
    ncells::Int
    numnodes::Int
    nelmts::Int
    function CutMesh(mesh::Mesh, cellsign::Vector{Int}, cutmeshnodeids::Matrix{Int})
        ncells = number_of_cells(mesh)
        nummeshnodes = number_of_nodes(mesh)
        @assert length(cellsign) == ncells
        @assert size(cutmeshnodeids) == (2, nummeshnodes)

        activecells = active_cells(cellsign)
        numnodes = maximum(cutmeshnodeids)
        nelmts = count(activecells)
        new(mesh, cellsign, activecells, cutmeshnodeids, ncells, numnodes, nelmts)
    end
end

function active_cells(cellsign)
    ncells = length(cellsign)
    activecells = zeros(Bool, 2, ncells)
    idx1 = findall((cellsign .== +1) .| (cellsign .== 0))
    idx2 = findall((cellsign .== -1) .| (cellsign .== 0))
    activecells[1, idx1] .= true
    activecells[2, idx2] .= true
    return activecells
end

function CutMesh(levelset::InterpolatingPolynomial, levelsetcoeffs, mesh; tol = 1e-4, perturbation = 1e-3)
    nodalconnectivity = nodal_connectivity(mesh)
    cellsign = cell_sign!(levelset, levelsetcoeffs, nodalconnectivity, tol, perturbation)

    posactivenodeids = active_node_ids(+1, cellsign, nodalconnectivity)
    negactivenodeids = active_node_ids(-1, cellsign, nodalconnectivity)

    totalnumnodes = number_of_nodes(mesh)
    cutmeshnodeids = cut_mesh_nodeids(posactivenodeids, negactivenodeids, totalnumnodes)
    return CutMesh(mesh, cellsign, cutmeshnodeids)
end

function Base.show(io::IO, cutmesh::CutMesh)
    ncells = number_of_cells(cutmesh)
    nnodes = number_of_nodes(cutmesh)
    nelmts = number_of_elements(cutmesh)
    str = "CutMesh\n\tNum. Cells: $ncells\n\tNum. Elements: $nelmts\n\tNum. Nodes: $nnodes"
    print(io, str)
end

function nodal_connectivity(cutmesh::CutMesh, s, cellid)
    row = cell_sign_to_row(s)
    ncells = cutmesh.mesh.ncells
    @assert 1 <= cellid <= ncells
    cs = cell_sign(cutmesh, cellid)
    @assert cs == s || cs == 0

    nc = nodal_connectivity(cutmesh.mesh)
    ids = nc[:, cellid]
    nodeids = cutmesh.cutmeshnodeids[row, ids]
    return nodeids
end

function nodal_connectivity(cutmesh::CutMesh,cellid)
    return nodal_connectivity(cutmesh.mesh,cellid)
end

function active_cells(cutmesh::CutMesh)
    return cutmesh.activecells
end

function dimension(cutmesh::CutMesh)
    return dimension(cutmesh.mesh)
end

function number_of_nodes(cutmesh::CutMesh)
    return cutmesh.numnodes
end

function reference_corner(cutmesh::CutMesh)
    return reference_corner(cutmesh.mesh)
end

function nodes_per_mesh_side(cutmesh::CutMesh)
    return nodes_per_mesh_side(cutmesh.mesh)
end

function widths(cutmesh::CutMesh)
    return widths(cutmesh.mesh)
end

function number_of_degrees_of_freedom(cutmesh::CutMesh)
    dim = dimension(cutmesh)
    numnodes = number_of_nodes(cutmesh)
    return dim * numnodes
end

function number_of_cells(cutmesh::CutMesh)
    return cutmesh.ncells
end

function background_mesh(cutmesh::CutMesh)
    return cutmesh.mesh
end

function number_of_elements(cutmesh::CutMesh)
    return cutmesh.nelmts
end

function cell_sign(cutmesh::CutMesh)
    return cutmesh.cellsign
end

function cell_id(cutmesh::CutMesh,x)
    return cell_id(cutmesh.mesh,x)
end

function cell_sign(cutmesh::CutMesh, cellid)
    return cutmesh.cellsign[cellid]
end

function cell_map(cutmesh::CutMesh, cellid)
    return cell_map(cutmesh.mesh, cellid)
end

function cell_map(cutmesh::CutMesh)
    return cell_map(cutmesh,1)
end

function cell_map(cutmesh::CutMesh, s, cellid)
    return cell_map(cutmesh, cellid)
end

function jacobian(cutmesh::CutMesh)
    return jacobian(cell_map(cutmesh))
end

function inverse_jacobian(cutmesh::CutMesh)
    return 1.0 ./ jacobian(cutmesh)
end

function face_determinant_jacobian(cutmesh::CutMesh)
    return face_determinant_jacobian(cutmesh.mesh)
end

function is_interior_cell(cutmesh::CutMesh)
    return is_interior_cell(cutmesh.mesh)
end

function cell_connectivity(cutmesh::CutMesh)
    return cell_connectivity(cutmesh.mesh)
end

function nodal_coordinates(cutmesh::CutMesh)
    return nodal_coordinates(cutmesh.mesh)
end

function cell_sign!(levelset, levelsetcoeffs, nodalconnectivity, tol, perturbation)
    ncells = size(nodalconnectivity)[2]
    cellsign = zeros(Int, ncells)
    xL,xR = [-1.,-1.],[1.,1.]
    for cellid = 1:ncells
        nodeids = nodalconnectivity[:, cellid]
        update!(levelset, levelsetcoeffs[nodeids])

        s = sign(levelset, xL, xR, tol=tol)
        if (s == +1 || s == 0 || s == -1)
            cellsign[cellid] = s
        else
            @warn "Perturbing levelset function by perturbation = $perturbation"
            levelsetcoeffs[nodeids] .+= perturbation
            update!(levelset,levelsetcoeffs[nodeids])

            s = sign(levelset, xL, xR, tol=tol)
            if (s == +1 || s == 0 || s == -1)
                cellsign[cellid] = s
            else
                error("Could not determine cell sign after perturbation = $perturbation")
            end
        end
    end
    return cellsign
end

function active_node_ids(s, cellsign, nodalconnectivity)
    @assert (s == -1 || s == +1)
    activenodeids = Int[]
    numcells = size(nodalconnectivity)[2]
    for cellid = 1:numcells
        if cellsign[cellid] == s || cellsign[cellid] == 0
            append!(activenodeids, nodalconnectivity[:, cellid])
        end
    end
    sort!(activenodeids)
    unique!(activenodeids)
    return activenodeids
end

function cut_mesh_nodeids(posactivenodeids, negactivenodeids, totalnumnodes)
    cutmeshnodeids = zeros(Int, 2, totalnumnodes)
    counter = 1
    for nodeid in posactivenodeids
        cutmeshnodeids[1, nodeid] = counter
        counter += 1
    end
    for nodeid in negactivenodeids
        cutmeshnodeids[2, nodeid] = counter
        counter += 1
    end
    return cutmeshnodeids
end


text:
# SPDX-License-Identifier: X11
# 2020-11-15 (Live)
# Travel (300pt)

function main()
  N, K = parse.(Int, (split ∘ readline)())
  T = Matrix{Int}(undef, N, N)
  for i ∈ 1:N
    t = parse.(Int, (split ∘ readline)())
    for j ∈ 1:N
      T[i, j] = t[j]
    end
  end

  (println ∘ length ∘
    (x -> filter!(y -> y == K, x)))(map(x -> pathl(T, x),
                                             (permutations ∘ collect)(2:N)))
end

function pathl(T::Matrix{Int}, p::Vector{Int})::Int
  sum = T[1, p[1]]
  for i ∈ 1:lastindex(p) - 1
    sum += T[p[i], p[i + 1]]
  end
  sum += T[p[end], 1]
end

function permutations(a::Vector{Int})::Vector{Vector{Int}}
  length(a) == 1 && return [a]

  ret = Vector{Vector{Int}}()
  for i ∈ 1:lastindex(a)
    for j ∈ permutations([a[1:i - 1]; a[i + 1:end]])
      push!(ret, [a[i]; j])
    end
  end

  return ret
end

main()

ext: jl

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
## -*-Julia-*-
## Test suite for Julia's WAV module
reload("WAV.jl")

# These float array comparison functions are from dists.jl
function absdiff{T<:Real}(current::AbstractArray{T}, target::AbstractArray{T})
    @assert all(size(current) == size(target))
    maximum(abs(current - target))
end

function reldiff{T<:Real}(current::T, target::T)
    abs((current - target)/(bool(target) ? target : 1))
end

function reldiff{T<:Real}(current::AbstractArray{T}, target::AbstractArray{T})
    @assert all(size(current) == size(target))
    maximum([reldiff(current[i], target[i]) for i in 1:numel(target)])
end

## example from README, modified to use an IO buffer
let
    x = [0:7999]
    y = sin(2 * pi * x / 8000)
    io = IOBuffer()
    WAV.wavwrite(y, io, Fs=8000)
    seek(io, 0)
    y, Fs = WAV.wavread(io)
end

## default arguments, GitHub Issue #10
let
    tmp=rand(Float32,(10,2))
    io = IOBuffer()
    WAV.wavwrite(tmp, io; nbits=32)
    seek(io, 0)
    y, fs, nbits, extra = WAV.wavread(io; format="native")
    @assert typeof(y) == Array{Float32, 2}
    @assert fs == 8000.0
    @assert nbits == 32
end

let
    tmp=rand(Float64,(10,2))
    io = IOBuffer()
    WAV.wavwrite(tmp, io; nbits=64)
    seek(io, 0)
    y, fs, nbits, extra = WAV.wavread(io; format="native")
    @assert typeof(y) == Array{Float64, 2}
    @assert fs == 8000.0
    @assert nbits == 64
end

let
    tmp=rand(Float32,(10,2))
    io = IOBuffer()
    WAV.wavwrite(tmp, io; compression=WAV.WAVE_FORMAT_IEEE_FLOAT)
    seek(io, 0)
    y, fs, nbits, extra = WAV.wavread(io; format="native")
    @assert typeof(y) == Array{Float32, 2}
    @assert fs == 8000.0
    @assert nbits == 32
end

## Test wavread and wavwrite
## Generate some wav files for writing and reading
for fs = (8000,11025,22050,44100,48000,96000,192000), nbits = (1,7,8,9,12,16,20,24,32,64), nsamples = convert(Array{Int}, [0, logspace(1, 4, 4)]), nchans = 1:4
    ## Test wav files
    ## The tolerance is based on the number of bits used to encode the file in wavwrite
    tol = 2.0 / (2.0^nbits - 1)

    in_data = rand(nsamples, nchans)
    if nsamples > 0
        @assert maximum(in_data) <= 1.0
        @assert minimum(in_data) >= -1.0
    end
    io = IOBuffer()
    WAV.wavwrite(in_data, io, Fs=fs, nbits=nbits, compression=WAV.WAVE_FORMAT_PCM)
    file_size = position(io)

    ## Check for the common header identifiers
    seek(io, 0)
    @assert read(io, Uint8, 4) == b"RIFF"
    @assert WAV.read_le(io, Uint32) == file_size - 8
    @assert read(io, Uint8, 4) == b"WAVE"

    ## Check that wavread works on the wavwrite produced memory
    seek(io, 0)
    sz = WAV.wavread(io, format="size")
    @assert sz == (nsamples, nchans)

    seek(io, 0)
    out_data, out_fs, out_nbits, out_extra = WAV.wavread(io)
    @assert length(out_data) == nsamples * nchans
    @assert size(out_data, 1) == nsamples
    @assert size(out_data, 2) == nchans
    @assert typeof(out_data) == Array{Float64, 2}
    @assert out_fs == fs
    @assert out_nbits == nbits
    @assert out_extra == None
    if nsamples > 0
        @assert absdiff(out_data, in_data) < tol
    end

    ## test the "subrange" option.
    if nsamples > 0
        seek(io, 0)
        # Don't convert to Int, test if passing a float (nsamples/2) behaves as expected
        subsamples = min(10, int(nsamples / 2))
        out_data, out_fs, out_nbits, out_extra = WAV.wavread(io, subrange=subsamples)
        @assert length(out_data) == subsamples * nchans
        @assert size(out_data, 1) == subsamples
        @assert size(out_data, 2) == nchans
        @assert typeof(out_data) == Array{Float64, 2}
        @assert out_fs == fs
        @assert out_nbits == nbits
        @assert out_extra == None
        @assert absdiff(out_data, in_data[1:int(subsamples), :]) < tol

        seek(io, 0)
        sr = convert(Int, min(5, int(nsamples / 2))):convert(Int, min(23, nsamples - 1))
        out_data, out_fs, out_nbits, out_extra = WAV.wavread(io, subrange=sr)
        @assert length(out_data) == length(sr) * nchans
        @assert size(out_data, 1) == length(sr)
        @assert size(out_data, 2) == nchans
        @assert typeof(out_data) == Array{Float64, 2}
        @assert out_fs == fs
        @assert out_nbits == nbits
        @assert out_extra == None
        @assert absdiff(out_data, in_data[sr, :]) < tol
    end
end

## Test native encoding of 8 bits
for nchans = (1,2,4)
    in_data_8 = uint8(reshape(int(typemin(Uint8):typemax(Uint8)), (int(256 / nchans), nchans)))
    io = IOBuffer()
    WAV.wavwrite(in_data_8, io)

    seek(io, 0)
    out_data_8, fs, nbits, extra = WAV.wavread(io, format="native")
    @assert fs == 8000
    @assert nbits == 8
    @assert extra == None
    @assert in_data_8 == out_data_8
end

## Test native encoding of 16 bits
for nchans = (1,2,4)
    in_data_16 = int16(reshape(int(typemin(Int16):typemax(Int16)), (int(65536 / nchans), nchans)))
    io = IOBuffer()
    WAV.wavwrite(in_data_16, io)

    seek(io, 0)
    out_data_16, fs, nbits, extra = WAV.wavread(io, format="native")
    @assert fs == 8000
    @assert nbits == 16
    @assert extra == None
    @assert in_data_16 == out_data_16
end

## Test native encoding of 24 bits
for nchans = (1,2,4)
    in_data_24 = convert(Array{Int32}, reshape(-63:64, int(128 / nchans), nchans))
    io = IOBuffer()
    WAV.wavwrite(in_data_24, io)

    seek(io, 0)
    out_data_24, fs, nbits, extra = WAV.wavread(io, format="native")
    @assert fs == 8000
    @assert nbits == 24
    @assert extra == None
    @assert in_data_24 == out_data_24
end

## Test encoding 32 bit values
for nchans = (1,2,4)
    in_data_single = convert(Array{Float32}, reshape(linspace(-1.0, 1.0, 128), int(128 / nchans), nchans))
    io = IOBuffer()
    WAV.wavwrite(in_data_single, io)

    seek(io, 0)
    out_data_single, fs, nbits, extra = WAV.wavread(io, format="native")
    @assert fs == 8000
    @assert nbits == 32
    @assert extra == None
    @assert in_data_single == out_data_single
end

## Test encoding 32 bit values outside the valid range
for nchans = (1,2,4)
    nsamps = int(128 / nchans)
    in_data_single = convert(Array{Float32}, reshape(-63:64, nsamps, nchans))
    io = IOBuffer()
    WAV.wavwrite(in_data_single, io)

    seek(io, 0)
    out_data_single, fs, nbits, extra = WAV.wavread(io, format="native")
    @assert fs == 8000
    @assert nbits == 32
    @assert extra == None
    @assert [clamp(in_data_single[i, j], float32(-1), float32(1)) for i = 1:nsamps, j = 1:nchans] == out_data_single
end

### Test A-Law and Mu-Law
for nbits = (8, 16), nsamples = convert(Array{Int}, [0, logspace(1, 4, 4)]), nchans = 1:2, fmt=(WAV.WAVE_FORMAT_ALAW, WAV.WAVE_FORMAT_MULAW)
    const fs = 8000.0
    const tol = 2.0 / (2.0^6)
    in_data = rand(nsamples, nchans)
    if nsamples > 0
        @assert maximum(in_data) <= 1.0
        @assert minimum(in_data) >= -1.0
    end
    io = IOBuffer()
    WAV.wavwrite(in_data, io, Fs=fs, nbits=nbits, compression=fmt)
    file_size = position(io)

    ## Check for the common header identifiers
    seek(io, 0)
    @assert read(io, Uint8, 4) == b"RIFF"
    @assert WAV.read_le(io, Uint32) == file_size - 8
    @assert read(io, Uint8, 4) == b"WAVE"

    ## Check that wavread works on the wavwrite produced memory
    seek(io, 0)
    sz = WAV.wavread(io, format="size")
    @assert sz == (nsamples, nchans)

    seek(io, 0)
    out_data, out_fs, out_nbits, out_extra = WAV.wavread(io)
    @assert length(out_data) == nsamples * nchans
    @assert size(out_data, 1) == nsamples
    @assert size(out_data, 2) == nchans
    @assert typeof(out_data) == Array{Float64, 2}
    @assert out_fs == fs
    @assert out_nbits == 8
    @assert out_extra == None
    if nsamples > 0
        @assert absdiff(out_data, in_data) < tol
    end

    ## test the "subrange" option.
    if nsamples > 0
        seek(io, 0)
        # Don't convert to Int, test if passing a float (nsamples/2) behaves as expected
        subsamples = min(10, int(nsamples / 2))
        out_data, out_fs, out_nbits, out_extra = WAV.wavread(io, subrange=subsamples)
        @assert length(out_data) == subsamples * nchans
        @assert size(out_data, 1) == subsamples
        @assert size(out_data, 2) == nchans
        @assert typeof(out_data) == Array{Float64, 2}
        @assert out_fs == fs
        @assert out_nbits == 8
        @assert out_extra == None
        @assert absdiff(out_data, in_data[1:int(subsamples), :]) < tol

        seek(io, 0)
        sr = convert(Int, min(5, int(nsamples / 2))):convert(Int, min(23, nsamples - 1))
        out_data, out_fs, out_nbits, out_extra = WAV.wavread(io, subrange=sr)
        @assert length(out_data) == length(sr) * nchans
        @assert size(out_data, 1) == length(sr)
        @assert size(out_data, 2) == nchans
        @assert typeof(out_data) == Array{Float64, 2}
        @assert out_fs == fs
        @assert out_nbits == 8
        @assert out_extra == None
        @assert absdiff(out_data, in_data[sr, :]) < tol
    end
end

### Test float formatting
for nbits = (32, 64), nsamples = convert(Array{Int}, [0, logspace(1, 4, 4)]), nchans = 1:2, fmt=(WAV.WAVE_FORMAT_IEEE_FLOAT)
    const fs = 8000.0
    const tol = 1e-6
    in_data = rand(nsamples, nchans)
    if nsamples > 0
        @assert maximum(in_data) <= 1.0
        @assert minimum(in_data) >= -1.0
    end
    io = IOBuffer()
    WAV.wavwrite(in_data, io, Fs=fs, nbits=nbits, compression=fmt)
    file_size = position(io)

    ## Check for the common header identifiers
    seek(io, 0)
    @assert read(io, Uint8, 4) == b"RIFF"
    @assert WAV.read_le(io, Uint32) == file_size - 8
    @assert read(io, Uint8, 4) == b"WAVE"

    ## Check that wavread works on the wavwrite produced memory
    seek(io, 0)
    sz = WAV.wavread(io, format="size")
    @assert sz == (nsamples, nchans)

    seek(io, 0)
    out_data, out_fs, out_nbits, out_extra = WAV.wavread(io)
    @assert length(out_data) == nsamples * nchans
    @assert size(out_data, 1) == nsamples
    @assert size(out_data, 2) == nchans
    @assert typeof(out_data) == Array{Float64, 2}
    @assert out_fs == fs
    @assert out_nbits == nbits
    @assert out_extra == None
    if nsamples > 0
        @assert absdiff(out_data, in_data) < tol
    end

    ## test the "subrange" option.
    if nsamples > 0
        seek(io, 0)
        # Don't convert to Int, test if passing a float (nsamples/2) behaves as expected
        subsamples = min(10, int(nsamples / 2))
        out_data, out_fs, out_nbits, out_extra = WAV.wavread(io, subrange=subsamples)
        @assert length(out_data) == subsamples * nchans
        @assert size(out_data, 1) == subsamples
        @assert size(out_data, 2) == nchans
        @assert typeof(out_data) == Array{Float64, 2}
        @assert out_fs == fs
        @assert out_nbits == nbits
        @assert out_extra == None
        @assert absdiff(out_data, in_data[1:int(subsamples), :]) < tol

        seek(io, 0)
        sr = convert(Int, min(5, int(nsamples / 2))):convert(Int, min(23, nsamples - 1))
        out_data, out_fs, out_nbits, out_extra = WAV.wavread(io, subrange=sr)
        @assert length(out_data) == length(sr) * nchans
        @assert size(out_data, 1) == length(sr)
        @assert size(out_data, 2) == nchans
        @assert typeof(out_data) == Array{Float64, 2}
        @assert out_fs == fs
        @assert out_nbits == nbits
        @assert out_extra == None
        @assert absdiff(out_data, in_data[sr, :]) < tol
    end
end


text:
module GDAL

using Libdl
using CEnum

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("GDAL not installed properly, run Pkg.build(\"GDAL\"), restart Julia and try again")
end
include(depsjl_path)

const Ctm = Base.Libc.TmStruct

include("common.jl")

include("gdal_h.jl")
include("gdal_alg.jl")
include("gdal_vrt.jl")
include("gdal_utils.jl")
include("gdalwarper.jl")

include("ogr_api.jl")
include("ogr_srs_api.jl")
include("ogr_core.jl")

include("cpl_error.jl")
include("cpl_conv.jl")
include("cpl_minixml.jl")
include("cpl_progress.jl")
include("cpl_virtualmem.jl")
include("cpl_vsi.jl")

include("error.jl")

const GDALVERSION = Ref{VersionNumber}()
const GDAL_DATA = Ref{String}()
const PROJ_LIB = Ref{String}()

function __init__()
    # Always check your dependencies from `deps.jl`
    check_deps()

    # register custom error handler
    funcptr = @cfunction(gdaljl_errorhandler, Ptr{Cvoid}, (CPLErr, Cint, Cstring))
    cplseterrorhandler(funcptr)

    # get GDAL version number
    versionstring = gdalversioninfo("RELEASE_NAME")
    GDALVERSION[] = VersionNumber(versionstring)

    # set GDAL_DATA location, this overrides setting the environment variable
    GDAL_DATA[] = abspath(@__DIR__, "..", "deps", "usr", "share", "gdal")
    cplsetconfigoption("GDAL_DATA", GDAL_DATA[])

    # set PROJ_LIB location, this overrides setting the environment variable
    PROJ_LIB[] = abspath(@__DIR__, "..", "deps", "usr", "share", "proj")
    osrsetprojsearchpaths([PROJ_LIB[]])
end

end # module


text:
module OAuth

using HTTP, MbedTLS, Base64, Random

export
    oauth_timestamp,
    oauth_nonce,
    oauth_sign_hmac_sha1,
    oauth_signing_key,
    oauth_signature_base_string,
    oauth_percent_encode_keys!,
    oauth_serialize_url_parameters,
    encodeURI!,
    oauth_body_hash_file,
    oauth_body_hash_data,
    oauth_body_hash_encode,
    oauth_header,
    oauth_request_resource,
    OAuth1a,
    OAuth2

abstract type AbstractOAuth end

include("common_functions.jl")
include("oauth1a.jl")
include("oauth2.jl")

end


text:
"""
    minishow(key, value)
Show an element of tuple as simple format.

# Examples
```julia
julia> A = (a = "1", b = 1, c = 1.0, d = -1.0)
julia> julia> Equate.minishow.(keys(A), values(A))
         a :          1
         b :        1
         c : 1.00000
         d : -1.00000
(nothing, nothing, nothing, nothing)
```
"""

function minishow(io::IO, k, v)
    if v isa Float64
        @printf io "%10s : %3.5f \n" k v
    elseif v isa Int64
        @printf io "%10s : %8i \n" k v
    elseif v isa String
        @printf io "%10s : %10s \n" k v
    else
        @printf io "%10s : can't be shown! \n" k
    end
end

function minishow(io::IO, A)
    minishow.(Ref(io), keys(A), values(A))
end

function Base.show(io::IO, x::FreqTab)
    println(io, "Frequency table stats.")
    minishow(io, x.stats)
    return
end

function Base.show(io::IO, x::NEATFreqTab)
    println(io, "Frequency table stats.")
    println(io, "* X (independent part)")
    minishow(io, x.statsX)
    println(io, "* V (common part)")
    minishow(io, x.statsV)
    return
end

function Base.show(io::IO, x::NEATEquateResult)
    println(io, "Equating design: NEAT\nEquated method: $(x.method).\nTo show the table, extract `table` element.")
    println(io, x.synthetic)
    return
end

function Base.show(io::IO, x::SGEquateResult)
    println(io, "Equating design: EG\nEquated method: $(x.method).\nTo show the table, extract `table` element.")
    return
end


text:
module DiffEqPhysics

using Reexport
@reexport using DiffEqBase, RecursiveArrayTools
using ForwardDiff, StaticArrays, RecipesBase
using Random, Printf, LinearAlgebra

include("nbody_problem.jl")
include("hamiltonian.jl")
include("plot.jl")

export HamiltonianProblem, NBodyProblem, plot_orbits

export NBodySimulation
export MassBody, ChargedParticle, MagneticParticle
export PotentialParameters, LennardJonesParameters, GravitationalParameters,
       ElectrostaticParameters, MagnetostaticParameters, SPCFwParameters
export PotentialNBodySystem, ChargedParticles, GravitationalSystem, WaterSPCFw
export PeriodicBoundaryConditions, CubicPeriodicBoundaryConditions, InfiniteBox
export AndersenThermostat, BerendsenThermostat, NoseHooverThermostat, LangevinThermostat
export run_simulation, get_position, get_velocity, get_masses, temperature,
       initial_energy, kinetic_energy, potential_energy, total_energy, rdf, msd,
       generate_bodies_in_cell_nodes, run_simulation_sde

end # module


text:
#=
Copyright 2017 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

# This file contains all the functions used within `openBF` to convert
# quantities. In particular, the following functions are provided
#
# *   [`waveSpeed`](converter.html#waveSpeed): $A \rightarrow c$;
# *   [`pressure`](converter.html#pressure): $A \rightarrow P$;
# *   [`areaFromPressure`](converter.html#areaFromPressure):
#       $P \rightarrow A$;
# *   [`riemannInvariants`](converter.html#riemannInvariants):
#       $(u, c) \rightarrow (W_1, W_2)$;
# *   [`rI2uc`](converter.html#rI2uc): $(W_1, W_2) \rightarrow (u, c)$.
#

# Wave speed
# $$
#   c = \sqrt{\frac{3}{2}\gamma \sqrt{A}}, \quad
#     \gamma = \frac{\beta}{3 \rho R_0 \sqrt{\pi}}.
# $$

# *function* __`waveSpeed`__ $\rightarrow$ `c::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `A`              `::Float` cross sectional area.
#
# `gamma`          `::Float` elastic constant.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ------------- --------------------------------------------------------------
# `c`           `::Float` wave speed given cross sectional area A.
# ----------------------------------------------------------------------------
#
function waveSpeed(A :: Float64, gamma :: Float64)

  return sqrt(3*gamma*sqrt(A)*0.5)

end

# *function* __`waveSpeed`__ $\rightarrow$ `c::Array{Float, 1}`
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `A`              `::Array{Float, 1}` cross sectional area.
#
# `gamma`          `::Float` elastic constant.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# To compute the wave speed along the entire vessel, `waveSpeed` is
# re-defined to handle also a vector of `A`s. `waveSpeed` is called
# recursively with the same `gamma` along the vessel.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ------------- --------------------------------------------------------------
# `c`           `::Arrary{Float, 1}` wave speed given cross sectional area A.
# ----------------------------------------------------------------------------
#
function waveSpeed(A :: Array{Float64, 1}, gamma :: Array{Float64, 1}, c :: Array{Float64, 1})

  for i in 1:length(A)
    c[i] = waveSpeed(A[i], gamma[i])
  end

  return c
end

# Trans-mural pressure
# $$
#   P = P_{ext} + \beta \left[\left(\frac{A}{A_0}\right)^{1/2} -1 \right], \quad
#     \beta = \sqrt{\frac{\pi}{A_0}} \frac{h_0 E}{1 - \sigma^2}.
# $$

# *function* __`pressure`__ $\rightarrow$ `P::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `A`              `::Float` cross sectional area.
#
# `A0`             `::Float` unstressed cross sectional area.
#
# `beta`           `::Float` elastic constant.
#
# `Pext`           `::Float` constant external pressure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ------------- --------------------------------------------------------------
# `P`           `::Float` trans-mural pressure given the cross sectional area
#               `A`.
# ----------------------------------------------------------------------------
#
function pressure(A    :: Float64, A0   :: Float64,
                  beta :: Float64, Pext :: Float64)

  return Pext + beta*(sqrt(A/A0) - 1.)

end

# *function* __`pressure`__ $\rightarrow$ `P::Array{Float, 1}`
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `A`              `::Array{Float, 1}` cross sectional area.
#
# `A0`             `::Float` unstressed cross sectional area.
#
# `beta`           `::Float` elastic constant.
#
# `Pext`           `::Float` constant external pressure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# To compute the pressure along the entire vessel, `pressure` is
# re-defined to handle also a vector of `A`s. `pressure` is called
# recursively with the same `beta` along the vessel.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ------------- --------------------------------------------------------------
# `P`           `::Array{Float, 1}` trans-mural pressure given
#               the cross sectional area `A`.
# ----------------------------------------------------------------------------
# function pressure(A    :: Array{Float64, 1}, A0   :: Float64,
#                   beta :: Float64,           Pext :: Float64)

#   p = zeros(Float64, length(A))

#   for i in 1:length(A)
#     p[i] = Pext + beta*(sqrt(A[i]/A0) - 1.)

#   end

#   return p
# end

function pressure(A    :: Array{Float64, 1}, A0   :: Array{Float64, 1},
                  beta :: Array{Float64,1},  Pext :: Float64,
                  p :: Array{Float64, 1})

  for i in 1:length(A)
    p[i] = Pext + beta[i]*(sqrt(A[i]/A0[i]) - 1.)
  end

  return p
end

# Trans-mural pressure can be inverted to find cross sectional area as
# $$
#   A = A_0 \left( \frac{P-P_{ext}}{\beta} +1 \right)^2.
# $$

# *function* __`areaFromPressure`__ $\rightarrow$ `A::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `P`              `::Float` cross sectional area.
#
# `A0`             `::Float` unstressed cross sectional area.
#
# `beta`           `::Float` elastic constant.
#
# `Pext`           `::Float` constant external pressure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ------------- --------------------------------------------------------------
# `A`           `::Float` cross sectional area given the trans-mural pressure
#               `P`.
# ----------------------------------------------------------------------------
#
function areaFromPressure(P    :: Float64, A0   :: Float64,
                          beta :: Float64, Pext :: Float64)

   return A0 * ((P-Pext)/beta + 1)*((P-Pext)/beta + 1)

end

# *function* __`areaFromPressure`__ $\rightarrow$ `A::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `P`              `::Float` cross sectional area.
#
# `A0`             `::Float` unstressed cross sectional area.
#
# `beta`           `::Float` elastic constant.
#
# `Pext`           `::Float` constant external pressure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# To compute the cross sectional area along the entire vessel,
# `areaFromPressure` is re-defined to handle also a vector of `P`s.
# `areaFromPressure` is called recursively with the same `beta` along the
# vessel.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ------------- --------------------------------------------------------------
# `A`           `::Float` cross sectional area given the trans-mural pressure
#               `P`.
# ----------------------------------------------------------------------------
# function areaFromPressure(P    :: Array{Float64, 1}, A0   :: Float64,
#                           beta :: Float64,           Pext :: Float64)

#   a = zeros(Float64, length(P))

#   for i in 1:length(P)
#     a[i] = areaFromPressure(P[i], A0, beta)

#   end

#   return a
# end

# Riemann invariants are computed from `u` and `c` values as
# $$
#   W_1 = u - 4c, \quad W_2 = u + 4c.
# $$

# *function* __`riemannInvariants`__ $\rightarrow$ `W1::Float`, `w2::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `i`              `::Float` cell index.
#
# `v`              `::Vessel` vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ------------- --------------------------------------------------------------
# `W1`          `::Float` backward Riemann invariant.
#
# `W2`          `::Float` forward Riemann invariant.
# ----------------------------------------------------------------------------
#
function riemannInvariants(i :: Int64, v :: Vessel)

  W1 = v.u[i] - 4*v.c[i]
  W2 = v.u[i] + 4*v.c[i]

  return W1, W2
end

# `u` and `c` can be computed from Riemann invariants as
# $$
#   u = \frac{1}{2} (W_1 + w_2), \quad c = \frac{W_2 - W_1}{8}.
# $$

# *function* __`rI2uc`__ $\rightarrow$ `W1::Float`, `w2::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `i`              `::Float` cell index.
#
# `v`              `::Vessel` vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ------------- --------------------------------------------------------------
# `u`           `::Float` longitudinal velocity.
#
# `c`           `::Float` wave speed.
# ----------------------------------------------------------------------------
#
function rI2uc(W1 :: Float64, W2 :: Float64)
  u = 0.5*(W1 + W2)
  c = (W2 - W1)*0.125

  return u, c
end

# *function* __`calculatePrimitives`__
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `v`              `::Vessel` vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# This function updates all primitives variables at once.
# ----------------------------------------------------------------------------
#
# function calculatePrimitives(v :: Vessel)

#   v.P = pressure(v.A, v.A0, v.beta)
#   v.u = v.Q./v.A
#   v.c = waveSpeed(v.A, v.gamma)

# end


text:
using Distributions
using Base.Test

my_tests = ["test/distributions.jl",
            "test/utils.jl",
            "test/wisharts.jl",
            "test/fit.jl",
            "test/truncate.jl",
            "test/multivariate.jl"]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
using Distributions: cdf, FDist, TDist, Chisq, Bernoulli
using StatsBase: sample
using DataFrames: DataFrame
using StatsModels: @formula
using GLM: glm, loglikelihood, coef, confint, predict

"""
    conditional_inv_test_chow(X, y, env, n_env; α=0.01)

Test H0: `y` | `X` is invariant under environments by combining leave-one-out chow test.

* `X`:                 n x p matrix of covariates
* `y`:                 vector of response
* `env`:               vector of environments (1:n_env), corresponding to rows of X
* `n_env`:             number of environments
* `α`:                 significance level
* `n_max_for_exact`:   maximum number of observations of an environment for exact testing;
                      otherwise a subsample of n_max rows will be used

Return: `rej`, `p_value`, `conf_intervals`
* `rej`: false if invariant
* `conf_intervals`: p x 2 matrix, cols = (min, max)
"""
function conditional_inv_test_chow(X::Matrix{Float64}, y::Vector{Float64}, env::Vector{Int64}, n_env::Int64;
                                     α=0.01, n_max_for_exact=5000)
    @assert n_env >= 1
    p_values = ones(n_env)
    p = size(X, 2)
    n_all = size(X, 1)
    for i in 1:n_env
        @assert sum(env.==i) > 0
        # fit on (-env), test on (env)
        idx_in = (1:n_all)[env.==i]
        idx_out = (1:n_all)[env.!=i]
        if length(idx_in) > n_max_for_exact
            idx_in = sample(idx_in, n_max_for_exact, replace=false)
        end
        p_values[i] = two_sample_chow(X[idx_out, :], X[idx_in, :], y[idx_out], y[idx_in];
                                      α=α, add_intercept=true)
        if p_values[i] < min(α / n_env)
            break   # early termination
        end
    end
    # Bonferroni correction
    p_value = min(minimum(p_values) * n_env, 1.)
    if p_value < α
        # reject
        reject = true
        conf_intervals = zeros(p, 2)  # min = max = 0.
    else
        # not reject
        reject = false
        # pool all data and run regression
        X = hcat(X, ones(size(X,1)))
        β = (X' * X) \ (X' * y)
        σ = sqrt(sum((y - X * β).^2) / (n_all - p - 1))
        prec_X = diag(inv(X' * X))
        qt = quantile(TDist(n_all - p - 1), 1 - α / (2 * p))
        conf_left = β - qt * σ * sqrt.(prec_X)
        conf_right = β + qt * σ * sqrt.(prec_X)
        # note: the union of (1-α) conf ints have coverage (1-2α) (Theorem 2 of PBM)
        conf_intervals = hcat(conf_left[1:p], conf_right[1:p])
    end
    return reject, p_value, conf_intervals
end

"""
    two_sample_chow(X1, X2, y1, y2; [α=0.01, add_intercept=true])

Two-sample Chow's test of H0: two linear regressions `y1 ~ X1` and `y2 ~ X2` have
the same linear coefficients, assuming Gaussian errors with equal variances.

Will fit on `X1`` and test residual on `X2`.
Choose `X2` to be the sample with smaller n since n2 x n2 matrix inversion is performed.
"""
function two_sample_chow(X1::Matrix{Float64}, X2::Matrix{Float64},
                         y1::Vector{Float64}, y2::Vector{Float64}; α=0.01, add_intercept=true)
    n1 = size(X1, 1)
    n2 = size(X2, 1)
    if add_intercept
        X1 = hcat(X1, ones(n1))
        X2 = hcat(X2, ones(n2))
    end
    p = size(X1, 2)
    # fit on X1
    β = zeros(p)
    try
        β = (X1' * X1) \ (X1' * y1)
    catch _err
        printstyled("encountered $_err in least square for Chow's test\n", color=:light_red)
        return 0.
    end
    res2 = y2 - X2 * β
    Σ_res = diagm(ones(n2)) + X2 * ((X1' * X1) \ X2')   # inv(A) * B = A \ B
    σ2 = var(y1 - X1 * β) * (n1 - 1) / (n1 - p)  # should use dof = (n - p) as denominator
    chow_stat = res2' * (Σ_res \ res2) / (σ2 * n2)   # inv(A) * B = A \ B
    # F distribution
    ν1 = n2
    ν2 = n1 - size(X1, 2)
    p_value = 1 - cdf(FDist(ν1, ν2), chow_stat)
    return p_value
end

"""
    conditional_inv_test_logistic(df, :target, S, env, n_env; α=0.01, add_intercept=true, method="logistic-LR")

Test `H0: y | X` is invariant under environments by combining leave-one-out likelihood ratio test.
The model is logistic regression specified by `fmla`.

* `df`:                a DataFrames.DataFrame
* `target`:            target variable in df (Symbol)
* `S`:                 covariates to condition on (`X = df[:,S]`)
* `env`:               vector of environments (1:n_env), corresponding to rows of X
* `n_env`:             number of environments
* `α`:                 significance level
* `add_intercept`:     add `+ 1` or not in formula
* `method`:
    + `logistic-LR`: likelihood ratio test
    + `logistic-SF`: test equal mean and variance of prediction errors with Sukhatme-Fisher

Return: `rej`, `p_value`, `conf_intervals`
* `rej`: false if invariant
* `conf_intervals`: p x 2 matrix, cols = (min, max)
"""
function conditional_inv_test_logistic(df::DataFrame, target::Symbol, S::Vector{Symbol},
                                          env::Vector{Int64}, n_env::Int64; α=0.01, add_intercept=true, method="logistic-LR")
    @assert n_env >= 1
    @assert !(target in S)
    p = length(S) + add_intercept
    n_all = size(df, 1)
    p_values = ones(n_env)
    fmla = get_formula(df, target, S, add_intercept=add_intercept)
    # fit on pooled data
    fit0 = glm(fmla, df, Bernoulli())
    # iterate over environments
    for i in 1:n_env
        @assert sum(env.==i) > 0
        if method == "logistic-LR"
            # fit separately
            fit1 = glm(fmla, df[env.==i, :], Bernoulli())
            fit2 = glm(fmla, df[env.!=i, :], Bernoulli())
            # log likelihood ratio = 2 log (p(bigger model) / p(smaller model))
            lr = 2 * (loglikelihood(fit1) + loglikelihood(fit2) - loglikelihood(fit0))
            p_values[i] = 1 - cdf(Chisq(p), lr)
        elseif method == "logistic-SF"
            # predict and test equal mean of residuals
            if p > 1
                p_hat = predict(fit0, df)
            else
                p_hat = predict(fit0)
            end
            res = (df[!, target] - p_hat) ./ sqrt.(p_hat .* (1 .- p_hat))
            p_values[i] = sukhatme_fisher_test(res[env.==i], res[env.!=i])
        else
            error("method undefined")
        end
        if p_values[i] < α / n_env
            break   # early termination
        end
    end
    # Bonferroni correction
    p_value = min(minimum(p_values) * n_env, 1.)
    if p_value < α
        # reject
        reject = true
        conf_intervals = zeros(length(S), 2)  # min = max = 0.
    else
        # not reject
        reject = false
        # the pooled fit is accepted
        if length(S) == 0
            conf_intervals = zeros(0, 2)
        elseif add_intercept
            conf_intervals = confint(fit0)[2:end, :]  # intercept is always the 1st row
        else
            conf_intervals = confint(fit0)
        end
    end
    return reject, p_value, conf_intervals
end

function conditional_inv_test_logistic(df::DataFrame, target::Symbol, S::Vector{Int64},
    env::Vector{Int64}, n_env::Int64; α=0.01, add_intercept=true, method="logistic-LR")
    conditional_inv_test_logistic(df, target, propertynames(df)[S], env, n_env;
        α=α, add_intercept=add_intercept, method=method)
end

"""
    get_formula(df::DataFrame, target::Symbol, S=setdiff(names(df), [target]); add_intercept=true)

Generate formula of `y ~ .` type for `df`, where `y` is specified by `target`.
"""
function get_formula(df::DataFrame, target::Symbol, S=setdiff(names(df), [target]); add_intercept=true)
    if length(S)==0
        return @eval @formula($target ~ 1)
    end
    if add_intercept
        return @eval @formula($target ~ +(1, $(S...)))
    else
        return @eval @formula($target ~ +($(S...), 0))
    end
end

"""
    sukhatme_fisher_test(x, y)

Sukhatme-Fisher test of H0: vectors x and y are two independent normal samples with equal mean and variance

See Perng, S. K., and Ramon C. Littell. "A test of equality of two normal population means and variances."
    Journal of the American Statistical Association 71.356 (1976): 968-971.

Return: p-value
"""
function sukhatme_fisher_test(x::Vector{Float64}, y::Vector{Float64})
    m = length(x)
    n = length(y)
    ss1 = var(x) * (m - 1)
    ss2 = var(y) * (n - 1)
    T = (mean(y) - mean(x)) / sqrt((m + n) * (ss1 + ss2) / ((m + n - 2) * (m * n)))
    F = (ss2 * (m-1)) / (ss1 * (n-1))
    g = FDist(n-1, m-1)
    if F >= median(g)
        H = 2 * (1 - cdf(g, F))
    else
        H = 2 * cdf(g, F)
    end
    W = - 2 * log(H)
    t = TDist(m + n - 2)
    w = Chisq(2)
    Q = - 2 * log(2 * (1 - cdf(t, abs(T)))) - 2 * log(1 - cdf(w, W))
    q = Chisq(4)
    return 1 - cdf(q, Q)
end

max_issues_repo_name: richardkwo/InvariantCausal

max_forks_repo_name: richardkwo/InvariantCausal


text:
module DiffEqSensitivity

using DiffEqBase, ForwardDiff, Tracker, FiniteDiff, Statistics
using DiffEqCallbacks, QuadGK, RecursiveArrayTools, LinearAlgebra
using Adapt
using Parameters: @unpack
using Requires
using StochasticDiffEq
using SharedArrays
import DiffEqNoiseProcess
import RandomNumbers: Xorshifts
using Random
import ZygoteRules, Zygote, ReverseDiff
import ArrayInterface
import Enzyme

using Cassette, DiffRules
using Core: CodeInfo, SlotNumber, SSAValue, ReturnNode, GotoIfNot

using Reexport
import ChainRulesCore: @thunk, NoTangent, @not_implemented
@reexport using GlobalSensitivity
abstract type SensitivityFunction end
abstract type TransformedFunction end

include("require.jl")
include("hasbranching.jl")
include("sensitivity_algorithms.jl")
include("derivative_wrappers.jl")
include("sensitivity_interface.jl")
include("forward_sensitivity.jl")
include("adjoint_common.jl")
include("lss.jl")
include("nilss.jl")
include("nilsas.jl")
include("backsolve_adjoint.jl")
include("interpolating_adjoint.jl")
include("quadrature_adjoint.jl")
include("callback_tracking.jl")
include("concrete_solve.jl")
include("second_order.jl")
include("steadystate_adjoint.jl")
include("sde_tools.jl")

export extract_local_sensitivities

export ODEForwardSensitivityFunction, ODEForwardSensitivityProblem, SensitivityFunction,
       ODEAdjointSensitivityProblem, ODEAdjointProblem, AdjointSensitivityIntegrand,
       SDEAdjointProblem, RODEAdjointProblem, SensitivityAlg,
       adjoint_sensitivities, adjoint_sensitivities_u0,
       ForwardLSSProblem, AdjointLSSProblem,
       NILSSProblem, NILSASProblem,
       shadow_forward, shadow_adjoint

export BacksolveAdjoint, QuadratureAdjoint, InterpolatingAdjoint,
       TrackerAdjoint, ZygoteAdjoint, ReverseDiffAdjoint,
       ForwardSensitivity, ForwardDiffSensitivity,
       ForwardDiffOverAdjoint,
       SteadyStateAdjoint,
       ForwardLSS, AdjointLSS, NILSS, NILSAS

export second_order_sensitivities, second_order_sensitivity_product

export TrackerVJP, ZygoteVJP, EnzymeVJP, ReverseDiffVJP

export StochasticTransformedFunction

function gsa(f, method::GlobalSensitivity.GSAMethod, args...; kwargs...)
    @warn "Global Sensitivity Analysis has been moved to separate package GlobalSensitivity.jl (https://github.com/SciML/GlobalSensitivity.jl)"
    GlobalSensitivity.gsa(f, method, args...; kwargs...)
end
end # module


text:
using DataArrays,DataFrames,StatsBase

function main()

  path = "../ml-100k/u.data"
  df = readrel(path,'\t')

  indices = shuffle([1:size(df,1)])

  # number of unique users
  totalusers = size(aggregate(df,:user_id,length),1)
  # total number of ratings
  totalratings = length(indices)

  # training

  trainingsetlimit = int(totalratings*(8/10))
  trainingindices  = indices[1:trainingsetlimit]
  trainingset      = df[trainingindices,:]

  # aggregating by ITEMS
  avgratingsbyitem  = aggregate(trainingset,:item_id,mean)[ [:item_id,:rating_mean] ]
  # and by USERS
  avgratingsbyuser  = aggregate(trainingset,:user_id,mean)[ [:user_id,:rating_mean] ]

  # global average for all items
  globalavgbyitem = mean(avgratingsbyitem[:rating_mean])

  # global average for all users
  globalavgbyuser = mean(avgratingsbyuser[:rating_mean])

  # test

  testsetlimit   = trainingsetlimit + (totalratings - trainingsetlimit)
  testindices    = indices[trainingsetlimit:testsetlimit]
  testset        = df[testindices,:]

  global sumerror = 0.0

  testratings = testset[:,[:user_id,:item_id,:rating]]

  totaltestinstances = size(testratings,1)

  for row in eachrow(testratings)

    # guess based upon the item itself
    item_id = row[:item_id]

    if item_id in avgratingsbyitem[:item_id]
      guessbyitem = avgratingsbyitem[:rating_mean][avgratingsbyitem[:item_id] .== item_id]
    else
      guessbyitem = globalavgbyitem
    end

    # guess based upon the user
    user_id = row[:user_id]

    if user_id in avgratingsbyuser[:user_id]
      guessbyuser = avgratingsbyuser[:rating_mean][avgratingsbyuser[:user_id] .== user_id]
    else
      guessbyuser = globalavgbyuser
    end

    actual = row[:rating]

    guess = (guessbyitem + guessbyuser)/2

    error = abs(actual-guess)
    sumerror += error

#     print("user_id: ")
#     print(user_id)
#     print(" guess: ")
#     print(guess)
#     print(" actual: ")
#     print(actual)
#     print(" error: ")
#     println(abs(actual-guess))

  end

  # dividing the error by the number of instances in the test set
  mae = sumerror/totaltestinstances

  mae

end

function readrel(relpath,delimiter=',')

  base = dirname(@__FILE__)

  if relpath[1] == '/'
    relpathwslash = relpath
  else
    relpathwslash = string("/",relpath)
  end

  path = string(base,relpathwslash)

  df = readtable(path,separator = delimiter)

  names!(df,[:user_id,:item_id,:rating,:timestamp])

  df

end

mae = main()
println(mae)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dash_graph2d

"""
    dash_graph2d(;kwargs...)

A Graph2D component.

Keyword arguments:
- `id` (String; required): The ID used to identify this component in Dash callbacks.
- `active` (Bool; optional): whether or not session is active. Used to enable or disable warning browser dialog when closing
- `backgroundColor` (String; optional): Getter/setter for the canvas background color, default transparent
- `centerAtZoom` (Dict with Strings as keys and values of type Real; optional): calls centerAt, then zoom. Takes an object with keys "k", "x", "y"
- `controlType` (String; optional): Which type of control to use to control the camera on 3D mode. Choice between trackball, orbit or fly.
- `cooldownTicks` (Real; optional): How many build-in frames to render before stopping and freezing the layout engine.
- `cooldownTime` (Real; optional): How long (ms) to render for before stopping and freezing the layout engine.
- `currentZoomPan` (Dict; optional)
- `emitParticle` (Dict; optional): An alternative mechanism for generating particles, this method emits a non-cyclical single particle within a specific link. The emitted particle shares the styling (speed, width, color) of the regular particle props. A valid link object that is included in graphData should be passed as a single parameter.
- `enableNavigationControls` (Bool; optional): Whether to enable the trackball navigation controls used to move the camera using mouse interactions (rotate/zoom/pan).
- `enableNodeDrag` (Bool; optional): Whether to enable the user interaction to drag nodes by click-dragging. If enabled, every time a node is dragged the simulation is re-heated so the other nodes react to the changes. Only applicable if enablePointerInteraction is true.
- `enablePointerInteraction` (Bool; optional): Whether to enable the mouse tracking events. This activates an internal tracker of the canvas mouse position and enables the functionality of object hover/click and tooltip labels, at the cost of performance. If you're looking for maximum gain in your graph performance it's recommended to switch off this property.
- `enableZoomPanInteraction` (Bool; optional): Whether to enable zooming and panning user interactions on a 2D canvas.
- `externalobjectInput` (Real | String | Bool | Array | Dict; optional): externalobject_source:
- `externalobject_source` (String; optional): externalobject_source:
- `fixNodes` (Bool; optional): Whether to fix node coordinates after simulation has cooled
- `forceEngine` (String; optional): Which force-simulation engine to use (d3 or ngraph).
- `forceRefresh` (Real; optional)
- `getGraphBbox` (Bool; optional): Returns the current bounding box of the nodes in the graph, formatted as { x: [<num>, <num>], y: [<num>, <num>], z: [<num>, <num>] }.
If no nodes are found, returns null. Accepts an optional argument to define a custom node filter: node => <boolean>, which should return a truthy value if the node is to be included.
This can be useful to calculate the bounding box of a portion of the graph.
Bounding box is saved as the graphBbox prop
- `graphData` (Dict; required): Graph data structure. Prop which is provided by user.
Can also be used to apply incremental updates. Format {nodes:{}, links:{}}
- `heightRatio` (Real; optional): height of component as proportion of container
- `key` (String; optional): The key used to identify this component in React
- `linkAutoColor` (Bool; optional): Automatically color link with inverse of background color
- `linkAutoColorBy` (String; optional): Link object accessor function or attribute to automatically group colors by. Only affects links without a color attribute.
- `linkClicked` (Dict; optional): clicked link
- `linkColor` (String; optional): String giving link color
- `linkCurvature` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the curvature radius of the link line. A value of 0 renders a straight line. 1 indicates a radius equal to half of the line length, causing the curve to approximate a semi-circle. For self-referencing links (source equal to target) the curve is represented as a loop around the node, with length proportional to the curvature value.
- `linkDesc` (String; optional): For VR only. Link object accessor function or attribute for description (shown under label).
- `linkDirectionalArrowColor` (String; optional): Link object accessor function or attribute for the color of the arrow head.
- `linkDirectionalArrowLength` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the length of the arrow head indicating the link directionality. The arrow is displayed directly over the link line, and points in the direction of source > target. A value of 0 hides the arrow.
- `linkDirectionalArrowRelPos` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the longitudinal position of the arrow head along the link line, expressed as a ratio between 0 and 1, where 0 indicates immediately next to the source node, 1 next to the target node, and 0.5 right in the middle.
- `linkDirectionalArrowResolution` (Real; optional): Geometric resolution of the arrow head, expressed in how many slice segments to divide the cone base circumference. Higher values yield smoother arrows.
- `linkDirectionalParticleColor` (String; optional): Link object accessor function or attribute for the directional particles color.
- `linkDirectionalParticleResolution` (Real; optional): Geometric resolution of each 3D directional particle, expressed in how many slice segments to divide the circumference. Higher values yield smoother particles.
- `linkDirectionalParticleSpeed` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the directional particles speed, expressed as the ratio of the link length to travel per frame. Values above 0.5 are discouraged.
- `linkDirectionalParticleWidth` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the directional particles width. Values are rounded to the nearest decimal for indexing purposes.
- `linkDirectionalParticles` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the number of particles (small spheres) to display over the link line. The particles are distributed equi-spaced along the line, travel in the direction source > target, and can be used to indicate link directionality.
- `linkHoverPrecision` (Real; optional): Whether to display the link label when gazing the link closely (low value) or from far away (high value).
- `linkId` (String; optional): The link attribute containing the unique link id
- `linkIdsHighlightUser` (Array of Strings; optional): ids of highlighted links (through search)
- `linkIdsInvisibleAuto` (Array of Strings; optional)
- `linkIdsInvisibleUser` (Array of Strings; optional): ids of visible links
- `linkLabel` (String; optional): Link object accessor function or attribute for name (shown in label). Supports plain text or HTML content (except in VR).
- `linkOpacity` (Real; optional): Line opacity of links, between [0,1]. 3D, VR, AR
- `linkResolution` (Real; optional): Geometric resolution of each link 3D line, expressed in how many radial segments to divide the cylinder. Higher values yield smoother cylinders. Applicable only to links with positive width. 3D, VR, AR
- `linkRightClicked` (Dict; optional): right-clicked link
- `linkRightClickedViewpointCoordinates` (Dict with Strings as keys and values of type Real; optional)
- `linkSource` (String; optional): Link object accessor attribute referring to id of source node.
- `linkTarget` (String; optional): Link object accessor attribute referring to id of target node.
- `linkWidth` (Real | String; optional): Link object accessor function, attribute or a numeric constant for the link line width.
- `linksSelected` (Array of Dicts; optional): selected (clicked) links
- `maxDepth_neighbours_select` (Real; optional): max levels of neighbourhood selection around a node by repeat clicking
- `maxZoom` (Real; optional)
- `minZoom` (Real; optional): ids of links highlighted due to being dragged
- `n_linkRightClicks` (Real; optional)
- `n_nodeRightClicks` (Real; optional)
- `ngraphPhysics` (Dict; optional): Specify custom physics configuration for ngraph, according to its configuration object syntax. Only applicable if using the ngraph simulation engine.
- `nodeAutoColorBy` (String; optional): Node object accessor function or attribute to automatically group colors by. Only affects nodes without a color attribute.
- `nodeClicked` (Dict; optional): left-clicked node
- `nodeColor` (String; optional): The node attribute whose value should be used for coloring nodes
- `nodeCoordinates` (String | Array of Reals; optional): Ratio of node circle area (square px) [2D] or sphere volume (cubic px) [3D] per value unit.
- `nodeDesc` (String; optional): For VR only. Node object accessor function or attribute for description (shown under label).
- `nodeIcon` (String; optional): The node attribute containing object with icon to display for each individual node.
- `nodeIconSizeFactor` (Real; optional): controls nodeIcon size, relative to nodeRelSize
- `nodeId` (String; optional): Node object accessor attribute for unique node id (used in link objects source/target).
- `nodeIdsHighlightUser` (Array of Strings; optional): ids of highlighted nodes (through search)
- `nodeIdsInvisibleAuto` (Array of Strings; optional): ids of visible nodes. Not to be supplied by user. Available to allow for saving state
- `nodeIdsInvisibleUser` (Array of Strings; optional): ids of invisible nodes supplied by user as prop
- `nodeImg` (String; optional): The node attribute containing url to image to display for each individual node
- `nodeImgSizeFactor` (Real; optional): controls nodeImg size, reltive to nodeRelSize
- `nodeLabel` (String; optional): Node object accessor function or attribute for name (shown in label). Supports plain text or HTML content (except in VR).
2D, 3D and VR
- `nodeLabelRelSize` (Real; optional): controls node label size
- `nodeOpacity` (Real; optional): Nodes sphere opacity, between [0,1]. 3D, VR, AR
- `nodeRelSize` (Real; optional): Ratio of node circle area (square px) [2D] or sphere volume (cubic px) [3D] per value unit.
- `nodeResolution` (Real; optional): Geometric resolution of each node's sphere, expressed in how many slice segments to divide the circumference. Higher values yield smoother spheres. Only applicable to 3D modes.
3D, VR, AR
- `nodeRightClicked` (Dict; optional): right-clicked node
- `nodeRightClickedViewpointCoordinates` (Dict with Strings as keys and values of type Real; optional): screen coordinates of right-clicked node
- `nodeTextAutoColor` (Bool; optional): Automatically color node text with inverse of backgroundColor
- `nodeTextColor` (String; optional): Node text color
- `nodeVal` (Real | String; optional): Ratio of node circle area (square px) [2D] or sphere volume (cubic px) [3D] per value unit.
- `nodeZoomId` (String; optional): id of node to zoom to
- `nodesSelected` (Array of Dicts; optional): selected (clicked) nodes
- `numDimensions` (Real; optional): Not applicable to 2D mode. Number of dimensions to run the force simulation on. 3D, VR, AR
- `pauseAnimation` (Bool; optional): Pauses the rendering cycle of the component, effectively freezing the current view and cancelling all user interaction. This method can be used to save performance in circumstances when a static image is sufficient.
- `pixelUnitRatio` (Real; optional): pixelUnitRatio: if node attribute (in some unit of measurement) is used as coordinates, pixel:unit scale
- `refresh` (Bool; optional): Redraws all the nodes/links. 3D, VR, AR
- `rendererConfig` (Dict; optional): Configuration parameters to pass to the ThreeJS WebGLRenderer constructor. This prop only has an effect on component mount. 3D only
- `scripts` (Array of Strings; optional)
- `showCoordinates` (Bool; optional): showCoordinates: whether or not to show pointer coordinates in hover tooltip (not yet used)
- `showGUI` (Bool; optional)
- `size` (Dict; optional): provided react-sizeme. Contains an object with "width" and "height" attributes
- `sortRels1Descend` (Bool; optional): sort in descending order?
- `sortRels2Descend` (Bool; optional): sort in descending order?
- `sortRelsBy1` (String; optional): in zoom view, node attribute to sort relations by first
- `sortRelsBy2` (String; optional): in zoom view, node attribute to sort relations by after first sort
- `sortRoleplayers1Descend` (Bool; optional): sort in descending order?
- `sortRoleplayers2Descend` (Bool; optional): sort in descending order?
- `sortRoleplayersBy1` (String; optional): in zoom view, node attribute to sort role players by first
- `sortRoleplayersBy2` (String; optional): in zoom view, node attribute to sort role players by after first sort
- `updateNeighbours` (Bool; optional)
- `useCoordinates` (Bool; optional): useCoordinates: whether to set nodeCoordinates for node coordinates
- `warmupTicks` (Real; optional): Number of layout engine cycles to dry-run at ignition before starting to render.
- `zoom` (Array of Reals; optional): Calls zoom() method. ([number], [ms])
- `zoomToFit` (Real; optional)
"""
function dash_graph2d(; kwargs...)
        available_props = Symbol[:id, :active, :backgroundColor, :centerAtZoom, :controlType, :cooldownTicks, :cooldownTime, :currentZoomPan, :emitParticle, :enableNavigationControls, :enableNodeDrag, :enablePointerInteraction, :enableZoomPanInteraction, :externalobjectInput, :externalobject_source, :fixNodes, :forceEngine, :forceRefresh, :getGraphBbox, :graphData, :heightRatio, :key, :linkAutoColor, :linkAutoColorBy, :linkClicked, :linkColor, :linkCurvature, :linkDesc, :linkDirectionalArrowColor, :linkDirectionalArrowLength, :linkDirectionalArrowRelPos, :linkDirectionalArrowResolution, :linkDirectionalParticleColor, :linkDirectionalParticleResolution, :linkDirectionalParticleSpeed, :linkDirectionalParticleWidth, :linkDirectionalParticles, :linkHoverPrecision, :linkId, :linkIdsHighlightUser, :linkIdsInvisibleAuto, :linkIdsInvisibleUser, :linkLabel, :linkOpacity, :linkResolution, :linkRightClicked, :linkRightClickedViewpointCoordinates, :linkSource, :linkTarget, :linkWidth, :linksSelected, :maxDepth_neighbours_select, :maxZoom, :minZoom, :n_linkRightClicks, :n_nodeRightClicks, :ngraphPhysics, :nodeAutoColorBy, :nodeClicked, :nodeColor, :nodeCoordinates, :nodeDesc, :nodeIcon, :nodeIconSizeFactor, :nodeId, :nodeIdsHighlightUser, :nodeIdsInvisibleAuto, :nodeIdsInvisibleUser, :nodeImg, :nodeImgSizeFactor, :nodeLabel, :nodeLabelRelSize, :nodeOpacity, :nodeRelSize, :nodeResolution, :nodeRightClicked, :nodeRightClickedViewpointCoordinates, :nodeTextAutoColor, :nodeTextColor, :nodeVal, :nodeZoomId, :nodesSelected, :numDimensions, :pauseAnimation, :pixelUnitRatio, :refresh, :rendererConfig, :scripts, :showCoordinates, :showGUI, :size, :sortRels1Descend, :sortRels2Descend, :sortRelsBy1, :sortRelsBy2, :sortRoleplayers1Descend, :sortRoleplayers2Descend, :sortRoleplayersBy1, :sortRoleplayersBy2, :updateNeighbours, :useCoordinates, :warmupTicks, :zoom, :zoomToFit]
        wild_props = Symbol[]
        return Component("dash_graph2d", "Graph2D", "dash_react_force_graph", available_props, wild_props; kwargs...)
end

max_stars_repo_name: nodelab-aps/dash-react-force-graph

max_issues_repo_name: nodelab-aps/dash-react-force-graph

max_forks_repo_name: nodelab-aps/dash-react-force-graph


text:
module Master

export Instance_utils, Life, Startup_utils, Types

include("types.jl")
using Master.Types

include("instance_utils.jl")
using Master.Instance_info

include("startup_utils.jl")
using Master.Startup_utils

include("life.jl")
using Master.Life

end # module

max_stars_repo_name: Gollor/JuliaAutomation

max_issues_repo_name: Gollor/JuliaAutomation

max_forks_repo_name: Gollor/JuliaAutomation


text:
using CSFML.LibCSFML

mode = sfVideoMode(1280, 720, 32)

window = sfRenderWindow_create(mode, "SFML window", sfResize | sfClose, C_NULL)
@assert window != C_NULL

texture = sfTexture_createFromFile(joinpath(@__DIR__, "julia-tan.png"), C_NULL)
@assert texture != C_NULL

sprite = sfSprite_create()
sfSprite_setTexture(sprite, texture, sfTrue)

font = sfFont_createFromFile(joinpath(@__DIR__, "Roboto-Bold.ttf"))
@assert font != C_NULL

text = sfText_create()
sfText_setString(text, "Hello SFML")
sfText_setFont(text, font)
sfText_setCharacterSize(text, 50)

music = sfMusic_createFromFile(joinpath(@__DIR__, "Chrono_Trigger.ogg"))
@assert music != C_NULL

sfMusic_play(music)

event_ref = Ref{sfEvent}()

while Bool(sfRenderWindow_isOpen(window))
    # process events
    while Bool(sfRenderWindow_pollEvent(window, event_ref))
        # close window : exit
        event_ref[].type == sfEvtClosed && sfRenderWindow_close(window)
        event_ref[].type == sfEvtResized && println("Trigger sfEvtResized.")
        event_ref[].type == sfEvtLostFocus && println("Trigger sfEvtLostFocus.")
        event_ref[].type == sfEvtGainedFocus && println("Trigger sfEvtGainedFocus.")
        event_ref[].type == sfEvtTextEntered && println("Trigger sfEvtTextEntered: $(event_ref[].text.unicode)")
        event_ref[].type == sfEvtKeyPressed && println("Trigger sfEvtKeyPressed: $(event_ref[].key.code)")
        event_ref[].type == sfEvtKeyReleased && println("Trigger sfEvtKeyReleased: $(event_ref[].key.code)")
        event_ref[].type == sfEvtMouseWheelMoved && println("Trigger sfEvtMouseWheelMoved: $(event_ref[].mouseWheel.sfMouseVerticalWheel)")
        event_ref[].type == sfEvtMouseWheelScrolled && println("Trigger sfEvtMouseWheelScrolled: $(event_ref[].mouseWheel.sfMouseVerticalWheel)")
        event_ref[].type == sfEvtMouseButtonPressed && println("Trigger sfEvtMouseButtonPressed: $(event_ref[].mouseButton.button)")
        event_ref[].type == sfEvtMouseButtonReleased && println("Trigger sfEvtMouseButtonReleased: $(event_ref[].mouseButton.x), $(event_ref[].mouseButton.y)")
        event_ref[].type == sfEvtMouseMoved && println("Trigger sfEvtMouseMoved: $(event_ref[].mouseMove.x), $(event_ref[].mouseMove.y)")
    end
    # clear the screen
    sfRenderWindow_clear(window, sfColor_fromRGBA(0,0,0,1))
    # draw the sprite
    sfRenderWindow_drawSprite(window, sprite, C_NULL)
    # draw the text
    sfRenderWindow_drawText(window, text, C_NULL)
    # update the window
    sfRenderWindow_display(window)
end

sfMusic_destroy(music)
sfText_destroy(text)
sfFont_destroy(font)
sfSprite_destroy(sprite)
sfTexture_destroy(texture)
sfRenderWindow_destroy(window)


text:
# # Simulations with TAMU data and [PowerSimulations.jl](https://github.com/NREL/PowerSimulations.jl)

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# This is a basic simulation example using the [TAMU Cases](https://electricgrids.engr.tamu.edu/).

# ## Dependencies
using SIIPExamples
pkgpath = dirname(dirname(pathof(SIIPExamples)))
using Dates

# ### Modeling Packages
using PowerSystems
using PowerSimulations
using PowerSystemCaseBuilder

# ### Optimization packages
# For this simple example, we can use the Cbc solver with a relatively relaxed tolerance.
using Cbc #solver
solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.5)

# ### Create a `System` from TAMU data
# We have included some of the TAMU cases (with truncated time series info)
# in the PowerSystemCaseBuilder testing, so we can just use that.
sys = build_system(PSYTestSystems, "tamu_ACTIVSg2000_sys")

horizon = 24;
interval = Dates.Hour(24);
transform_single_time_series!(sys, horizon, interval);

# ## Run a PCM
# note that the TAMU data doesn't contain startup and shutdown costs, or minimum up/down
# time limits, so a UC problem merely respects minmum generation levels.

sim_folder = mkpath(joinpath(pkgpath, "TAMU-sim"))
problems = SimulationProblems(
    UC = OperationsProblem(
        template_unit_commitment(transmission = CopperPlatePowerModel),
        sys,
        optimizer = solver,
    ),
)
intervals = Dict("UC" => (Hour(24), Consecutive()))
DA_sequence = SimulationSequence(
    problems = problems,
    intervals = intervals,
    ini_cond_chronology = IntraProblemChronology(),
)

# ### Define and build a simulation
sim = Simulation(
    name = "TAMU-test",
    steps = 3,
    problems = problems,
    sequence = DA_sequence,
    simulation_folder = sim_folder,
)

build!(sim)

# ### Execute the simulation
execute!(sim)

# ### Load and analyze results
sim_results = SimulationResults(sim);
uc_results = get_problem_results(sim_results, "UC")
read_realized_variables(uc_results, names = [:On__ThermalStandard])[:On__ThermalStandard]


text:
using LinearAlgebra: norm

export RandomSearch, NullSearch

"""
    Search(p::DSProblem{T})::IterationOutcome  where T

Generate search points and call evaluate on them.
"""
function Search(p::DSProblem{T})::IterationOutcome  where T
    p.status.search_time_total += @elapsed points = GenerateSearchPoints(p, p.config.search)
    p.full_output && OutputSearchStep(p, points)
    p.status.directions = nothing
    return EvaluatePoint!(p, points)
end

"""
    NullSearch()

Return no trial points for a search stage (ie, skips the
search stage from running)
"""
struct NullSearch <: AbstractSearch end

"""
    GenerateSearchPoints(p::DSProblem)

Search method that returns an empty vector.

Use when no search method is desired.
"""
(GenerateSearchPoints(p::DSProblem{T}, ::NullSearch)::Vector{Vector{T}}) where T = Vector{T}[]

"""
    GenerateSearchPoints(p::DSProblem{T})::Vector{Vector{T}} where T

Calls `GenerateSearchPoints` for the search step within `p`.
"""
(GenerateSearchPoints(p::DSProblem{T})::Vector{Vector{T}}) where T = GenerateSearchPoints(p, p.config.search)

"""
    RandomSearch(M::Int)

Return `M` randomly chosen trial points from the current mesh.
"""
mutable struct RandomSearch <: AbstractSearch
    M::Int #Number of points to generate
end

"""
    GenerateSearchPoints(p::DSProblem{T}, ::RandomSearch)::Vector{Vector{T}} where T

Finds points that are Δᵐ distance from any point in the mesh in a uniformly random direction.
"""
(GenerateSearchPoints(p::DSProblem{T}, s::RandomSearch)::Vector{Vector{T}}) where T =
                    RandomPointsFromCache(p.N, p.cache, p.config.mesh.δ, s)

function RandomPointsFromCache(N::Int, c::PointCache{T}, dist::Vector{T}, s::RandomSearch
                              )::Vector{Vector{T}} where T
    mesh_points = CacheRandomSample(c, s.M)

    if length(mesh_points) == s.M
        for i in 1:s.M
            dir = rand(N)
            dir .*= dist./norm(dir)
            mesh_points[i] += dir
        end
    end

    return mesh_points
end


text:
import Base: +,-

typealias clockid_t Int32
@enum(CLOCK_ID,
  CLOCK_REALTIME = 0,
  CLOCK_MONOTONIC = 1
)

const BILLION = 1000000000

if Base.Sys.WORD_SIZE==32
  typealias time_t Clong
elseif Base.Sys.WORD_SIZE==64
  typealias time_t Clonglong
else
  error("Unsupported architecture")
end

type timespec
  sec::time_t
  nsec::Clong
end

@inline function normalize!(t::timespec)
  if t.nsec >= BILLION
    (ss,ns)=divrem(t.nsec,BILLION)
    t.sec += ss
    t.nsec = ns
  end
  return t
end

function -(t1::timespec,t2::timespec)
  if t1.sec < 0 || t1.nsec < 0 || t2.sec < 0 || t2.nsec < 0
    error("Subtracting timespec is only defined for positive time values")
  end
  # find the larger time value
  if t1.sec-t2.sec > 0
    tend = t1; tbegin = t2; dsign=+1
  elseif t1.sec-t2.sec < 0
    tend = t2; tbegin = t1; dsign=-1
  else # t1.sec == t2.sec
    if t1.nsec > t2.nsec
      tend = t1; tbegin = t2; dsign=+1
    elseif t1.nsec < t2.nsec
      tend = t2; tbegin = t1; dsign=-1
    else # a tie!
      tend = t1; tbegin = t2; dsign=+1
    end
  end
  dnsec = tend.nsec - tbegin.nsec; dsec = tend.sec - tbegin.sec
  dnsec < 0 && (dnsec+=BILLION; dsec-=1)
  d::Int64 = dsign*(dnsec+BILLION*dsec)
end

module TIMER_FLAG
  const RELTIME = Val{0}
  const ABSTIME = Val{1}
end

function gettime!(tx::timespec,clockid::CLOCK_ID)
  s = ccall((:clock_gettime,librt),Int32,(clockid_t,Ptr{timespec}),
    clockid,pointer_from_objref(tx))
  s!=0 && error("Error in gettime()")
  return tx
end

function nanosleep(clockid::CLOCK_ID, tx::timespec, ::Type{TIMER_FLAG.ABSTIME})
  f = pointer_from_objref(tx) # hack, to avoid unnecessary memory allocation
  s = ccall((:clock_nanosleep,librt),Int32,(clockid_t,Int32,Ptr{timespec},Ptr{timespec}),
    clockid,1,f,f)
  s!=0 && error("Error in nanosleep()")
  return nothing
end
nanosleep(clockid::CLOCK_ID, t::timespec) = nanosleep(clockid, t, TIMER_FLAG.ABSTIME)

function nanosleep(clockid::CLOCK_ID, t::timespec, ::Type{TIMER_FLAG.RELTIME})
  error("Relative time nanosleep not supported yet!")
  return nothing
end

@inline function nanosleep!(tx::timespec,nanosec::Clong,clockid::CLOCK_ID)
  tx.nsec += nanosec
  normalize!(tx)
  nanosleep(clockid,tx)
  return tx
end
nanosleep!(tx::timespec,nanosec::Clong) = nanosleep!(tx,nanosec,CLOCK_MONOTONIC)


text:
function dump_state_and_aux_init(dgngrp, currtime)
    if isnothing(dgngrp.interpol)
        @warn """
            Diagnostics ($dgngrp.name): requires interpolation!
            """
    end
end

function dump_state_and_aux_collect(dgngrp, currtime)
    interpol = dgngrp.interpol
    if isnothing(interpol)
        @warn """
            Diagnostics ($dgngrp.name): requires interpolation!
            """
        return nothing
    end

    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q
    FT = eltype(Q.data)
    bl = dg.balance_law
    mpirank = MPI.Comm_rank(mpicomm)

    # filename (may also want to take out)
    dprefix = @sprintf(
        "%s_%s-%s-num%04d",
        dgngrp.out_prefix,
        dgngrp.name,
        Settings.starttime,
        dgngrp.num
    )
    statefilename = joinpath(Settings.output_dir, dprefix)
    auxfilename = joinpath(Settings.output_dir, "$(dprefix)_aux")

    statenames = flattenednames(vars_state_conservative(bl, FT))
    auxnames = flattenednames(vars_state_auxiliary(bl, FT))

    istate = similar(Q.data, interpol.Npl, number_state_conservative(bl, FT))
    interpolate_local!(interpol, Q.data, istate)

    if interpol isa InterpolationCubedSphere
        # TODO: get indices here without hard-coding them
        _ρu, _ρv, _ρw = 2, 3, 4
        project_cubed_sphere!(interpol, istate, (_ρu, _ρv, _ρw))
    end

    all_state_data = accumulate_interpolated_data(mpicomm, interpol, istate)

    iaux = similar(
        dg.state_auxiliary.data,
        interpol.Npl,
        number_state_auxiliary(bl, FT),
    )
    interpolate_local!(interpol, dg.state_auxiliary.data, iaux)

    all_aux_data = accumulate_interpolated_data(mpicomm, interpol, iaux)

    if mpirank == 0
        dims = dimensions(interpol)
        dim_names = tuple(collect(keys(dims))...)

        statevarvals = OrderedDict()
        for i in 1:number_state_conservative(bl, FT)
            statevarvals[statenames[i]] =
                (dim_names, all_state_data[:, :, :, i], Dict())
        end
        write_data(dgngrp.writer, statefilename, dims, statevarvals, currtime)

        auxvarvals = OrderedDict()
        for i in 1:number_state_auxiliary(bl, FT)
            auxvarvals[auxnames[i]] =
                (dim_names, all_aux_data[:, :, :, i], Dict())
        end
        write_data(dgngrp.writer, auxfilename, dims, auxvarvals, currtime)
    end

    MPI.Barrier(mpicomm)
    return nothing
end

function dump_state_and_aux_fini(dgngrp, currtime) end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
let
    T = 50
    N = 50

    pomdp = BabyPOMDP()

    rsum = 0.0
    policy = solve(BlindPolicySolver(verbose=true), pomdp)
    for i in 1:N
        sim = RolloutSimulator(max_steps=T, rng=MersenneTwister(i))
        rsum += simulate(sim, pomdp, policy)
    end

    @show rsum/N
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(PyPlot.find_backend), PyCall.PyObject})
    precompile(Tuple{getfield(PyPlot, Symbol("##draw#40")), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, typeof(PyPlot.draw)})
    precompile(Tuple{getfield(PyPlot, Symbol("##figure#7")), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, typeof(PyPlot.figure)})
    precompile(Tuple{typeof(PyPlot.isdisplayok)})
    precompile(Tuple{getfield(PyPlot, Symbol("##ioff#68")), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, typeof(PyPlot.ioff)})
    precompile(Tuple{typeof(PyPlot.__init__)})
    precompile(Tuple{typeof(PyPlot.gcf)})
    precompile(Tuple{typeof(PyPlot.init_colormaps)})
    precompile(Tuple{typeof(PyPlot.figure)})
    precompile(Tuple{typeof(PyPlot.ioff)})
end

max_stars_repo_name: jdeldre/BookExamples

max_issues_repo_name: jdeldre/BookExamples

max_forks_repo_name: jdeldre/BookExamples


text:
module HSMath

    export powi, optimize, offset_func, offset_calc, cplx_floor, cplx_ceil, cplx_mod, ∏, ∑, ≆, ≊

    include("Primes.jl")

    powi(base, pow) = base < 0 ? -((-base) ^ pow) : base ^ pow

    """f -> Function in modifying g(x)
       g -> Internal Core Function
       o -> Operator function

       Returns required such that O(f(g(x))) = f(g(x) + C). Solves for C
    """

    ≆(x, p) = !isapprox(x, p)
    ≊(x, p) = isapprox(x, p)

    offset_func(f::Function, f_inv::Function, g::Function, operator::Function) = x -> offset_calc(x, f, f_inv, g, operator)

    offset_calc(x, f::Function, f_inv::Function, g::Function, operator::Function) = f_inv(operator(f(g(x)))) - g(x)

    cplx_floor(x) = Base.floor(real(x)) + Base.floor(imag(x)) * im

    cplx_ceil(x) = Base.ceil(real(x)) + Base.ceil(imag(x)) * im

    cplx_mod(base, num) = base - num * cplx_floor(base / num)

    function ∏(array::NTuple{N, T}) where {T, N}
        prod::T = 1
        for arg in array
            prod *= arg
        end
        prod
    end

    function ∏(array::AbstractArray{T}) where T
        prod::T = 1
        for arg in array
            prod *= arg
        end
        prod
    end

    function ∑(array::AbstractArray{T}) where T
        sum::T = 0
        for arg in array
            sum += arg
        end
        sum
    end

    function ∑(array::NTuple{N, T}) where {T, N}
        sum::T = 0
        for arg in array
            sum += arg
        end
        sum
    end
end


text:
CARDPK = 1526110
DOORPK = 20175123

SUBJECT = 7
SALT = 20201227

function find_loop_size(pk, subject)
    i = 0
    n = 1
    while n != pk
        i += 1
        n *= subject
        n = mod(n, SALT)
    end
    return i
end

function enc(subject, loops)
    n = 1
    for _ in 1:loops
        n *= subject
        n = mod(n, SALT)
    end
    return n
end

card = find_loop_size(CARDPK, SUBJECT)
door = find_loop_size(DOORPK, SUBJECT)
p1 = enc(DOORPK, card)

println("-----------------------------------------------------------------------")
println("combo breaker -- part one :: $p1")
println("-----------------------------------------------------------------------")

@assert(p1 == 10924063)


text:
using Completion: LPO

@testset "lpo" begin
    >ₗₚₒ = LPO((:i, 1), (:f, 2), (:e, 0))
    @test !(@term(x) >ₗₚₒ @term(y))
    @test @term(f(x, e())) >ₗₚₒ @term(x)
    @test @term(i(e())) >ₗₚₒ @term(e())
    @test @term(i(f(x, y))) >ₗₚₒ @term(f(i(y), i(x)))
    @test @term(f(f(x, y), z)) >ₗₚₒ @term(f(x, f(y, z)))
    @test @term(f(x, y)) >ₗₚₒ @term(x)

    @test_throws ArgumentError("(:f, 1) is not contained in order") LPO((:i, 1))(@term(f(x)), @term(i(x)))
    @test_throws ArgumentError("(:g, 1) is not contained in order") LPO((:i, 1))(@term(i(x)), @term(g(x)))
end


text:
mutable struct MixedSumAcc{T}
    s :: T
end

mixedSumAcc(::Type{Float32})                   = MixedSumAcc(Float64(0))
mixedSumAcc(::Type{Vec{N, Float32}}) where {N} = MixedSumAcc(vzero(Vec{N, Float64}))

function add!(acc::MixedSumAcc, x)
    @fastmath acc.s += f64(x)
end

function add!(acc::MixedSumAcc{T}, x::MixedSumAcc{T}) where {T}
    @fastmath acc.s += x.s
end

Base.sum(acc::MixedSumAcc{T}) where {T<:Vec}  = MixedSumAcc(vsum(acc.s))
Base.sum(acc::MixedSumAcc{T}) where {T<:Real} = acc.s


text:
#!/usr/bin/env -S julia -O3
# =======================================================================================
#
#      Author:   Thomas Gruber (tg), thomas.gruber@fau.de
#      Copyright (c) 2020 RRZE, University Erlangen-Nuremberg
#
#      Permission is hereby granted, free of charge, to any person obtaining a copy
#      of this software and associated documentation files (the "Software"), to deal
#      in the Software without restriction, including without limitation the rights
#      to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#      copies of the Software, and to permit persons to whom the Software is
#      furnished to do so, subject to the following conditions:
#
#      The above copyright notice and this permission notice shall be included in all
#      copies or substantial portions of the Software.
#
#      THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#      IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#      FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#      AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#      LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#      OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
#      SOFTWARE.
#
# =======================================================================================

using Printf

const NTIMES = Int64(5)
const VECTORSIZE = Int64(120000000)

struct Benchmark
    label::String
    words::Int
    flops::Int
end

@enum BenchmarkNum begin
    INIT = 1
    COPY = 2
    UPDATE = 3
    TRIAD = 4
    DAXPY = 5
    STRIAD = 6
    SDAXPY = 7
end

function init_func(A::Array{Float64,1}, scalar::Float64, vsize::Int64)
    @assert length(A) == vsize
    num_threads = Threads.nthreads()
    if num_threads == 1
        @simd for j in 1:vsize
            @inbounds A[j] = scalar
        end
    else
        Threads.@threads for j in 1:vsize
            @inbounds A[j] = scalar
        end
# Version according to https://slides.com/valentinchuravy/julia-parallelism#/5/6/1 with
# remainder loop but currently unused. Instead of work-sharing loop over the threads, we
# could also use Threads.@spawn, Threads.threadid() and wait().
#        len = div(vsize, num_threads)
#        remain = vsize - (len * num_threads)
#        Threads.@threads for t in 1:num_threads
#            st = (t-1)*len + 1
#            en = min(t*len, vsize)
#            dom = st:en
#            @inbounds @simd for j in dom
#                A[j] = scalar
#            end
#        end
#        if remain > 0
#            @inbounds @simd for j in (vsize-remain):vsize
#                A[j] = scalar
#            end
#        end
    end
end

function update_func(A::Array{Float64,1}, scalar::Float64, vsize::Int64)
    @assert length(A) == vsize
    if Threads.nthreads() == 1
        @simd for j in 1:vsize
            @inbounds A[j] = A[j] * scalar
        end
    else
        Threads.@threads for j in 1:vsize
            @inbounds A[j] = A[j] * scalar
        end
    end
end

function copy_func(A::Array{Float64,1}, B::Array{Float64,1}, vsize::Int64)
    @assert length(A) == length(B) == vsize
    if Threads.nthreads() == 1
        @simd for j in 1:vsize
            @inbounds A[j] = B[j]
        end
    else
        Threads.@threads for j in 1:vsize
            @inbounds A[j] = B[j]
        end
    end
end

function daxpy_func(A::Array{Float64,1}, B::Array{Float64,1}, scalar::Float64, vsize::Int64)
    @assert length(A) == length(B) == vsize
    if Threads.nthreads() == 1
        @simd for j in 1:vsize
            @inbounds A[j] = A[j] + scalar * B[j]
        end
    else
        Threads.@threads for j in 1:vsize
            @inbounds A[j] = A[j] + scalar * B[j]
        end
    end
end

function sdaxpy_func(A::Array{Float64,1}, B::Array{Float64,1}, C::Array{Float64,1}, vsize::Int64)
    @assert length(A) == length(B) == length(C) == vsize
    if Threads.nthreads() == 1
        @simd for j in 1:vsize
            @inbounds A[j] = A[j] + B[j] * C[j]
        end
    else
        Threads.@threads for j in 1:vsize
            @inbounds A[j] = A[j] + B[j] * C[j]
        end
    end
end

function triad_func(A::Array{Float64,1}, B::Array{Float64,1}, C::Array{Float64,1}, scalar::Float64, vsize::Int64)
    @assert length(A) == length(B) == length(C) == vsize
    if Threads.nthreads() == 1
        @simd for j in 1:vsize
            @inbounds A[j] = B[j] + scalar * C[j]
        end
    else
        Threads.@threads for j in 1:vsize
            @inbounds A[j] = B[j] + scalar * C[j]
        end
    end
end

function striad_func(A::Array{Float64,1}, B::Array{Float64,1}, C::Array{Float64,1}, D::Array{Float64,1}, vsize::Int64)
    @assert length(A) == length(B) == length(C) == length(D) == vsize
    if Threads.nthreads() == 1
        @simd for j in 1:vsize
            @inbounds A[j] = B[j] + C[j] * D[j]
        end
    else
        Threads.@threads for j in 1:vsize
            @inbounds A[j] = B[j] + C[j] * D[j]
        end
    end
end

function check(A::Array{Float64,1}, B::Array{Float64,1}, C::Array{Float64,1}, D::Array{Float64,1}, vsize::Int64)
    aj = 2.0
    bj = 2.0
    cj = 0.5
    dj = 1.0
    scalar = 3.0
    epsilon = 1.e-8;

    for k in 1:NTIMES
        bj = scalar
        cj = aj
        aj = aj * scalar
        aj = bj + scalar * cj
        aj = aj + scalar * bj
        aj = bj + cj * dj
        aj = aj + bj * cj
    end

    aj = aj * Float64(vsize)
    bj = bj * Float64(vsize)
    cj = cj * Float64(vsize)
    dj = dj * Float64(vsize)

    asum = sum(A)
    bsum = sum(B)
    csum = sum(C)
    dsum = sum(D)

    if abs(aj - asum) / asum > epsilon
        @printf("Failed Validation on array a[]\n\tExpected  : %f\n\tObserved  : %f\n", aj, asum)
    elseif abs(bj - bsum) / bsum > epsilon
        @printf("Failed Validation on array b[]\n\tExpected  : %f\n\tObserved  : %f\n", bj, bsum)
    elseif abs(cj - csum) / csum > epsilon
        @printf("Failed Validation on array c[]\n\tExpected  : %f\n\tObserved  : %f\n", cj, csum)
    elseif abs(dj - dsum) / dsum > epsilon
        @printf("Failed Validation on array d[]\n\tExpected  : %f\n\tObserved  : %f\n", dj, dsum)
    else
        @printf("Solution Validates\n")
    end
end

function bwBench()

    A = zeros(Float64, VECTORSIZE)
    B = zeros(Float64, VECTORSIZE)
    C = zeros(Float64, VECTORSIZE)
    D = zeros(Float64, VECTORSIZE)

    benchmarks = Vector{Benchmark}
    benchmarks = [
        Benchmark("Init", 1, 0),
        Benchmark("Copy", 2, 0),
        Benchmark("Update", 2, 1),
        Benchmark("Triad", 3, 2),
        Benchmark("Daxpy", 3, 2),
        Benchmark("STriad", 4, 2),
        Benchmark("SDaxpy", 2, 2),
    ]

    NUMBENCH = size(benchmarks, 1)

    avgtime = zeros(Float64, NUMBENCH)
    mintime = zeros(Float64, NUMBENCH)
    maxtime = zeros(Float64, NUMBENCH)
    for i in 1:NUMBENCH
        mintime[i] = typemax(Float64)
    end
    times = zeros(Float64, (NUMBENCH, NTIMES))

    scalar = 3.0

    if Threads.nthreads() > 1
        @printf("Multi-threading enabled, running with %d threads\n", Threads.nthreads());
    end

    for j in 1:VECTORSIZE
        A[j] = 2.0
        B[j] = 2.0
        C[j] = 0.5
        D[j] = 1.0
    end

    for j in 1:NTIMES
        times[Int(INIT), j] = @elapsed init_func(B, scalar, VECTORSIZE)
        times[Int(COPY), j] = @elapsed copy_func(C, A, VECTORSIZE)
        times[Int(UPDATE), j] = @elapsed update_func(A, scalar, VECTORSIZE)
        times[Int(TRIAD), j] = @elapsed triad_func(A, B, C, scalar, VECTORSIZE)
        times[Int(DAXPY), j] = @elapsed daxpy_func(A, B, scalar, VECTORSIZE)
        times[Int(STRIAD), j] = @elapsed striad_func(A, B, C, D, VECTORSIZE)
        times[Int(SDAXPY), j] = @elapsed sdaxpy_func(A, B, C, VECTORSIZE)
    end

    for j in 1:NUMBENCH
        for k in 1:NTIMES
            t = times[j, k]
            avgtime[j] = avgtime[j] + t
            mintime[j] = min(mintime[j], t)
            maxtime[j] = max(maxtime[j], t)
        end
        avgtime[j] / NTIMES
    end

    @printf("Function    Rate(MB/s)  Rate(MFlop/s)  Avg time     Min time     Max time\n");
    @printf("%s\n", "-"^80)
    for j in 1:NUMBENCH
        bytes = benchmarks[j].words * sizeof(Float64) * VECTORSIZE
        flops = benchmarks[j].flops * VECTORSIZE

        if flops > 0
            @printf("%-10s%11.2f %11.2f %11.4f  %11.4f  %11.4f\n",
                        benchmarks[j].label,
                        1.0E-06 * bytes/mintime[j],
                        1.0E-06 * flops/mintime[j],
                        avgtime[j],
                        mintime[j],
                        maxtime[j])
        else
            @printf("%-10s%11.2f         -   %11.4f  %11.4f  %11.4f\n",
                        benchmarks[j].label,
                        1.0E-06 * bytes/mintime[j],
                        avgtime[j],
                        mintime[j],
                        maxtime[j])
        end
    end

    @printf("%s\n", "-"^80)
    check(A, B, C, D, VECTORSIZE)

end

bwBench()

max_stars_repo_name: RRZE-HPC/TheBandwidthBenchmarkGame

max_issues_repo_name: RRZE-HPC/TheBandwidthBenchmarkGame

max_forks_repo_name: RRZE-HPC/TheBandwidthBenchmarkGame


text:
#=
*** rewrite for import purposes ***

This is a resnet implementation by Deniz Yuret (https://github.com/denizyuret/Knet.jl/blob/master/examples/resnet/resnet.jl) with minor changes:

Copyright (c) 2015: Deniz Yuret.

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

julia resnet.jl image-file-or-url

This example implements the ResNet-50, ResNet-101 and ResNet-152 models from
'Deep Residual Learning for Image Regocnition', Kaiming He, Xiangyu Zhang,
Shaoqing Ren, Jian Sun, arXiv technical report 1512.03385, 2015.

* Paper url: https://arxiv.org/abs/1512.03385
* Project page: https://github.com/KaimingHe/deep-residual-networks
* MatConvNet weights used here: http://www.vlfeat.org/matconvnet/pretrained
=#

# mode, 0=>train, 1=>test
function resnet50(w,x,ms; mode=1)
    # layer 1
    conv1  = conv4(w[1],x; padding=3, stride=2) .+ w[2]
    bn1    = batchnorm(w[3:4],conv1,ms; mode=mode)
    pool1  = pool(bn1; window=3, stride=2)

    # layer 2,3,4,5
    r2 = reslayerx5(w[5:34], pool1, ms; strides=[1,1,1,1], mode=mode)
    r3 = reslayerx5(w[35:73], r2, ms; mode=mode)
    r4 = reslayerx5(w[74:130], r3, ms; mode=mode) # 5
    r5 = reslayerx5(w[131:160], r4, ms; mode=mode)

    # fully connected layer
    pool5  = pool(r5; stride=1, window=7, mode=2)
    fc1000 = w[161] * mat(pool5) .+ w[162]
end

# mode, 0=>train, 1=>test
function resnet101(w,x,ms; mode=1)
    # layer 1
    conv1 = reslayerx1(w[1:3],x,ms; padding=3, stride=2, mode=mode)
    pool1 = pool(conv1; window=3, stride=2)

    # layer 2,3,4,5
    r2 = reslayerx5(w[4:33], pool1, ms; strides=[1,1,1,1], mode=mode)
    r3 = reslayerx5(w[34:72], r2, ms; mode=mode)
    r4 = reslayerx5(w[73:282], r3, ms; mode=mode)
    r5 = reslayerx5(w[283:312], r4, ms; mode=mode)

    # fully connected layer
    pool5  = pool(r5; stride=1, window=7, mode=2)
    fc1000 = w[313] * mat(pool5) .+ w[314]
end

# mode, 0=>train, 1=>test
function resnet152(w,x,ms; mode=1)
    # layer 1
    conv1 = reslayerx1(w[1:3],x,ms; padding=3, stride=2, mode=mode)
    pool1 = pool(conv1; window=3, stride=2)

    # layer 2,3,4,5
    r2 = reslayerx5(w[4:33], pool1, ms; strides=[1,1,1,1], mode=mode)
    r3 = reslayerx5(w[34:108], r2, ms; mode=mode)
    r4 = reslayerx5(w[109:435], r3, ms; mode=mode)
    r5 = reslayerx5(w[436:465], r4, ms; mode=mode)

    # fully connected layer
    pool5  = pool(r5; stride=1, window=7, mode=2)
    fc1000 = w[466] * mat(pool5) .+ w[467]
end

# Batch Normalization Layer
# works both for convolutional and fully connected layers
# mode, 0=>train, 1=>test
function batchnorm(w, x, ms; mode=1, epsilon=1e-5)
    mu, sigma = nothing, nothing
    if mode == 0
        d = ndims(x) == 4 ? (1,2,4) : (2,)
        s = prod(size(x,d...))
        mu = sum(x,d) / s
        x0 = x .- mu
        x1 = x0 .* x0
        sigma = sqrt(epsilon + (sum(x1, d)) / s)
    elseif mode == 1
        mu = popfirst!(ms)
        sigma = popfirst!(ms)
    end

    # we need value in backpropagation
    push!(ms, value(mu), value(sigma))
    xhat = (x.-mu) ./ sigma
    return w[1] .* xhat .+ w[2]
end

function reslayerx0(w,x,ms; padding=0, stride=1, mode=1)
    b  = conv4(w[1],x; padding=padding, stride=stride)
    bx = batchnorm(w[2:3],b,ms; mode=mode)
end

function reslayerx1(w,x,ms; padding=0, stride=1, mode=1)
    relu.(reslayerx0(w,x,ms; padding=padding, stride=stride, mode=mode))
end

function reslayerx2(w,x,ms; pads=[0,1,0], strides=[1,1,1], mode=1)
    ba = reslayerx1(w[1:3],x,ms; padding=pads[1], stride=strides[1], mode=mode)
    bb = reslayerx1(w[4:6],ba,ms; padding=pads[2], stride=strides[2], mode=mode)
    bc = reslayerx0(w[7:9],bb,ms; padding=pads[3], stride=strides[3], mode=mode)
end

function reslayerx3(w,x,ms; pads=[0,0,1,0], strides=[2,2,1,1], mode=1) # 12
    a = reslayerx0(w[1:3],x,ms; stride=strides[1], padding=pads[1], mode=mode)
    b = reslayerx2(w[4:12],x,ms; strides=strides[2:4], pads=pads[2:4], mode=mode)
    relu.(a .+ b)
end

function reslayerx4(w,x,ms; pads=[0,1,0], strides=[1,1,1], mode=1)
    relu.(x .+ reslayerx2(w,x,ms; pads=pads, strides=strides, mode=mode))
end

function reslayerx5(w,x,ms; strides=[2,2,1,1], mode=1)
    x = reslayerx3(w[1:12],x,ms; strides=strides, mode=mode)
    for k = 13:9:length(w)
        x = reslayerx4(w[k:k+8],x,ms; mode=mode)
    end
    return x
end

function get_params(params, atype)
    len = length(params["value"])
    ws, ms = [], []
    for k = 1:len
        name = params["name"][k]
        value = convert(Array{Float32}, params["value"][k])

        if endswith(name, "moments")
            push!(ms, reshape(value[:,1], (1,1,size(value,1),1)))
            push!(ms, reshape(value[:,2], (1,1,size(value,1),1)))
        elseif startswith(name, "bn")
            push!(ws, reshape(value, (1,1,length(value),1)))
        elseif startswith(name, "fc") && endswith(name, "filter")
            push!(ws, transpose(reshape(value,(size(value,3),size(value,4)))))
        elseif startswith(name, "conv") && endswith(name, "bias")
            push!(ws, reshape(value, (1,1,length(value),1)))
        else
            push!(ws, value)
        end
    end
    map(wi->convert(atype, wi), ws),
    map(mi->convert(atype, mi), ms)
end

# This allows both non-interactive (shell command) and interactive calls like:
# $ julia vgg.jl cat.jpg
# julia> ResNet.main("cat.jpg")
#PROGRAM_FILE=="resnet.jl" && main(ARGS)


text:
# check dimensions of inputs
checkdimensions(sizedest::Tuple{Int, Vararg{Int}}, sizesrc::Tuple{Int, Vararg{Int}}) =
    (sizedest == sizesrc ||
        throw(DimensionMismatch("tried to assign $(sizesrc) sized " *
                                "array to $(sizedest) destination")) )

checkdimensions(dest::AbstractVector, src::AbstractVector) =
    checkdimensions(size(dest), size(src))

checkdimensions(ldest::Int, src::AbstractVector) =
    checkdimensions((ldest, ), size(src))

checkdimensions(kr::AbstractRange, jr::AbstractRange, src::AbstractMatrix) =
    checkdimensions((length(kr), length(jr)), size(src))

# return the bandwidths of A*B
prodbandwidths(A) = bandwidths(A)
prodbandwidths() = (0,0)
prodbandwidths(A...) = broadcast(+, bandwidths.(A)...)

# helper functions in matrix addition routines
function sumbandwidths(A::AbstractMatrix, B::AbstractMatrix)
    max(bandwidth(A, 1), bandwidth(B, 1)), max(bandwidth(A, 2), bandwidth(B, 2))
end


text:
using CSV
using DataFrames
using Dates
using UnicodePlots

function load_data(filename)
    df = DataFrame(CSV.File(filename))
    df.dayname = Dates.dayname.(df.Day)
    df.dayofweek = Dates.dayofweek.(df.Day)
    return df
end

function plot(df)
    xticks = collect(df.Day[1]:Dates.Day(7):df.Day[lastindex(df.Day)])
    # plt = lineplot(df[:, 1], df[:, 2], width=80, title="Data is Beautiful. Hangovers Are Not", ylabel="Google Search Volume (Indexed, 100=Maximum)", color=:green)
    plt = lineplot(df[:, 1], df[:, 2], width=50, title="Data is Beautiful. Hangovers Are Not", ylabel="", color=:green)
    return plt
end

function plot(df, plt)
    xticks = collect(df.Day[1]:Dates.Day(7):df.Day[lastindex(df.Day)])
    if @isdefined plt
        plt = lineplot!(plt, df[:, 1], df[:, 2], color=:green)
    else
        plt = lineplot(df[:, 1], df[:, 2], width=80, title="Data is Beautiful. Hangovers Are Not", ylabel="Google Search Volume (Indexed, 100=Maximum)")
    end
    # lineplot!(plt, xticks=xticks)
    return plt
end

function weekends(df)
    weekends_df = filter(row -> row[:dayofweek] >= 6, df)
    return weekends_df
end

function vspan_xs(df)
    new_df = DataFrame()
    for row in Tables.rowtable(df)
        push!(new_df, row)
        push!(new_df, row)
    end
    return new_df
end

function vspan_ys(total)
    v = []
    for i in 1:total
        push!(v, 0)
        push!(v, 100)
    end
    return v
end

function plot_vspan(plt, df)
    weekend_df = weekends(df)
    xs = vspan_xs(weekend_df)
    ys = vspan_ys(nrow(weekend_df))
    for i in 1:2:nrow(xs)-1
        lineplot!(plt, xs[i:i+1, 1], ys[i:i+1], color=:red)
        # println(xs[i:i+1, 1], ys[i:i+1])
    end

    annotate!(plt, :r, 1, "🍺 Sat & Sun", color=:red)
    annotate!(plt, :r, 3, "Search Volume", color=:green)

    return plt
end

max_stars_repo_name: katnoria/data-books

max_issues_repo_name: katnoria/data-books

max_forks_repo_name: katnoria/data-books