text:
"""
    agents_plots_complete(property_plot::Array{Tuple}, model::AbstractModel)

Plots the agents_data_complete() results in your browser.

# Parameters

* property_plot: An array of tuples. The first element of the tuple is the agent property you would like to plot (as a symbol). The second element of the tuple is the type of plot you want on that property (as a symbol). Example: [(:wealth, :hist)]. Available plot types: :hist.

You can add more plot types by adding more if statements in the function.
"""
function agents_plots_complete(property_plot::Array, model::AbstractModel)
  properties = [i[1] for i in property_plot]
  data = agents_data_complete(properties, model)
  for (property, pt) in property_plot
    if pt == :hist
      # data |> @vlplot(:bar, x={property, bin=true, title=property}, y={"count()", title="Frequency"}) # This is a bug and doesn't work
      data |> @vlplot(:bar, x=property, y="count()") |> save("histogram.pdf") # this is temporary until the above bug is fixed
    end
  end
end

"""
    visualize_data(data::DataFrame)

sends data to Data Voyager.
"""
function visualize_data(data::DataFrame)
  v = Voyager(data)
end

"""
    plot_locs(g, dims::Tuple{Integer,Integer,Integer})

Return three arrays for x, y, z coordinates of each node
"""
function node_locs(g, dims::Tuple{Integer,Integer,Integer})
  coords = []
  for nn in 1:nv(g)
    push!(coords, vertex_to_coord(nn, dims))
  end
  locs_x = [Float64(i[1]) for i in coords]
  locs_y = [Float64(i[2]) for i in coords]
  locs_z = [Float64(i[3]) for i in coords]
  return locs_x, locs_y, locs_z
end

"""
    plot_locs(g, dims::Tuple{Integer,Integer})

Return arrays for x, y coordinates of each node
"""
function node_locs(g, dims::Tuple{Integer,Integer})
  coords = []
  for nn in 1:nv(g)
    push!(coords, vertex_to_coord(nn, dims))
  end
  locs_x = [Float64(i[1]) for i in coords]
  locs_y = [Float64(i[2]) for i in coords]
  return locs_x, locs_y
end

"""
    visualize_2D_agent_distribution(data::DataFrame, model::AbstractModel, position_column::Symbol; types::Symbol=:id, savename::AbstractString="2D_agent_distribution", saveloc::AbstractString="./", cc::Dict=Dict())

Plots the distribution of agents on a 2D grid.

Plots are saved as PDF files under the name given by the `savename` argument. Optionally, choose a path to save the plots using the `saveloc` argument. The default behavior is to save in the current directory, where the code is run.

* You should provide `position_colomn` which is the name of the column that holds agent positions.
* If agents have different types and you want each type to be a different color, provide types=<column name>. Use a dictionary (the `cc` argument) to pass colors for each type. You may choose any color name from the [list of colors on Wikipedia](https://en.wikipedia.org/wiki/Lists_of_colors).
"""
function visualize_2D_agent_distribution(data::DataFrame, model::AbstractModel, position_column::Symbol; types::Symbol=:id, savename::AbstractString="2D_agent_distribution", saveloc::AbstractString="./", cc::Dict=Dict())
  g = model.space.space
  locs_x, locs_y, = node_locs(g, model.space.dimensions)

  # base node color is light grey
  nodefillc = [RGBA(0.1,0.1,0.1,.1) for i in 1:gridsize(model.space.dimensions)]

  # change node color given the position of the agents. Automatically uses any columns with names: pos, or pos_{some number}
  # TODO a new plot where the alpha value of a node corresponds to the value of an individual on a node
  if types == :id  # there is only one type
    pos = position_column
    d = by(data, pos, N = pos => length)
    maxval = maximum(d[!, :N])
    nodefillc[d[pos]] .= [RGBA(0.1, 0.1, 0.1, i) for i in  (d[!, :N] ./ maxval) .- 0.001]
  else  # there are different types of agents based on the values of the "types" column
    dd = dropmissing(data[:, [position_column, types]])
    unique_types = sort(unique(dd[!, types]))
    pos = position_column
    if length(cc) == 0
      colors = colorrgb(length(unique_types))
      colordict = Dict{Any, Tuple}()
      colorvalues = collect(values(colors))
      for ut in 1:length(unique_types)
        colordict[unique_types[ut]] = colorvalues[ut]
      end
    else
      colors = colorrgb(collect(values(cc)))
      colordict = Dict{Any, Tuple}()
      for key in keys(cc)
        colordict[key] = colors[cc[key]]
      end
    end
    colorrev = Dict(v=>k for (k,v) in colors)
    for index in 1:length(unique_types)
      tt = unique_types[index]
      d = by(dd[dd[!, types] .== tt, :], pos, N = pos => length)
      maxval = maximum(d[!, :N])
      # colormapname = "L$(index+1)"  # a linear colormap
      # (cmapc, name, desc) = cmap(colormapname, returnname=true)
      # nodefillc[d[pos]] .= [cmapc[round(Int64, i*256)] for i in  (d[:N] ./ maxval) .- 0.001]
      # println("$tt: $name")
      nodefillc[d[!, pos]] .= [RGBA(colordict[tt][1], colordict[tt][2], colordict[tt][3], i) for i in  (d[!, :N] ./ maxval) .- 0.001]
      println("$tt: $(colorrev[colordict[tt]])")
    end
  end

  NODESIZE = 0.8/sqrt(gridsize(model))
  draw(PDF(joinpath(saveloc, "$savename.pdf")), gplot(g, locs_x, locs_y, nodefillc=nodefillc, edgestrokec=RGBA(0.1,0.1,0.1,.1), NODESIZE=NODESIZE))
end

"""
    colorrgb(color_names::Array)
Returns a dictionary of each colorname and its RGB values. See colors and names on [list of colors on Wikipedia](https://en.wikipedia.org/wiki/Lists_of_colors)
"""
function colorrgb(color_names::Array)
  script_path = splitdir(realpath(@__FILE__))[1]
  f = joinpath(script_path, "color_names.csv")
  ff = CSV.File(f)
  rgb_dict = Dict{AbstractString, Tuple}()
  for row in ff
    if row.cname in color_names || row.cname2 in color_names
      rgb_dict[row.cname] = (row.R/256, row.G/256, row.B/256)
    end
  end
  if length(rgb_dict) < length(color_names)
    simm = intersect(keys(rgb_dict), color_names)
    for ss in keys(rgb_dict)
      if !in(ss, color_names)
        println("$ss is not a valid color name!")
      end
    end
    throw("Provide valid colornames.")
  end
  return rgb_dict
end

"""
    colornames(n::Integer)

Returns n random colors as a dictionary (Dict{colorname=>rgb})
"""
function colorrgb(n::Integer)
  script_path = splitdir(realpath(@__FILE__))[1]
  f = joinpath(script_path, "color_names.csv")
  ff = CSV.File(f)
  randcolors = rand(1:length(ff), n)
  rgb_dict = Dict{AbstractString, Tuple}()
  for (index,row) in enumerate(ff)
    if index in randcolors
      rgb_dict[row.cname] = (row.R/256, row.G/256, row.B/256)
    end
  end
  return rgb_dict
end

"""
    visualize_1DCA(data::DataFrame, model::AbstractModel, position_column::Symbol, status_column::Symbol, nrows::Integer; savename::AbstractString="CA_1D", saveloc::AbstractString="./")

Visualizes data of a 1D cellular automaton and saves it in a PDF file under the name given by the `savename` argument. Optionally provide a location for the plots to be saved using the `saveloc` argument. The default behavior is to save them in the currently active directory.

* `data` are the result of multiple runs of the simulation.
* `position_column` is the field of the agent that holds their position.
* `status_column` is the field of the agents that holds their status.
* `nrows` is the number of times the model was run.
"""
function visualize_1DCA(data::DataFrame, model::AbstractModel, position_column::Symbol, status_column::Symbol, nrows::Integer; savename::AbstractString="CA_1D", saveloc::AbstractString="./")
  dims = (model.space.dimensions[1], nrows)
  g = Agents.grid2D(dims[1], dims[2])
  locs_x, locs_y = node_locs(g, dims)

  # base node color is light grey
  nodefillc = [RGBA(0.1,0.1,0.1,.1) for i in 1:Agents.gridsize(dims)]

  for row in 1:nrows
    pos = Symbol(string(position_column)*"_$row")
    status = Symbol(string(status_column)*"_$row")
    newcolors = [RGBA(0.1, 0.1, 0.1, 0.01) for i in 1:dims[1]]
    for ll in 1:dims[1]
      if data[!, status][ll] == "1"
        newcolors[ll] = RGBA(0.1, 0.1, 0.1, 1.0)
      end
    end
    nodefillc[(dims[1]*row)-(dims[1]-1):dims[1]*row] .= newcolors
  end

  NODESIZE = 1/sqrt(gridsize(dims))
  draw(PDF(joinpath(saveloc, "$savename.pdf")), gplot(g, locs_x, locs_y, nodefillc=nodefillc, edgestrokec=RGBA(0.1,0.1,0.1,0.01), NODESIZE=NODESIZE))

end

"""
    visualize_2DCA(data::DataFrame, model::AbstractModel, position_column::Symbol, status_column::Symbol, runs::Integer; savename::AbstractString="CA_2D")

Visualizes data of a 2D cellular automaton and saves it in PNG format under the name given by the `savename` argument. Optionally provide a location for the plots to be saved using the `saveloc` argument. The default behavior is to save them in the currently active directory.

`data` are the result of multiple runs of the simulation. `position_column` is the field of the agent that holds their position. `status_column` is the field of the agents that holds their status. `runs` is the number of times the simulation was run.
"""
function visualize_2DCA(data::DataFrame, model::AbstractModel, position_column::Symbol, status_column::Symbol, runs::Integer; savename::AbstractString="CA_2D", saveloc::AbstractString="./")
  dims = model.space.dimensions
  g = model.space.space
  locs_x, locs_y = node_locs(g, dims)
  NODESIZE = 0.8/sqrt(gridsize(dims))

  for r in 1:runs
    # base node color is light grey
    nodefillc = [RGBA(0.1,0.1,0.1,.1) for i in 1:Agents.gridsize(dims)]
    stat = Symbol(string(status_column)*"_$r")
    nonzeros = findall(a-> a =="1", data[!, stat])

    nodefillc[nonzeros] .= RGBA(0.1, 0.1, 0.1, 1)

    draw(PNG(joinpath(saveloc, "$(savename)_$r.png")), gplot(g, locs_x, locs_y, nodefillc=nodefillc, edgestrokec=RGBA(0.1,0.1,0.1,0.01), NODESIZE=NODESIZE))
  end
end


text:
@safetestset "isbic" begin
    using Validator

    @test Validator.isbic("SBICKEN1345") === true
    @test Validator.isbic("SBICKEN1") === true
    @test Validator.isbic("SBICKENY") === true
    @test Validator.isbic("SBICKEN1YYP") === true

    @test Validator.isbic("SBIC23NXXX") === false
    @test Validator.isbic("S23CKENXXXX") === false
    @test Validator.isbic("SBICKENXX") === false
    @test Validator.isbic("SBICKENXX9") === false
    @test Validator.isbic("SBICKEN13458") === false
    @test Validator.isbic("SBICKEN") === false
    @test Validator.isbic("SBI") === false
end

ext: jl


text:
module mouse

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
module GroundDataL
"""
Module for the download and processing of atmospheric data gathered by measuring stations located in Lombardia, Italy.
"""
#=
Si possono ottenere dati dal link (ALL CSV):
    Stazioni:
        https://www.dati.lombardia.it/resource/wkyn-7szs.csv
    Stazioni meteo provincia di Milano:
        https://www.dati.lombardia.it/resource/my9y-8ykb.csv

    Stazioni meteo:
        https://www.dati.lombardia.it/resource/nf78-nj6b.csv
    Dati sensori meteo:
        https://www.dati.lombardia.it/resource/647i-nhxk.csv
    Stazioni aria:
        https://www.dati.lombardia.it/resource/ib47-atvt.csv
    Dati sensori aria:
        https://www.dati.lombardia.it/resource/nicp-bhqi.csv
Si possono ottenere direttamente CVS dei dati.
Dati sensori meteo sembra fare al caso nostro ma ha un campo che indica cosa rappresenta il dato che è encoded
=#

using CSV
using DataFrames
using HTTP

export getData,
       getRegionAttributes, getRegionIds, getRegionStationInfo

"""
    getRegionAttributes( [ type::Symbol=:METEO ] )

Obtain the names of the columns of the region's dataframe required by `GroundData.createMap`'s `attributes` parameter to create `GroundData.standardize`'s `map` parameter.
"""
function getRegionAttributes( type::Symbol=:METEO )
    return type == :METEO ? [ :tipologia, :unit_dimisura, :valore, nothing, :data, :lng, :lat, :quota, :stato, nothing, :rmh ] :
        type == :AIRQUALITY ? [ :nometiposensore, :unitamisura, :valore, nothing, :data, :lng, :lat, :quota, :stato, nothing ] :
            throw(DomainError(type, "`type` must be either `:METEO` OR `:AIRQUALITY`."))
end

"""
    getRegionAttributes( [ type::Symbol=:METEO ] )

Obtain the names of the columns of the dataframe required for `GroundData.standardize`'s `bridge` parameter.
"""
function getRegionIds( type::Symbol=:METEO )
    if type != :METEO && type != :AIRQUALITY
        throw(DomainError(type, "`type` must be either `:METEO` OR `:AIRQUALITY`."))
    end
    return :idsensore
end

"""
    getRegionStationInfo( [ type::Symbol=:METEO  ] )

Obtain the names of the columns of the region's stations dataframe required by `GroundData.createMap`'s `attributes` parameter to be used in `GroundData.generateUuidsTable`.
"""
function getRegionStationsInfo( type::Symbol=:METEO )
    if type != :METEO && type != :AIRQUALITY
        throw(DomainError(type, "`type` must be either `:METEO` OR `:AIRQUALITY`."))
    end
    return [ :idstazione, :nomestazione, :lng, :lat ]
end

"""
    getData(; <keyword arguments> )

Obtain data of category `type` and source of category `kind`.

# Arguments
- `type::Symbol=:METEO`: defines the type of data to be downloaded, it must either be `:METEO` or `:AIRQUALITY`.
- `kind::Symbol=:STATIONS`: defines if the data to be downloaded has to regard information on the stations or their actual measurements, it must either be `:STATIONS` or `:SENSORS`.
"""
function getData(; type::Symbol=:METEO, kind::Symbol=:STATIONS )
    str = type == :METEO ? (
            kind == :STATIONS ? "nf78-nj6b" :
            kind == :SENSORS ? "647i-nhxk" :
            throw(DomainError(kind, "`kind` must either be `:STATIONS` or `:SENSORS`."))
        ) :
        type == :AIRQUALITY ? (
            kind == :STATIONS ? "ib47-atvt" :
            kind == :SENSORS ? "nicp-bhqi" :
            throw(DomainError(kind, "`kind` must either be `:STATIONS` or `:SENSORS`."))
        ) :
        throw(DomainError(type, "`type` must either be `:METEO` or `:AIRQUALITY`."))
    data = HTTP.get( "https://www.dati.lombardia.it/resource/$str.csv" )
    df = CSV.read( data.body, DataFrame )
    if type ==:METEO && kind == :SENSORS
        insertcols!( df, :rmh => "0m" )
    end
    return df
end

#   ressta = getData( type=:METEO, source=:STATIONS )
#   ressen = getData( type=:METEO, source=:SENSORS )
#   ressta = getData( type=:AIRQUALITY, source=:STATIONS )
#   ressen = getData( type=:AIRQUALITY, source=:SENSORS )

end # module

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
#  Defines the types "Population" and "PopList".
#  Defines the composite types "cont_var_result_type" and "variant_type"
export variant_type, fitness_location_type
const Population = Array{Int64,1}
const PopList = Array{Population,1}
try
  using Distributed
  using Random
  using Dates
catch
end

mutable struct variant_type
  fitness::Float64    # The fitness of this variant
  attributes::Vector{Float64}   # attributes of the variant
end

mutable struct cont_var_result_type
  N_list::Vector{Int64}
  num_attributes_list::Vector{Int64}
  mutation_stddev_list::Vector{Float64}  # list of mutation_stddev values
  N_mut_list::Vector{Float64}  # list of mutation_stddev values
  num_trials::Int64  # Number of times to repeat simulation for each setting of the parameters
  N::Int64   # meta-population size
  num_subpops::Int64   # number of subpopulations
  num_attributes::Int64  # number of attributes of a variant
  ngens::Int64  # number of generations after burn-in
  int_burn_in::Int64
  mutation_stddev::Float64  # standard deviation of mutation distribution of mutation perturbations
  ideal::Float64             #ideal value
  fit_slope::Float64         # fitness = 1.0/(fit_slope*distance(attributes,ideal)+1.0)
  neutral::Bool              # If true, fitness = 1
  fitness_mean::Float64      # average of fitnesses over subpopulations and generations
  fitness_median::Float64      # average of fitnesses over subpopulations and generations
  fitness_coef_var::Float64  # average coefficent of variation of fitnesss over subpopulations and generations
  attribute_mean::Float64    # average of attributes over attributes, subpopulations, and generations
  attribute_median::Float64    # average of attributes over attributes, subpopulations, and generations
  attribute_coef_var::Float64 # average coefficent of variation of attributes over attributes, subpopulations and generations
  neg_count::Int64        # Number of fitness differences < -1.0/N
  neg_neutral::Int64      # Number of fitness differences >= -1.0/N  and < 0.0
  pos_neutral::Int64     # Number of fitness differences >= 0.0 and < 1.0/N
  pos_count::Int64        # Number of fitness differences >= 1.0/N
end


text:
using LazyArrays, Test
import LazyArrays: Add, AddArray, MulAdd, materialize!, MemoryLayout, ApplyLayout

@testset "Add/Subtract" begin
    @testset "Add" begin
        @testset "Mul-Add" begin
            A = Add(randn(5,5), randn(5,5))
            b = randn(5)
            c = similar(b)
            fill!(c,NaN)
            @test (c .= @~ A*b) ≈ A.args[1]*b + A.args[2]*b
        end

        @testset "gemv Float64" begin
            for A in (AddArray(randn(5,5), randn(5,5)),
                    AddArray(randn(5,5), view(randn(9, 5), 1:2:9, :))),
                b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9))

                Ã = Array(A)
                c = similar(b)
                fill!(c,NaN)

                materialize!(MulAdd(1.0, A, b, 0.0, c))
                @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c))

                c = similar(b)
                fill!(c,NaN)
                c .= @~ A*b
                @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c))

                copyto!(c, @~ A*b)
                @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c))

                b̃ = Array(b)
                copyto!(b̃, Mul(A,b̃))
                @test c ≈ b̃

                c .= @~ 2.0 * A * b
                @test c ≈ BLAS.gemv!('N', 2.0, Ã, b, 0.0, similar(c))

                c = Array(b)
                c .= @~ A*b + c
                @test c ≈ BLAS.gemv!('N', 1.0, Ã, b, 1.0, Array(b))

                c = Array(b)
                c .= @~ A*b + 2.0 * c
                @test c ≈ BLAS.gemv!('N', 1.0, Ã, b, 2.0, Array(b))

                c = Array(b)
                c .= @~ 2.0 * A*b + c
                @test c ≈ BLAS.gemv!('N', 2.0, Ã, b, 1.0, Array(b))

                c = Array(b)
                c .= @~ 3.0 * A*b + 2.0 * c
                @test c ≈ BLAS.gemv!('N', 3.0, Ã, b, 2.0, Array(b))

                d = similar(c)
                c = Array(b)
                d .= @~ 3.0 * A*b + 2.0 * c
                @test d ≈ BLAS.gemv!('N', 3.0, Ã, b, 2.0, Array(b))
            end
        end

        @testset "gemm" begin
            for A in (AddArray(randn(5,5), randn(5,5)),
                    AddArray(randn(5,5), view(randn(9, 5), 1:2:9, :))),
                B in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                    view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5))

                Ã = Array(A)
                C = similar(B)

                C .= @~ A*B
                @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 0.0, similar(C))

                B .= @~ A*B
                @test C ≈ B

                C .= @~ 2.0 * A*B
                @test C ≈ BLAS.gemm!('N', 'N', 2.0, Ã, B, 0.0, similar(C))

                C = Array(B)
                C .= @~ A*B + C
                @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 1.0, Array(B))

                C = Array(B)
                C .= @~ A*B + 2.0 * C
                @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 2.0, Array(B))

                C = Array(B)
                C .= @~ 2.0 * A*B + C
                @test C ≈ BLAS.gemm!('N', 'N', 2.0, Ã, B, 1.0, Array(B))

                C = Array(B)
                C .= @~ 3.0 * A*B + 2.0 * C
                @test C ≈ BLAS.gemm!('N', 'N', 3.0, Ã, B, 2.0, Array(B))

                d = similar(C)
                C = Array(B)
                d .= @~ 3.0 * A*B + 2.0 * C
                @test d ≈ BLAS.gemm!('N', 'N', 3.0, Ã, B, 2.0, Array(B))
            end
        end
    end

    @testset "Subtract" begin
        @testset "Mul-Subtract" begin
            A = applied(-, randn(5,5), randn(5,5))
            b = randn(5)
            c = similar(b)
            fill!(c,NaN)
            @test (c .= @~ A*b) ≈ A.args[1]*b - A.args[2]*b
        end

        @testset "gemv Float64" begin
            for A in (ApplyArray(-, randn(5,5), randn(5,5)),
                    ApplyArray(-, randn(5,5), view(randn(9, 5), 1:2:9, :))),
                b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9))

                Ã = Array(A)
                c = similar(b)
                fill!(c,NaN)

                materialize!(MulAdd(1.0, A, b, 0.0, c))
                @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c))

                c = similar(b)
                fill!(c,NaN)
                c .= @~ A*b
                @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c))

                copyto!(c, @~ A*b)
                @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c))

                b̃ = Array(b)
                copyto!(b̃, Mul(A,b̃))
                @test c ≈ b̃

                c .= @~ 2.0 * A * b
                @test c ≈ BLAS.gemv!('N', 2.0, Ã, b, 0.0, similar(c))

                c = Array(b)
                c .= @~ A*b + c
                @test c ≈ BLAS.gemv!('N', 1.0, Ã, b, 1.0, Array(b))

                c = Array(b)
                c .= @~ A*b + 2.0 * c
                @test c ≈ BLAS.gemv!('N', 1.0, Ã, b, 2.0, Array(b))

                c = Array(b)
                c .= @~ 2.0 * A*b + c
                @test c ≈ BLAS.gemv!('N', 2.0, Ã, b, 1.0, Array(b))

                c = Array(b)
                c .= @~ 3.0 * A*b + 2.0 * c
                @test c ≈ BLAS.gemv!('N', 3.0, Ã, b, 2.0, Array(b))

                d = similar(c)
                c = Array(b)
                d .= @~ 3.0 * A*b + 2.0 * c
                @test d ≈ BLAS.gemv!('N', 3.0, Ã, b, 2.0, Array(b))
            end
        end

        @testset "gemm" begin
            for A in (ApplyArray(-, randn(5,5), randn(5,5)),
                    ApplyArray(-, randn(5,5), view(randn(9, 5), 1:2:9, :))),
                B in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                    view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5))

                Ã = Array(A)
                C = similar(B)

                C .= @~ A*B
                @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 0.0, similar(C))

                B .= @~ A*B
                @test C ≈ B

                C .= @~ 2.0 * A*B
                @test C ≈ BLAS.gemm!('N', 'N', 2.0, Ã, B, 0.0, similar(C))

                C = Array(B)
                C .= @~ A*B + C
                @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 1.0, Array(B))

                C = Array(B)
                C .= @~ A*B + 2.0 * C
                @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 2.0, Array(B))

                C = Array(B)
                C .= @~ 2.0 * A*B + C
                @test C ≈ BLAS.gemm!('N', 'N', 2.0, Ã, B, 1.0, Array(B))

                C = Array(B)
                C .= @~ 3.0 * A*B + 2.0 * C
                @test C ≈ BLAS.gemm!('N', 'N', 3.0, Ã, B, 2.0, Array(B))

                d = similar(C)
                C = Array(B)
                d .= @~ 3.0 * A*B + 2.0 * C
                @test d ≈ BLAS.gemm!('N', 'N', 3.0, Ã, B, 2.0, Array(B))
            end
        end
    end

    @testset "view" begin
        A = ApplyArray(+, randn(5,5), randn(5,5))
        Ã = Array(A)
        V = view(A,2:3,2:4)
        @test MemoryLayout(typeof(V)) == ApplyLayout{typeof(+)}()
        b = randn(3)
        c = randn(2)
        materialize!(MulAdd(1.0, V, b, 0.0, c))
        @test c ≈ view(Ã,2:3,2:4)*b

        A = ApplyArray(-, randn(5,5), randn(5,5))
        Ã = Array(A)
        V = view(A,2:3,2:4)
        @test MemoryLayout(typeof(V)) == ApplyLayout{typeof(-)}()
        b = randn(3)
        c = randn(2)
        materialize!(MulAdd(1.0, V, b, 0.0, c))
        @test c ≈ view(Ã,2:3,2:4)*b
    end
end


text:
#########################################################################
#  Physical and mathematical constants
module constants
include("user_settings.jl")

g = 9.8;   # gravitational acceleration: m/s^2
cp = 1004.0;  # specific heat at constant p: J/(kg K)
rd = 287.0;  # dry gas constant:  J/(kg K)
c_v = cp-rd;  # specific heat at constant V: J/(kg K)
xk = rd/cp;  # kappa, the ratio of rd/cp
p0 = 100000.0;  # reference pressure:  Pa
lv = 2500000.0;  # latent heat of vaporization:  J/kg
trigpi = pi; # trigonmetric pi
rhol = 1000.0; # Density of Rain drop (1000 kg/m^3)

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

## Initialize Variable Arrays
nz = user_settings.nz;
nx = user_settings.nx;

zw = zeros(1,nz);
zu = zeros(1,nz);
xu = zeros(1,nx);
xs = zeros(1,nx);
rhow=zeros(1,nz);
rhou=zeros(1,nz);
tbv=zeros(1,nz);
ub=zeros(1,nz);
tb=zeros(1,nz);
qb=zeros(1,nz);
pib=zeros(1,nz);
qbs=zeros(1,nz);
rhb=zeros(1,nz);
qvm=zeros(nx,nz);
qv=zeros(nx,nz);
qvp=zeros(nx,nz);
qvtot = zeros(nx,nz);
qcm=zeros(nx,nz);
qc=zeros(nx,nz);
qcp=zeros(nx,nz);
qrainm=zeros(nx,nz);
qrain=zeros(nx,nz);
qrainp=zeros(nx,nz);
thm=zeros(nx,nz);
th=zeros(nx,nz);
thp=zeros(nx,nz);
thvm=zeros(nx,nz);
thv=zeros(nx,nz);
thvp=zeros(nx,nz);
um=zeros(nx,nz);
u=zeros(nx,nz);
up=zeros(nx,nz);
wm=zeros(nx,nz);
w=zeros(nx,nz);
wp=zeros(nx,nz);
pim=zeros(nx,nz);
pic=zeros(nx,nz);
pip=zeros(nx,nz);
pprt=zeros(nx,nz);

end

max_stars_repo_name: jeffreycunn/toy-atmospheric-model

max_issues_repo_name: jeffreycunn/toy-atmospheric-model

max_forks_repo_name: jeffreycunn/toy-atmospheric-model


text:
"""
Struct for DomainwallFermion
"""
struct DomainwallFermion_5D_wing{NC,WilsonFermion} <: Abstract_DomainwallFermion_5D{NC,WilsonFermion} #<: AbstractFermionfields_5D{NC}
    w::Array{WilsonFermion,1}
    NC::Int64
    NX::Int64
    NY::Int64
    NZ::Int64
    NT::Int64
    L5::Int64
    Dirac_operator::String
    NWilson::Int64

    function DomainwallFermion_5D_wing(L5,NC::T,NX::T,NY::T,NZ::T,NT::T) where T<: Integer
        x = WilsonFermion_4D_wing(NC,NX,NY,NZ,NT)
        xtype = typeof(x)
        w = Array{xtype,1}(undef,L5)
        w[1] = x
        for i=2:L5
            w[i] = similar(x)
        end
        #println(w[2][1,1,1,1,1,1])
        NWilson = length(x)
        Dirac_operator = "Domainwall"
        return new{NC,xtype}(w,NC,NX,NY,NZ,NT,L5,Dirac_operator,NWilson)
    end

end

function Base.similar(x::DomainwallFermion_5D_wing{NC,WilsonFermion} ) where {NC,WilsonFermion}
    return DomainwallFermion_5D_wing(x.L5,NC,x.NX,x.NY,x.NZ,x.NT)
end

#=

function apply_J!(xout::DomainwallFermion_5D_wing{NC,WilsonFermion} ,
    x::DomainwallFermion_5D_wing{NC,WilsonFermion}) where  {NC,WilsonFermion}
    clear_fermion!(xout)

    L5 = xout.L5
    for i5=1:L5
        j5 = L5-i5+1
        substitute_fermion!(xout.w[i5],x.w[j5])
    end
end

function apply_P!(xout::DomainwallFermion_5D_wing{NC,WilsonFermion} ,
    x::DomainwallFermion_5D_wing{NC,WilsonFermion}) where  {NC,WilsonFermion}
    L5 = xout.L5
    clear_fermion!(xout)

    for i5=1:L5
        j5 = i5
        #P_- -> P_+ in this definition
        mul_1plusγ5x_add!(xout.w[i5],x.w[j5],1)
        set_wing_fermion!(xout.w[i5])

        #P_+ -> P_- in this definition
        j5 = i5+1
        j5 += ifelse(j5 > L5,-L5,0)
        mul_1minusγ5x_add!(xout.w[i5],x.w[j5],1)
        set_wing_fermion!(xout.w[i5])
    end
end

function apply_Pdag!(xout::DomainwallFermion_5D_wing{NC,WilsonFermion} ,
    x::DomainwallFermion_5D_wing{NC,WilsonFermion}) where  {NC,WilsonFermion}
    L5 = xout.L5
    clear_fermion!(xout)

    for i5=1:L5
        j5 = i5
        #P_- -> P_+ in this definition
        mul_1plusγ5x_add!(xout.w[i5],x.w[j5],1)
        set_wing_fermion!(xout.w[i5])

        #P_+ -> P_- in this definition
        j5 = i5-1
        j5 += ifelse(j5 < 1,L5,0)
        mul_1minusγ5x_add!(xout.w[i5],x.w[j5],1)
        set_wing_fermion!(xout.w[i5])
    end
end

function D5DWx!(xout::DomainwallFermion_5D_wing{NC,WilsonFermion} ,U::Array{G,1},
    x::DomainwallFermion_5D_wing{NC,WilsonFermion} ,m,A,L5) where  {NC,WilsonFermion,G <: AbstractGaugefields}

    #temp = temps[4]
    #temp1 = temps[1]
    #temp2 = temps[2]
    clear_fermion!(xout)
    ratio = 1
    #ratio = xout.L5/L5
    if L5 != xout.L5
        @assert L5 % 2 == 0
        irange = Int64[]
        irange_out = Int64[]
        #irange = 1:L5
        #irange_out = (L5+1):xout.L5

        for i5=1:xout.L5
            if i5 <= div(L5,2) || i5 >= xout.L5-div(L5,2)+1
                push!(irange,i5)
            else
                push!(irange_out,i5)
            end

        end

        #for i5 in irange_out
        #    axpy!(1,x.w[i5],xout.w[i5])
        #end
    else
        irange = 1:L5
    end

    for i5 in irange
        j5=i5
        D4x!(xout.w[i5],U,x.w[j5],A,4) #Dw*x
        #Dx!(xout.w[i5],U,x.w[j5],A) #Dw*x
        #Wx!(xout.w[i5],U,x.w[j5],temps) #Dw*x
        #1/(2*A.κ)
        massfactor = -(1/(2*A.κ) + 1)
        set_wing_fermion!(xout.w[i5])
        #add!(ratio,xout.w[i5],ratio,x.w[j5]) #D = x + Ddagw*x
        add!(ratio,xout.w[i5],ratio*massfactor,x.w[j5]) #D = x + Dw*x
        set_wing_fermion!(xout.w[i5])

        #println("xout ",xout.w[i5][1,1,1,1,1,1])

        j5=i5+1
        if 1 <= j5 <= xout.L5
            #-P_- -> - P_+ :gamma_5 of LTK definition
            if xout.L5 != 2
                #mul_1minusγ5x_add!(xout.w[i5],x.w[j5],-1*ratio)
                mul_1plusγ5x_add!(xout.w[i5],x.w[j5],ratio)
                set_wing_fermion!(xout.w[i5])

            end
        end

        j5=i5-1
        if 1 <= j5 <= xout.L5
            #-P_+ -> - P_- :gamma_5 of LTK definition
            if xout.L5 != 2
                #mul_1plusγ5x_add!(xout.w[i5],x.w[j5],-1*ratio)
                mul_1minusγ5x_add!(xout.w[i5],x.w[j5],ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

        if xout.L5 != 1
            if i5==1
                j5 = xout.L5
                #mul_1plusγ5x_add!(xout.w[i5],x.w[j5],m*ratio)
                mul_1minusγ5x_add!(xout.w[i5],x.w[j5],-m*ratio)
                set_wing_fermion!(xout.w[i5])
            end

            if i5== xout.L5
                j5 = 1
                #mul_1minusγ5x_add!(xout.w[i5],x.w[j5],m*ratio)
                mul_1plusγ5x_add!(xout.w[i5],x.w[j5],-m*ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

        #println("xout ",xout.w[i5][1,1,1,1,1,1])

    end

    if L5 != xout.L5
        for i5 in irange_out
            axpy!(1,x.w[i5],xout.w[i5])
        end
    end

    set_wing_fermion!(xout)

    return
end

function D5DWdagx!(xout::DomainwallFermion_5D_wing{NC,WilsonFermion} ,U::Array{G,1},
    x::DomainwallFermion_5D_wing{NC,WilsonFermion} ,m,A,L5) where  {NC,WilsonFermion,G <: AbstractGaugefields}

    #temp = temps[4]
    #temp1 = temps[1]
    #temp2 = temps[2]
    clear_fermion!(xout)
    ratio = 1
    #ratio = xout.L5/L5

    if L5 != xout.L5
        @assert L5 % 2 == 0
        irange = Int64[]
        irange_out = Int64[]
        #irange = 1:L5
        #irange_out = (L5+1):xout.L5

        for i5=1:xout.L5
            if i5 <= div(L5,2) || i5 >= xout.L5-div(L5,2)+1
                push!(irange,i5)
            else
                push!(irange_out,i5)
            end

        end

        #for i5 in irange_out
        #    axpy!(1,x.w[i5],xout.w[i5])
        #end
    else
        irange = 1:L5
    end

    for i5 in irange
        j5=i5
        #Ddagx!(xout.w[i5],U,x.w[j5],A) #Ddagw*x
        D4dagx!(xout.w[i5],U,x.w[j5],A,4) #Dw*x
        #Wx!(xout.w[i5],U,x.w[j5],temps) #Dw*x
        #1/(2*A.κ)
        massfactor = -(1/(2*A.κ) + 1)
        #println(massfactor)

        #Wdagx!(xout.w[i5],U,x.w[j5],temps) #Ddagw*x
        set_wing_fermion!(xout.w[i5])
        add!(ratio,xout.w[i5],ratio*massfactor,x.w[j5]) #D = x + Dw*x
        #add!(ratio,xout.w[i5],ratio,x.w[j5]) #D = x + Ddagw*x
        set_wing_fermion!(xout.w[i5])

        j5=i5+1
        if 1 <= j5 <= xout.L5
            #-P_-
            if xout.L5 != 2
                #mul_1plusγ5x_add!(xout.w[i5],x.w[j5],-1*ratio)
                mul_1minusγ5x_add!(xout.w[i5],x.w[j5],ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

        j5=i5-1
        if 1 <= j5 <= xout.L5
            #-P_+
            if xout.L5 != 2
                #mul_1minusγ5x_add!(xout.w[i5],x.w[j5],-1*ratio)
                mul_1plusγ5x_add!(xout.w[i5],x.w[j5],ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

        if L5 != 1
            if i5==1
                j5 = xout.L5
                #mul_1minusγ5x_add!(xout.w[i5],x.w[j5],m*ratio)
                mul_1plusγ5x_add!(xout.w[i5],x.w[j5],-m*ratio)
                set_wing_fermion!(xout.w[i5])
            end

            if i5==xout.L5
                j5 = 1
                #mul_1plusγ5x_add!(xout.w[i5],x.w[j5],m*ratio)
                mul_1minusγ5x_add!(xout.w[i5],x.w[j5],-m*ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

    end

    #if L5 != xout.L5
    #    for i5=L5+1:xout.L5
    #        axpy!(1,x.w[i5],xout.w[i5])
    #    end
    #end

    if L5 != xout.L5
        for i5 in irange_out
            axpy!(1,x.w[i5],xout.w[i5])
        end
    end

    set_wing_fermion!(xout)

    return
end

=#


text:
module TestOnlineStats

include("../preamble.jl")

using Distributed: @everywhere
@everywhere import OnlineStats

using OnlineStats: Mean, fit!, value
using Transducers: OSNonZeroNObsError

const parseint = Base.Fix1(parse, Int)
const square = Base.Fix2(^, 2)
const lessthan5 = Base.Fix2(<, 5)

@testset "$_fold" for _fold in [foldl, foldxt, foldxd]
    if _fold === foldl
        fold = _fold
    else
        fold(args...; kw...) = _fold(args...; basesize = 1, kw...)
    end
    @testset begin
        @test value(fold(Mean(), Map(identity), 1:9)) == 5
        @test value(fold(Mean(), Map(parseint), "123456789")) == 5
        @test value(fold(Mean(), Map(square), 1:11)) ≈ 46
        @test value(fold(Mean(), Filter(lessthan5), 1:9)) == 2.5
    end
end

@testset "non-zero `nobs`" begin
    m = Mean()
    fit!(m, 1)
    @testset "$fold" for fold in [foldxt, foldxd]
        @test_throws OSNonZeroNObsError fold(m, Map(identity), 2:9)
    end
    @testset "foldl" begin
        @test value(foldl(m, Map(identity), 2:9)) == 5
    end
end

end  # module


text:
"""
This is a test usecase of `CompilerBridge`, which is really not complete.
Here we define `BridgeChecker <: AbstractInterpreter` and overload the `Core.Compiler`
functions defined  within base/compiler/abstractinterpretation.jl
by literally copy-and-pasting their definitions, but within our own module context.
So if the coverage of "`Base.function(::Core.Compiler.Types)` bridge" is not enough,
the `typeinf(::BridgeChecker)` may result in an error somewhere.
"""
module TestBridge

const CC = Core.Compiler

for name in names(Core; all = true)
    isdefined(@__MODULE__, name) && continue
    Core.eval(@__MODULE__, :(import Core: $name))
end

for name in names(Core.Compiler; all = true)
    startswith(string(name), "#") && continue
    isdefined(@__MODULE__, name) && continue
    Core.eval(@__MODULE__, :(import .CC: $name))
end

struct BridgeChecker <: CC.AbstractInterpreter
    native::CC.NativeInterpreter
end
BridgeChecker() = BridgeChecker(NativeInterpreter())
InferenceParams(interp::BridgeChecker) = InferenceParams(interp.native)
OptimizationParams(interp::BridgeChecker) = OptimizationParams(interp.native)
get_world_counter(interp::BridgeChecker) = get_world_counter(interp.native)
get_inference_cache(interp::BridgeChecker) = get_inference_cache(interp.native)
code_cache(interp::BridgeChecker) = code_cache(interp.native)

let
    JULIA_DIR = normpath(Sys.BINDIR, Base.DATAROOTDIR, "julia")
    # JULIA_DIR = normpath(Sys.BINDIR, "..", "..")
    ABSTRACTINTERPRETATION = normpath(JULIA_DIR, "base", "compiler", "abstractinterpretation.jl")
    s = read(ABSTRACTINTERPRETATION, String)
    s = replace(s, "AbstractInterpreter" => "BridgeChecker")
    ex = Base.parse_input_line(s)
    @assert Meta.isexpr(ex, :toplevel)
    for x in ex.args
        isa(x, LineNumberNode) && continue
        occursin("BridgeChecker", string(x)) || continue # skip if this code block doesn't contain an overload
        Core.eval(@__MODULE__, x)
    end
end

using Test
code_typed(sin, (Int,); interp=BridgeChecker())
code_typed(println, (QuoteNode,); interp=BridgeChecker())
@test true

end # module TestBridge

max_stars_repo_name: aviatesk/CompilerBridge

max_issues_repo_name: aviatesk/CompilerBridge

max_forks_repo_name: aviatesk/CompilerBridge


text:
#
# chronologicalAcquisition.jl -
#
# This example demonstrates the usage of the `readChronological` method.  This
# method is always reading the image that is the oldest one in the loop
# buffers.
#
# Please refer to the "simpleAcquisition.jl" example if you're looking for a
# simpler example to acquire images.
#

module ChronologicalAcquisition

# In this example, we want to use the low-level interface as it is the closest
# to the C API.
import NuvuCameras: NC

function main(imageAcquisition::Integer = 40)

    # Opens the acquisition channel using the automatic detection and 4 loop
    # buffers (recommended).
    myCam = NC.open(NC.Cam, NC.AUTO_UNIT, NC.AUTO_CHANNEL, 4)

    # For this example we will not be checking for the available readout modes
    # on the camera and simply assume that the first one is valid.
    NC.setReadoutMode(myCam, 1)

    # Modifies the timeout value to "one second", meaning each image has up to
    # 1 second to reach the frame grabber before the function indicates a
    # timeout.
    NC.setTimeout(myCam, 1000)

    print("Press the Return to start an acquisition sequence of ",
          "$imageAcquisition images;\nas we simulate a long processing time, ",
          "some of the frames in the loop buffers will be overwritten.\n\n",
          "The `readChronological` method will step over these frames to ",
          "read the oldest ones in the loop buffers\n\n")
    readline(STDIN)

    # Open the shutter for the acquisition.
    NC.setShutterMode(myCam, NC.OPEN)

    # Launches 'imageCount' acquisitions by the framegrabber (this function
    # does not wait for the acquisition to be complete before returning).
    NC.start(myCam, imageCount)

    imageCount = imageAcquisition
    i = 0
    while i < imageCount

	# Reads the image received, if a timeout occurs an error code will be
	# returned.
	myImage, nbrImagesSkipped = NC.readChronological(myCam)
        i += 1

	# If you wish to proceed with some processing on each image read, it
	# should be done here (in this example we'll only save the image).

	# Saves each image acquired, at the end of each name the loop index
	# will be added.
	imageName = "Image_$i"

	# Saves the image acquired.
	NC.saveImage(myCam, myImage, imageName, NC.TIF,
                     "This is one of the images grabbed in the chronological example",
                     true)

	println("Number of images skipped to get image $i: ", nbrImagesSkipped)

	# Images that get overwritten are lost from the acquisition sequence
	# Skip that number of lost reads at the end of the sequence.
	imageCount -= nbrImagesSkipped

	# For the purpose of this example we will simulate a processing time
	# that is longer then the acquisition rate (this will depend on the
	# frame rate you currently have set on your camera).  The application
	# will here for a while, allowing for the images coming from the camera
	# to overwrite the previous ones in your loop buffer.
	sleep(0.1)
    end

    println("Total number of images lost from acquisition sequence: ",
            imageAcquisition  - imageCount)

    # Close the shutter, now that the acquistion is complete.
    NC.setShutterMode(myCam, NC.CLOSE)

    # Closes the acquisition channel no longer in use.
    NC.close(myCam)

    if is_windows()
	print("\nHit the Return key to close the command window.")
        flush(STDOUT)
        readline(STDIN)
    else
	println("The program ended.")
    end
end

end # module


text:
using StatsBase: var
function solver_default(
    X,
    data::AbstractArray{T,2};
    stderror = false,
) where {T<:Union{Missing,<:Number}}
    # likely much larger matrix, using lsqr
    minfo = []
    beta = Array{Float64}(undef, size(data, 1), size(X, 2))

    @showprogress 0.1 for ch = 1:size(data, 1)
        ix = .!ismissing.(data[ch, :])
        beta[ch, :], h = lsmr(X[ix, :], data[ch, ix], log = true)
        push!(minfo, h)
    end

    if stderror
        stderror = calculate_stderror(X, data, beta)
        modelfit = LinearModelFit(beta, ["lsmr"], stderror)
    else
        modelfit = LinearModelFit(beta, ["lsmr"])
    end
    return modelfit
end

function calculate_stderror(Xdc, data::Matrix{T}, beta) where {T<:Union{Missing,<:Number}}

    # remove missings
    ix = any(.!ismissing.(data), dims = 1)[1, :]
    if length(ix) != size(data, 2)
        @warn(
            "Limitation: Missing data are calculated over all channels for standard error"
        )
    end
    data = data[:, ix]
    Xdc = Xdc[ix, :]

    # Hat matrix only once
    hat_prime = inv(Matrix(Xdc' * Xdc))
    # Calculate residual variance
    @warn(
        "Autocorrelation was NOT taken into account. Therefore SE are UNRELIABLE. Use at your own discretion"
    )

    se = Array{Float64}(undef, size(data, 1), size(Xdc, 2))
    for ch = 1:size(data, 1)
        residualVar = var(data[ch, :] .- Xdc * beta[ch, :])
        @assert(!isnan(residualVar), "residual Variance was NaN")
        hat = hat_prime .* residualVar
        #se = sqrt(diag(cfg.contrast(:,:)*hat*cfg.contrast(:,:)'));
        se[ch, :] = sqrt.(diag(hat))
    end
    return se
end
function calculate_stderror(
    X,
    data::AbstractArray{T,3},
    beta,
) where {T<:Union{Missing,<:Number}}
    #function calculate_stderror(Xdc,data::AbstractArray{T,2},beta) where {T<:Union{Missing, <:Number}}

    # Hat matrix
    hat_prime = inv(Matrix(X' * X))
    # Calculate residual variance
    @warn(
        "Autocorrelation was NOT taken into account. Therefore SE are UNRELIABLE. Use at your own discretion"
    )

    se = Array{Float64}(undef, size(data, 1), size(data, 2), size(X, 2))
    for ch = 1:size(data, 1)
        for t = 1:size(data, 2)
            ix = .!ismissing.(data[ch, t, :])
            residualVar = var(data[ch, t, ix] .- X[ix, :] * beta[ch, t, :])
            @assert(!isnan(residualVar), "residual Variance was NaN")
            hat = hat_prime .* residualVar
            #se = sqrt(diag(cfg.contrast(:,:)*hat*cfg.contrast(:,:)'));
            se[ch, t, :] = sqrt.(diag(hat))
        end
    end
    return se
end
function solver_default(
    X,
    data::AbstractArray{T,3};
    stderror = false,
) where {T<:Union{Missing,<:Number}}
    beta = Array{Union{Missing,Number}}(undef, size(data, 1), size(data, 2), size(X, 2))
    @showprogress 0.1 for ch = 1:size(data, 1)
        for t = 1:size(data, 2)
            @debug("$(ndims(data,)),$t,$ch")
            ix = .!ismissing.(data[ch, t, :])
            beta[ch, t, :] = X[ix, :] \ data[ch, t, ix]
        end
    end

    if stderror
        stderror = calculate_stderror(X, data, beta)
        modelfit = LinearModelFit(beta, ["solver_default"], stderror)
    else
        modelfit = LinearModelFit(beta, ["solver_default"])
    end
    return modelfit
end

solver_b2b(X, data, cross_val_reps) = solver_b2b(X, data, cross_val_reps = cross_val_reps)
function solver_b2b(
    X,
    data::AbstractArray{T,3};
    cross_val_reps = 10,
) where {T<:Union{Missing,<:Number}}

    X, data = dropMissingEpochs(X, data)

    E = zeros(size(data, 2), size(X, 2), size(X, 2))
    W = Array{Float64}(undef, size(data, 2), size(X, 2), size(data, 1))

    prog = Progress(size(data, 2) * cross_val_reps, 0.1)
    for t = 1:size(data, 2)

        for m = 1:cross_val_reps
            k_ix = collect(Kfold(size(data, 3), 2))
            Y1 = data[:, t, k_ix[1]]
            Y2 = data[:, t, k_ix[2]]
            X1 = X[k_ix[1], :]
            X2 = X[k_ix[2], :]

            G = (Y1' \ X1)
            H = X2 \ (Y2' * G)

            E[t, :, :] = E[t, :, :] + Diagonal(H[diagind(H)])
            ProgressMeter.next!(prog; showvalues = [(:time, t), (:cross_val_rep, m)])
        end
        E[t, :, :] = E[t, :, :] ./ cross_val_reps
        W[t, :, :] = (X * E[t, :, :])' / data[:, t, :]

    end

    # extract diagonal
    beta = mapslices(diag, E, dims = [2, 3])
    # reshape to conform to ch x time x pred
    beta = permutedims(beta, [3 1 2])
    modelinfo = Dict("W" => W, "E" => E, "cross_val_reps" => cross_val_reps) # no history implemented (yet?)
    return LinearModelFit(beta, modelinfo)
end


text:
"""
    envstats(signal,epochlength,overlap,fs; pruned = 0.5)

Calculate the mean and coefficient of variation of the envelope for each epoch.

# Arguments
* `epochlength`: The length of the epoch in seconds
* `overlap`: The overlap between epochs in seconds
* `fs`: The sampling frequency
* `pruned`: The number of seconds cut from both tails of the envelope before calculating the CVE
"""
function envstats(signal,epochlength,overlap,fs; pruned = 0.5)
    signal = arraysplit(signal,epochlength*fs,overlap*fs)
    signal = broadcast.(abs,hilbert.(signal))
    signal = map(x->x[Int(fs*pruned)+1:end-Int(fs*pruned)],signal)
    return (variation.(signal),mean.(signal))
end


text:
module GeneMisc

#using ParallelAccelerator
using JLD

const data_dir = joinpath(dirname(dirname(@__FILE__)), "data")
const hs_fl    = joinpath(data_dir, "Homo_sapiens.gene_info")
const genecode_fl = joinpath(data_dir, "gencode.v19.annotation.gtf")

#index files

### synonym index
const gene_id_fl  = joinpath(data_dir, "gene_id_dict.jld")
const id_genes_fl = joinpath(data_dir, "id_genes_dict.jld")

### location index
const chr_sted_gene_fl = joinpath(data_dir, "chr_sted_gene_dict.jld")
const gene_chrsted_fl  = joinpath(data_dir, "gene_chrsted.jld")

### exon index

global id_genes,gene_id,chr_sted_gene,gene_chrsted

export query_gene, query_geneloc, query_exon

function __init__()
    load_index()
end

include("utils/utils.jl")
include("index/index.jl")
include("query/query.jl")

end # module


text:
@testset "Geo interface" begin
    pt = LibGEOS.Point(1.0,2.0)
    @test GeoInterface.coordinates(pt) ≈ [1,2] atol=1e-5
    @test GeoInterface.geotype(pt) == :Point

    pt = LibGEOS.Point(1, 2)
    @test GeoInterface.coordinates(pt) ≈ [1,2] atol=1e-5
    @test GeoInterface.geotype(pt) == :Point

    pt = LibGEOS.readgeom("POINT EMPTY")
    @test GeoInterface.coordinates(pt) ≈ Float64[] atol=1e-5
    @test GeoInterface.geotype(pt) == :Point

    mpt = LibGEOS.readgeom("MULTIPOINT(0 0, 10 0, 10 10, 11 10)")
    @test GeoInterface.coordinates(mpt) == Vector{Float64}[[0,0],[10,0],[10,10],[11,10]]
    @test GeoInterface.geotype(mpt) == :MultiPoint

    coords = Vector{Float64}[[8,1],[9,1],[9,2],[8,2]]
    ls = LibGEOS.LineString(coords)
    @test GeoInterface.coordinates(ls) == coords
    @test GeoInterface.geotype(ls) == :LineString

    ls = LibGEOS.readgeom("LINESTRING EMPTY")
    @test GeoInterface.coordinates(ls) == []
    @test GeoInterface.geotype(ls) == :LineString

    mls = LibGEOS.readgeom("MULTILINESTRING ((5 0, 10 0), (0 0, 5 0))")
    @test GeoInterface.coordinates(mls) == [[[5,0],[10,0]],[[0,0],[5,0]]]
    @test GeoInterface.geotype(mls) == :MultiLineString

    coords = Vector{Float64}[[8,1],[9,1],[9,2],[8,2],[8,1]]
    lr = LibGEOS.LinearRing(coords)
    @test GeoInterface.coordinates(lr) == coords
    @test GeoInterface.geotype(lr) == :LineString

    coords = Vector{Vector{Float64}}[Vector{Float64}[[0,0],[10,0],[10,10],[0,10],[0,0]],
                                     Vector{Float64}[[1,8],[2,8],[2,9],[1,9],[1,8]],
                                     Vector{Float64}[[8,1],[9,1],[9,2],[8,2],[8,1]]]
    polygon = LibGEOS.Polygon(coords)
    @test GeoInterface.coordinates(polygon) == coords
    @test GeoInterface.geotype(polygon) == :Polygon

    polygon = LibGEOS.readgeom("POLYGON EMPTY")
    @test GeoInterface.coordinates(polygon) == [[]]
    @test GeoInterface.geotype(polygon) == :Polygon

    multipolygon = LibGEOS.readgeom("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))")
    @test GeoInterface.coordinates(multipolygon) == Vector{Vector{Vector{Float64}}}[Vector{Vector{Float64}}[Vector{Float64}[[0,0],[0,10],[10,10],[10,0],[0,0]]]]
    @test GeoInterface.geotype(multipolygon) == :MultiPolygon

    geomcollection = LibGEOS.readgeom("GEOMETRYCOLLECTION (POLYGON ((8 2, 10 10, 8.5 1, 8 2)), POLYGON ((7 8, 10 10, 8 2, 7 8)), POLYGON ((3 8, 10 10, 7 8, 3 8)), POLYGON ((2 2, 8 2, 8.5 1, 2 2)), POLYGON ((2 2, 7 8, 8 2, 2 2)), POLYGON ((2 2, 3 8, 7 8, 2 2)), POLYGON ((0.5 9, 10 10, 3 8, 0.5 9)), POLYGON ((0.5 9, 3 8, 2 2, 0.5 9)), POLYGON ((0 0, 2 2, 8.5 1, 0 0)), POLYGON ((0 0, 0.5 9, 2 2, 0 0)))")
    collection = GeoInterface.geometries(geomcollection)
    coords = Vector{Vector{Vector{Float64}}}[Vector{Vector{Float64}}[Vector{Float64}[[8.0,2.0],[10.0,10.0],[8.5,1.0],[8.0,2.0]]],
                                             Vector{Vector{Float64}}[Vector{Float64}[[7.0,8.0],[10.0,10.0],[8.0,2.0],[7.0,8.0]]],
                                             Vector{Vector{Float64}}[Vector{Float64}[[3.0,8.0],[10.0,10.0],[7.0,8.0],[3.0,8.0]]],
                                             Vector{Vector{Float64}}[Vector{Float64}[[2.0,2.0],[ 8.0, 2.0],[8.5,1.0],[2.0,2.0]]],
                                             Vector{Vector{Float64}}[Vector{Float64}[[2.0,2.0],[ 7.0, 8.0],[8.0,2.0],[2.0,2.0]]],
                                             Vector{Vector{Float64}}[Vector{Float64}[[2.0,2.0],[ 3.0, 8.0],[7.0,8.0],[2.0,2.0]]],
                                             Vector{Vector{Float64}}[Vector{Float64}[[0.5,9.0],[10.0,10.0],[3.0,8.0],[0.5,9.0]]],
                                             Vector{Vector{Float64}}[Vector{Float64}[[0.5,9.0],[ 3.0, 8.0],[2.0,2.0],[0.5,9.0]]],
                                             Vector{Vector{Float64}}[Vector{Float64}[[0.0,0.0],[ 2.0, 2.0],[8.5,1.0],[0.0,0.0]]],
                                             Vector{Vector{Float64}}[Vector{Float64}[[0.0,0.0],[ 0.5, 9.0],[2.0,2.0],[0.0,0.0]]]]
    @test map(GeoInterface.coordinates,collection) == coords
    @test GeoInterface.geotype(geomcollection) == :GeometryCollection

    geomcollection = LibGEOS.readgeom("GEOMETRYCOLLECTION(MULTIPOINT(0 0, 0 0, 1 1),LINESTRING(1 1, 2 2, 2 2, 0 0),POLYGON((5 5, 0 0, 0 2, 2 2, 5 5)))")
    collection = GeoInterface.geometries(geomcollection)
    coords = Vector[Vector{Float64}[[0.0,0.0],[0.0,0.0],[1.0,1.0]],
                    Vector{Float64}[[1.0,1.0],[2.0,2.0],[2.0,2.0],[0.0,0.0]],
                    Vector{Vector{Float64}}[Vector{Float64}[[5.0,5.0],[0.0,0.0],[0.0,2.0],[2.0,2.0],[5.0,5.0]]]]
    geotypes = [:MultiPoint,:LineString,:Polygon]
    for (i,item) in enumerate(collection)
        @test GeoInterface.coordinates(item) == coords[i]
        @test GeoInterface.geotype(item) == geotypes[i]
    end
    @test GeoInterface.geotype(geomcollection) == :GeometryCollection

    geomcollection = LibGEOS.readgeom("GEOMETRYCOLLECTION EMPTY")
    collection = GeoInterface.geometries(geomcollection)
    @test length(collection) == 0
    @test GeoInterface.geotype(geomcollection) == :GeometryCollection
end


text:
## Primitive datatypes
# These get special handling only in that they have different HDF5 type
# representations than ordinary opaque types

# This construction prevents these methods from getting called on type unions
const SignedTypes        = Union{Type{Int8}, Type{Int16}, Type{Int32}, Type{Int64}, Type{Int128}}
const UnsignedTypes      = Union{Type{UInt8}, Type{UInt16}, Type{UInt32}, Type{UInt64}, Type{UInt128}}
const FloatTypes         = Union{Type{Float16}, Type{Float32}, Type{Float64}}
const PrimitiveTypeTypes = Union{SignedTypes, UnsignedTypes, FloatTypes, Type{Bool}}
const PrimitiveTypes     = Union{Bool, Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32,
                                 UInt64, UInt128, Float16, Float32, Float64}

for T in Base.uniontypes(SignedTypes)
    @eval h5fieldtype(::JLDFile, ::$T, ::$T, ::Initialized) =
        FixedPointDatatype($(T.parameters[1].size), true)
end
for T in Base.uniontypes(UnsignedTypes)
    @eval h5fieldtype(::JLDFile, ::$T, ::$T, ::Initialized) =
        FixedPointDatatype($(T.parameters[1].size), false)
end

function jltype(f::JLDFile, dt::FixedPointDatatype)
    signed = dt.bitfield1 == 0x08 ? true : dt.bitfield1 == 0x00 ? false : throw(UnsupportedFeatureException("Invalid Bit field 1"))
    ((dt.bitfield2 == 0x00) & (dt.bitfield3 == 0x00) & (dt.bitoffset == 0) & (dt.bitprecision == dt.size*8)) ||
        throw(UnsupportedFeatureException())
    if dt.size == 8
        return signed ? ReadRepresentation{Int64,Int64}() : ReadRepresentation{UInt64,UInt64}()
    elseif dt.size == 1
        return signed ? ReadRepresentation{Int8,Int8}() : ReadRepresentation{UInt8,UInt8}()
    elseif dt.size == 4
        return signed ? ReadRepresentation{Int32,Int32}() : ReadRepresentation{UInt32,UInt32}()
    elseif dt.size == 2
        return signed ? ReadRepresentation{Int16,Int16}() : ReadRepresentation{UInt16,UInt16}()
    elseif dt.size == 16
        return signed ? ReadRepresentation{Int128,Int128}() : ReadRepresentation{UInt128,UInt128}()
    else
        throw(UnsupportedFeatureException())
    end
end

# Special handling for booleans as they are not considered <: Integer in HDF5
h5fieldtype(::JLDFile, ::Type{Bool}, ::Type{Bool}, ::Initialized) =BitFieldDatatype(1)
jltype(::JLDFile, ::BitFieldDatatype) = ReadRepresentation{Bool, Bool}()

h5fieldtype(::JLDFile, ::Type{Float16}, ::Type{Float16}, ::Initialized) =
    FloatingPointDatatype(DT_FLOATING_POINT, 0x20, 0x0f, 0x00, 2, 0, 16, 10, 5, 0, 10, 0x0000000f)
h5fieldtype(::JLDFile, ::Type{Float32}, ::Type{Float32}, ::Initialized) =
    FloatingPointDatatype(DT_FLOATING_POINT, 0x20, 0x1f, 0x00, 4, 0, 32, 23, 8, 0, 23, 0x0000007f)
h5fieldtype(::JLDFile, ::Type{Float64}, ::Type{Float64}, ::Initialized) =
    FloatingPointDatatype(DT_FLOATING_POINT, 0x20, 0x3f, 0x00, 8, 0, 64, 52, 11, 0, 52, 0x000003ff)
# JLD3.FloatingPointDatatype(0x11, 0x20, 0x3f, 0x00, 0x00000008, 0x0000, 0x0040, 0x34, 0x0b, 0x00, 0x34, 0x000003ff)

function Base.:(==)(fp1::FloatingPointDatatype, fp2::FloatingPointDatatype)
    ret = true
    ret &= (fp1.class << 4) == (fp2.class << 4) # compare only class and not version
    ret &= fp1.bitfield1 == fp2.bitfield1
    ret &= fp1.bitfield2 == fp2.bitfield2
    ret &= fp1.bitfield3 == fp2.bitfield3
    ret &= fp1.size == fp2.size
    ret &= fp1.bitoffset == fp2.bitoffset
    ret &= fp1.bitprecision == fp2.bitprecision
    ret &= fp1.exponentlocation == fp2.exponentlocation
    ret &= fp1.exponentsize == fp2.exponentsize
    ret &= fp1.mantissalocation == fp2.mantissalocation
    ret &= fp1.mantissasize == fp2.mantissasize
    ret &= fp1.exponentbias == fp2.exponentbias
    ret
end

function jltype(f::JLDFile, dt::FloatingPointDatatype)
    if dt == h5fieldtype(f, Float64, Float64, Val{true})
        return ReadRepresentation{Float64,Float64}()
    elseif dt == h5fieldtype(f, Float32, Float32, Val{true})
        return ReadRepresentation{Float32,Float32}()
    elseif dt == h5fieldtype(f, Float16, Float16, Val{true})
        return ReadRepresentation{Float16,Float16}()
    else
        throw(UnsupportedFeatureException())
    end
end

function h5fieldtype(f::JLDFile, writeas::PrimitiveTypeTypes, readas::Type, init::Initialized)
    @lookup_committed f readas
    commit(f, h5fieldtype(f, writeas, writeas, init), writeas, readas)
end
h5type(f::JLDFile, writeas::PrimitiveTypeTypes, x) =
    h5fieldtype(f, writeas, typeof(x), Val{true})

# Used only for custom serialization
constructrr(f::JLDFile, T::PrimitiveTypeTypes, dt::Union{FixedPointDatatype,FloatingPointDatatype},
            ::Vector{ReadAttribute}) =
    dt == h5fieldtype(f, T, T, Val{true}) ? (ReadRepresentation{T,T}(), true) :
                                            throw(UnsupportedFeatureException())


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Styles

const text_formats = Dict(
    :black   => ("\e[30m", "\e[39m"),
    :red     => ("\e[31m", "\e[39m"),
    :green   => ("\e[32m", "\e[39m"),
    :yellow  => ("\e[33m", "\e[39m"),
    :blue    => ("\e[34m", "\e[39m"),
    :magenta => ("\e[35m", "\e[39m"),
    :cyan    => ("\e[36m", "\e[39m"),
    :white   => ("\e[37m", "\e[39m"),
    :reset   => ("\e[0m", "\e[0m"),
    :bold    => ("\e[1m", "\e[22m"),
    :underline => ("\e[4m", "\e[24m"),
    :blink     => ("\e[5m", "\e[25m"),
    :negative  => ("\e[7m", "\e[27m"))

function with_output_format(f::Function, formats::Vector{Symbol}, io::IO, args...)
    Base.have_color && for format in formats
        haskey(text_formats, format) &&
            print(io, text_formats[format][1])
    end
    try f(io, args...)
    finally
        Base.have_color && for format in formats
            haskey(text_formats, format) &&
                print(io, text_formats[format][2])
        end
    end
end

with_output_format(f::Function, format::Symbol, args...) =
    with_output_format(f, [format], args...)

with_output_format(format, f::Function, args...) =
    with_output_format(f, format, args...)

function print_with_format(format, io::IO, x)
    with_output_format(format, io) do io
        print(io, x)
    end
end

function println_with_format(format, io::IO, x)
    print_with_format(format, io, x)
    println(io)
end

# Wrapping

function ansi_length(s)
    replace(s, r"\e\[[0-9]+m", "") |> length
end

words(s) = split(s, " ")
lines(s) = split(s, "\n")

# This could really be more efficient
function wrapped_lines(s::AbstractString; width = 80, i = 0)
    if ismatch(r"\n", s)
        return vcat(map(s->wrapped_lines(s, width = width, i = i), split(s, "\n"))...)
    end
    ws = words(s)
    lines = AbstractString[ws[1]]
    i += ws[1] |> ansi_length
    for word in ws[2:end]
        word_length = ansi_length(word)
        if i + word_length + 1 > width
            i = word_length
            push!(lines, word)
        else
            i += word_length + 1
            lines[end] *= " " * word
        end
    end
    return lines
end

wrapped_lines(f::Function, args...; width = 80, i = 0) =
    wrapped_lines(sprint(f, args...), width = width, i = 0)

function print_wrapped(io::IO, s...; width = 80, pre = "", i = 0)
    lines = wrapped_lines(s..., width = width, i = i)
    println(io, lines[1])
    for line in lines[2:end]
        println(io, pre, line)
    end
    length(lines), length(pre) + ansi_length(lines[end])
end

print_wrapped(f::Function, io::IO, args...; kws...) = print_wrapped(io, f, args...; kws...)

function print_centred(io::IO, s...; columns = 80, width = columns)
    lines = wrapped_lines(s..., width = width)
    for line in lines
        print(io, " "^(div(columns-ansi_length(line), 2)))
        println(io, line)
    end
    length(lines), length(pre) + length(lines[end])
end

function centred(s, columns)
    pad = div(columns - ansi_length(s), 2)
    " "^pad * s
end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
module SpringCollab2020HeatWaveNoColorGrade

using ..Ahorn, Maple

@mapdef Effect "SpringCollab2020/HeatWaveNoColorGrade" HeatWaveNoColorGrade(only::String="*", exclude::String="")

placements = HeatWaveNoColorGrade

function Ahorn.canFgBg(effect::HeatWaveNoColorGrade)
    return true, true
end

end


text:
include("boundconstraint_enum.jl")
include("misc.jl")

include("abstract_bc_quadsolver_conf.jl")
include("abstract_bc_quadsolver_result.jl")
include("bc_quadsolver_interface.jl")

include("Kunisch-Rendl.jl")


text:
module QueensAlwaysMakeProblems

# For some background see: https://wp.me/paipV7-E
# There is also a 'distributed' version:
# https://github.com/PeterLuschny/IntegerSequences.jl/blob/master/demos/DistributedQueens.jl

function solve!(profile, level, size, start, cols, diag4, diag1)

    if level > 0
        for i in start:size
            save = cols & (1 << i) +
            diag1 & (1 << (i + level)) +
            diag4 & (1 << (32 + i - level))

            if save == 0

                cols  = xor(cols,  1 << i)
                diag1 = xor(diag1, 1 << (i + level))
                diag4 = xor(diag4, 1 << (32 + i - level))

                solve!(profile, level - 1, size, 0, cols, diag4, diag1)

                cols  = xor(cols,  1 << i)
                diag1 = xor(diag1, 1 << (i + level))
                diag4 = xor(diag4, 1 << (32 + i - level))

                profile[level + 1] += 1
            end
        end
    else
        for i in 0:size
            save = cols & (1<<i) + diag1 & (1<<i) + diag4 & (1<<(32+i))
            save == 0 && (profile[1] += 1)
        end
    end
end

function search(n::Int)

    profile = zeros(Int, n + 1)
    cols = diag4 = diag1 = Int(0)
    solve!(profile, n - 1, n - 1, 0, cols, diag4, diag1)
    return profile
end

function queens(n::Int)

    n == 0 && return [1]
    profile = search(n)
    profile[n+1] = 1  # add the root
    [profile[n-i+1] for i = 0:n]
end

function demo(up_to)

    for n in 0:up_to
        print("elapsed: ")
        @time profile = queens(n)
        println("size:      ", n)
        println("profile:   ", profile)
        println("nodes:     ", sum(profile))
        println("solutions: ", profile[n+1])
        println()
    end
end

demo(10)

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
#=
remls = [530.1445193510292, -30.67455875307806, 425.4465604691695, 314.22176883261096, -74.87997833266812, 530.1445193182162, 1387.0928273412144, 2342.5993980030553, 2983.26033032097, -16.41729812792036,
250.94514897106058, 1140.3816624784859, 2087.481017283834 , 1012.351698923092, 2087.481017283834 , 323.99766668546584, 77.56902301272578, 904.8743799636109 , 782.9395904949903 ,796.3124436472704,
470.59083255259935,248.99027587947566,119.80621157945501,274.30636135925795,660.0465401679254,433.84147581860896,1123.6556434756412,329.2574937705332,26.96606070210349,26.316526650535426]

remlsb = [536.2011487949176,-28.588171019268458,436.5116765557768,316.05056715993294,-74.54174147034092,536.2011487949178,1388.5367906718716,2345.962740738334,2985.724169788945,-15.89382595425026,
253.0510786936172,1143.3630490411597, 2090.103960238728,1050.3094549379975,2090.103960238728,326.53614303770223,78.84472797596732,937.9855867469028,815.2250801964899,827.9845298841103,
474.08826545614767,249.2537119063981,129.16329544316406,283.59476106417014,678.7453061894956,435.44528875067806,1125.0834323641677,331.8084510486417,45.304914448494124,32.41804784450319]
=#
remlsb = zeros(Float64, 30)

remlsc = zeros(Float64, 30)

mblcic = zeros(Float64, 30)
mbucic = zeros(Float64, 30)
mclcic = zeros(Float64, 30)
mcucic = zeros(Float64, 30)

remlsrb = [536.20114880,-28.58817102,436.51167656,316.05056716,-74.54174147,
           536.20114880,1388.53679067,2345.96274074,2985.72416979,-15.89382595,
           253.05107869,1143.36304904,2090.10396024,1050.30945494,2090.10396024,
           326.53614304,78.84472798,937.98558675,815.22508020,827.98452988,
           474.08826546,249.25371191,129.16329544,283.59476106,678.74530619,
           435.44528875,1125.09238307,331.80845105,45.30491445,32.41804784]

remlsrc= [530.14451859,-30.67456491,425.44656047,314.22176883,-74.87997833,
          530.14451859,1387.09282701,2342.59939800,2983.26033032,-16.41729813,
           250.94514897,1140.38166248,2087.48101397,1012.35169892,2087.48101397,
           323.99766668,77.56902301,904.87437996,782.93959049,796.31244365,
           470.59083255,248.99027588,119.80621158,274.30636136,660.04654017,
           433.84147582,1123.66435096,329.25749377,26.96606070,26.3165]

commentsc1 = [1,3,6,13,15,16,24,25]
commentsc2 = [2,5,7,10,22,30]

c1 =
"The final Hessian matrix is not positive definite although all convergence criteria
are satisfied. The MIXED procedure continues despite this warning. Validity of subsequent
results cannot be ascertained."
c2 =
"Iteration was terminated but convergence has not been achieved. The MIXED procedure
continues despite this warning. Subsequent results produced are based on the last iteration.
Validity of the model fit is uncertain."
#8, 9, 12
#15?
#! 16, 27

#=
SPSS
REML G - CSH, R - DIAG
13 - 2087.4810139727 / 2087.481017533285
15 - / 2087.481017533285 / 2087.4810192819373
16 - 323.99766668073 / 323.9976737508006
29 - / 26.96606070210349
30 - 14.944038069158

REML G - SI, R - SI
4 - 316.050567
27 - 1125.092383
30 - 32.418048

=#
@testset "  RDS Test                                                 " begin
for i = 1:30
    dfrds        = CSV.File(joinpath(path, "csv", "berds", "rds"*string(i)*".csv"), types = Dict(:PK => Float64, :subject => String, :period => String, :sequence => String, :treatment => String )) |> DataFrame
    dropmissing!(dfrds)
    dfrds.lnpk = log.(dfrds.PK)
    @testset "  RDS Test $(i)                                               " begin
        atol=1E-6
        if i ∈ [27] atol = 1E-2 end
        lmm = Metida.LMM(@formula(lnpk~sequence+period+treatment), dfrds;
        random = Metida.VarEffect(Metida.@covstr(1|subject), Metida.SI),
        )
        Metida.fit!(lmm)
        @test lmm.result.reml ≈ remlsrb[i] atol=atol
        remlsb[i] = Metida.m2logreml(lmm)
        mblcic[i] =  coef(lmm)[end] - stderror(lmm)[end]*quantile(TDist(Metida.dof_satter(lmm)[end]), 1.0-0.1/2)
        mbucic[i] =  coef(lmm)[end] + stderror(lmm)[end]*quantile(TDist(Metida.dof_satter(lmm)[end]), 1.0-0.1/2)

        atol=1E-6
        if i ∈ [2, 13, 15] atol = 1E-4 end
        if i ∈ [27, 30] atol = 1E-2 end
        lmm = Metida.LMM(@formula(lnpk~sequence+period+treatment), dfrds;
        random = Metida.VarEffect(Metida.@covstr(treatment|subject), Metida.CSH),
        repeated = Metida.VarEffect(Metida.@covstr(treatment|subject), Metida.DIAG),
        )
        Metida.fit!(lmm)
        @test lmm.result.reml ≈ remlsrc[i] atol=atol
        remlsc[i] = Metida.m2logreml(lmm)
        mclcic[i] =  coef(lmm)[end] - stderror(lmm)[end]*quantile(TDist(Metida.dof_satter(lmm)[end]), 1.0-0.1/2)
        mcucic[i] =  coef(lmm)[end] + stderror(lmm)[end]*quantile(TDist(Metida.dof_satter(lmm)[end]), 1.0-0.1/2)

        @test mclcic[i] ≈ confint(lmm; level =0.9, ddf = :satter)[end][1]
        @test mcucic[i] ≈ confint(lmm; level =0.9, ddf = :satter)[end][2]
    end
end
end
remlsrc[30] = 14.94403807
dftable = DataFrame(n = collect(1:30), a = remlsb, b = remlsrb, c = Vector{Any}(undef, 30), d = remlsc, e = remlsrc, f = Vector{Any}(undef, 30), g = fill!(Vector{String}(undef, 30),""))

for i = 1:30
    if isapprox(dftable.a[i], dftable.b[i]; atol=1E-6)  dftable.c[i] = "OK" else dftable.c[i] = dftable.a[i] - dftable.b[i] end
    if isapprox(dftable.d[i], dftable.e[i]; atol=1E-6)  dftable.f[i] = "OK" else dftable.f[i] = dftable.d[i] - dftable.e[i] end
    if i in commentsc1
        dftable.g[i] = "*"
    end
    if i in commentsc2
        dftable.g[i] = "**"
    end
end

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

dfrds        = CSV.File(joinpath(path, "csv", "berds", "rds"*string(27)*".csv"), types = Dict(:PK => Float64, :subject => String, :period => String, :sequence => String, :treatment => String )) |> DataFrame
dropmissing!(dfrds)
dfrds.lnpk = log.(dfrds.PK)
fm2 = @formula(lnpk~sequence+period+treatment+(1|subject))
mm  = fit(MixedModel, fm2, dfrds, REML=true)

println("")
println("Bioequivalence Reference Datasets - REML - type B")
header = ["RDS" "Metida" " SPSS " "DIFF" ;
          " N " "REML B" "REML B" "    " ]
header = tuple([header[i, :] for i = 1:size(header, 1)]...)
pretty_table(dftable[!, [1,2,3,4]]; header = header, tf = tf_ascii_rounded)
println("")
println("Bioequivalence Reference Datasets - REML - type C")
header = ["RDS" "Metida" " SPSS " "DIFF" "Comm.";
          " N " "REML C" "REML C" "    " "     "]
header = tuple([header[i, :] for i = 1:size(header, 1)]...)
pretty_table(dftable[!,[1,5,6,7,8]]; header = header, tf = tf_ascii_rounded)
println("")
println("*  - ", c1)
println("")
println("** - ", c2)
println("")
println("DataSet 27 - MixedModels.jl result:")
println("")
println(mm)
println("")
println("Bioequivalence Reference Datasets - 90% Confidence intervals")
lci = [107.17, 97.32, 113.31, 117.90, 103.82,  80.02,  86.46, 75.69, 75.69, 96.27,
        80.64, 90.35,  72.87,  69.21,  72.87,  69.54, 115.97, 59.13, 53.85, 50.92,
       111.72, 77.98,  97.13,  87.24,  77.93, 133.51,  78.86, 87.86, 88.43, 79.58]
uci = [124.97, 107.46, 136.73, 159.69, 112.04,  93.31,  92.81,  87.60,  87.60, 107.59,
       100.38, 157.88,  85.51, 121.27,  85.51,  89.37, 155.09, 107.20,  98.77,  95.62,
       127.73, 106.09, 128.41, 109.85,  98.10, 171.42,  89.30, 100.07, 121.59, 108.07]

citable = DataFrame(n = collect(1:30),
a = round.(exp.(mblcic) .* 100.0, digits = 2),
b = round.(exp.(mbucic) .* 100.0, digits = 2),
c = round.(exp.(mclcic) .* 100.0, digits = 2),
d = round.(exp.(mcucic) .* 100.0, digits = 2),
e = lci, f = uci)
header = ["RDS" "Metida B" "Metida B" "Metida C" "Metida C" "REF*" "REF*";
          " N " "LCI"      "UCI"      "LCI"      " UCI "    "LCI" "UCI"]
header = tuple([header[i, :] for i = 1:size(header, 1)]...)
pretty_table(citable; header = header , tf = tf_ascii_rounded)
println("")
println("* Reference: Schütz H, Labes D, Tomashevskiy M, la Parra MG, Shitova A,
Fuglsang A. Reference Datasets for Studies in a Replicate Design Intended for
Average Bioequivalence with Expanding Limits. AAPS J. 2020 Feb 7;22(2):44.
doi: 10.1208/s12248-020-0427-6. PMID: 32034551.
Table II, Method B, SPSS section.")


text:
"""
    mleresults(model, θ)

Do maximum likelihood estimation, where model is a function that gives
the vector of log likelihoods of the observations, and θ is the
parameter vector. call edit(mleresults()) to see a coded example.

"""
function mleresults()
    println("execute edit(mleresults,()) to examine the example code")
    x = [ones(10) 1:10]
    y = [1:2;3;4;5;6;7;8;9;10]
    β = zeros(2)
    model = β -> poisson(β, y, x)
    βhat, junk, junk, junk = mleresults(model, β, "simple ML example")
    βhat
end

function mleresults(model, θ, title="", names=""; vc=1)
    n = size(model(θ),1)
    thetahat, objvalue, V, converged = mle(model, θ, vc)
    k = size(V,1)
    if names==""
        names = 1:k
        names = names'
    end
    se = sqrt.(diag(V))
    t = thetahat ./ se
    p = 2.0 .- 2.0*cdf.(Ref(TDist(n-k)),abs.(t))
    PrintDivider()
    if title !="" printstyled(title, color=:yellow); println() end
    print("MLE Estimation Results    Convergence: ")
    printstyled(converged, color=:green)
    println()
    println("Average log-ℒ: ", round(objvalue; digits=5), "   Observations: ", n)
    if vc==1
        println("Sandwich form covariance estimator")
    elseif vc==2
        println("Inverse Hessian form covariance estimator")
    else
        println("OPG form covariance estimator")
    end
    a =[thetahat se t p]
    println("")
    PrintEstimationResults(a, names)
    println()
    println("Information Criteria")
    caic = -2.0*n*objvalue + k*(log(n)+1.0)
    bic = -2.0*n*objvalue + k*log(n)
    aic = -2.0*n*objvalue + 2.0*k
    infocrit = [caic; bic; aic]
    infocrit = round.([infocrit infocrit/n],digits=5)
    clabels = ["Crit.", "Crit/n"]
    rlabels = ["CAIC ", "BIC ", "AIC "]
    prettyprint(infocrit, clabels, rlabels)
    PrintDivider()
    return thetahat, objvalue, V, converged
end

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
function cumulate_grad!(∇m, m, neigs)
    p = 1.
    @inbounds for si in neigs
        J, i = sign(si), abs(si)
        p *= (1-J*m[i]) / 2
    end
    @inbounds for si in neigs
        J, i = sign(si), abs(si)
        ∇m[i] += J * p/(1-J*m[i]) / (1-p)
        @assert isfinite(∇m[i]) "$(∇m[i]) $(m[i]) p=$p "
    end
end

clip(h::Vector{Float64}) = Int.(sign.(h))

function solve_gd(cnf::CNF;
                    maxiters = 5000, ϵ::Float64 = 1e-4,
                    γ::Float64 = 0.1, γstep::Float64=0.01, γmax = 0.5,
                    η = 0.01, # gd learning rate
                    seed::Int = -1,
                    infotime = 100)

    seed > 0 && srand(seed)
    N, M = cnf.N, cnf.M
    h = randn(N)
    Δ = 0.

    σ = clip(h)
    E = energy(cnf, σ)
    @printf("it=0 γ=%.3f\t  E=%d   \tΔ=%f\n", γ, E, Δ)
    for it=1:maxiters
        m = tanh.(γ .* h)
        ∇m = zeros(N)
        for μ=1:M
            cumulate_grad!(∇m, m, cnf.clauses[μ])
        end

        h .+= γ .* η .* (1 .- m.^2) .* ∇m
        h ./= norm(h)/√N

        σ = clip(h)
        E = energy(cnf, σ)
        # it%infotime == 0 && @printf("\rit=%d γ=%.3f\t  E=%d   \tΔ=%f", it, γ, E, Δ)
        if it % infotime == 0
            @printf("it=%d γ=%.3f\t  E=%d \n", it, γ, E)
            γ += γstep
        end

        E == 0 && break
    end
    @printf("γ=%.3f\t  E=%d\n", γ, E)

    return clip(h)
end


text:
"""This module implements the projection algorithm."""
# Author: Henri Gérard <hgerard.proy@gmail.com>
# License: MIT

function proj_pi(f::Function,
                 t::Array{Float64,1},
                 x::Array{Float64,1},
                 i::Int,
                 data::Array{Float64,2},
                 regressionModel::RegressionModel)
    """
    Apply on vector x a projection step on function f.

    Parameters
    ----------
    f::Function
        The function on which we want to project.
    t::Array{Float64,1}
        Gradient of the function f.
    x::Array{Float64,1}
        Vector we want to project
    i::Int

    data::Array{Float64,2}

    regressionModel::RegressionModel
        Define the regression model used.

    Returns
    -------
    X_leaves : array_like, shape = [n_samples, n_estimators]
        For each datapoint x in X and for each tree in the forest,
        return the index of the leaf x ends up in.
    """
    if f(x) <= 0
        return x
    else
        step = f(x)*t/(t'*t)
        if isnan(sum(step))
            if i > size(data)[1]
                x[1] = max(0,x[1])
                return x
            else
                step = equivalent(x, i, data, regressionModel)
            end
        end
        return x - step
    end
end

# Aggregator of projection
function QR(x0::Array{Float64,1},
            xn::Array{Float64,1},
            zn::Array{Float64,1})
    """
    Aggregator of projection

    Parameters
    ----------
    x0::Array{Float64,1},

    xn::Array{Float64,1},

    zn::Array{Float64,1}

    Returns
    -------
    nothing
    """
    π = (x0-xn)'*(xn-zn)
    μ =(x0-xn)'*(x0-xn)
    ν = (xn-zn)'*(xn-zn)
    ρ = μ*ν-π^2
    if ((ρ==0.0) && (π >=0))
        return zn
    end
    if (ρ>0)
         if (π*ν>=ρ)
            return x0+(1+π/ν)*(zn-xn)
         else
            return xn+ν/ρ*(π*(x0-xn)+μ*(zn-xn))
         end
    end
    if ρ <= 0.0
        return zn
    end
    if isnan(ρ)
        # println("Warning, rho = NaN")
        # println("pi = ", π)
        # println("nu = ", ν)
        # println("mu = ", μ)
        # println("zn = ", zn)
        # println("x0 = ", x0)
        # println("xn = ", xn)
        if (π*ν>=ρ)
           return x0+(1+π/ν)*(zn-xn)
        else
           return xn+ν/ρ*(π*(x0-xn)+μ*(zn-xn))
        end
    end
    println("rho = ", ρ)
    println("error in QR")
    return nothing
end

function typeof_constraint(i::Int, robustModel::RobustModel)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------
    Name of constraint
    """
    if i == length(robustModel.I0)+1
        return PositiveConstraint()
    else
        return robustModel.name
    end
end

#  Test if a point x is inside a set of constraint
function x_in_inter(x,
                     data,
                     robustModel,
                     I,
                     para::ParallelParam)

    error("'error_loss' not defined for $(typedof(s))")
end

function x_in_inter(x::Array{Float64,1},
                    data::Array{Float64,2},
                    robustModel::RobustModel,
                    I::Array{Int64,1},
                    para::Sequential)
     """
     Aggregator of projection

     Parameters
     ----------
     i::Int

     robustModel::RobustModel

     Returns
     -------
     Name of constraint
     """

    aux::Float64 = 0.0
    for k in 1:length(I)
        i::Int = I[k]
        name::GeneralConstraint = typeof_constraint(i, robustModel)
        aux = max(aux,fconstraint(x, data, robustModel.regressionModel, i, name))
    end
    isnan(aux) && (aux = 1.0)
    return aux
end

function x_in_inter(x::Array{Float64,1},
                     data::Array{Float64,2},
                     robustModel::RobustModel,
                     I::Array{Int64,1}, para::Parallel)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------
    Name of constraint
    """
    aux = @distributed (max) for k in 1:length(I)
        max(0,fconstraint(x, data, robustModel.regressionModel, I[k], typeof_constraint(I[k], robustModel)))
    end
    isnan(aux) && (aux = 1.0)
    return aux
end

# Projection onto a set of constraint
function proj_In(xn::Array{Float64,1},
                  data::Array{Float64,2},
                  robustModel::RobustModel,
                  I::Array{Int64,1},
                  ω::Array{Float64,1},
                  para::Sequential)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------
    Name of constraint
    """
    n = length(I)
    pin::Array{Float64,1} = zeros(length(xn))
    sum_op = zeros(length(xn))
    Ln = 0
    for k in 1:n
        i = I[k]
        name = typeof_constraint(i, robustModel)
        subg = subgradient(xn, data, robustModel.regressionModel, i, name)
        pin = proj_pi(x->fconstraint(x, data, robustModel.regressionModel, i, name), subg, xn, i, data, robustModel.regressionModel)
        sum_op = sum_op + ω[k]*pin
        Ln = Ln + ω[k]*((pin-xn)'*(pin-xn))
    end
    un = xn-sum_op
    Ln = Ln/(un'*un)
    return un, Ln
end

# auxilaray function to speed uo proj parallel
function auxf(i::Int,
                xn::Array{Float64,1},
                data::Array{Float64,2},
                robustModel::RobustModel,
                ω::Float64)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------
    Name of constraint
    """

    name = typeof_constraint(i, robustModel)
    subg = subgradient(xn, data, robustModel.regressionModel, i, name)
    pin::Array{Float64,1} = proj_pi(x->fconstraint(x, data, robustModel.regressionModel, i, name), subg, xn, i, data, robustModel.regressionModel)
    aux::Array{Float64,1} =  ω*vcat(pin, dot((pin-xn),(pin-xn)))

    return aux
end

function proj_In(xn::Array{Float64,1},
                  data::Array{Float64,2},
                  robustModel::RobustModel,
                  I::Array{Int64,1},
                  ω::Array{Float64,1},
                  para::Parallel)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------
    Name of constraint
    """
    n::Int = length(I)
    ul::Array{Float64,1} = @distributed (+) for k in 1:n
        auxf(I[k], xn, data, robustModel, ω[k])
    end
    un::Array{Float64,1} = xn - ul[1:end-1]
    Ln::Float64 = ul[end]
    Ln = Ln/dot(un, un)
    return un, Ln
end

# Projection algorithm
function algo_proj(x0::Array{Float64,1},
                    data::Array{Float64,2},
                    robustModel::RobustModel,
                    projParams::ProjParams)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------
    Name of constraint
    """
    dist_mem::Float64 = projParams.precision+1
    xn::Array{Float64,1} = copy(x0)
    xn1::Array{Float64,1} = copy(xn)+ones(length(xn))
    iter = 0
    while (iter < Int(round(length(robustModel.I0)/projParams.sample))+1) || ((dist_mem > projParams.precision) && (iter < projParams.ITER_MAX))
        # println("iter = ", iter, " dm = ", dist_mem)
        iter = iter +1
        I = StatsBase.sample(robustModel.I0, projParams.sample, replace=false, ordered=true)
        push!(I, length(robustModel.I0)+1)
        ω = ones(length(I))/(length(I))
        un, Ln = proj_In(xn, data, robustModel, I, ω, projParams.para_proj)
        dist_mem = x_in_inter(xn, data, robustModel, I, projParams.para_inter)
        if (dist_mem <= 0.0) || (dot(un,un) == 0)
            Ln = 1
            # dist_mem = (xn1-xn)'*(xn1-xn)
            I = StatsBase.sample(robustModel.I0, projParams.sample, replace=false, ordered=true)
            dist_mem = x_in_inter(xn, data, robustModel, I, projParams.para_inter)
        end
        zn = xn - Ln*un
        xn1 = xn
        xn = QR(x0,xn,zn)
    end
    return xn, dist_mem, iter
end

function aux_proj(i::Int,
                    xn::Array{Float64,1},
                    data::Array{Float64,2},
                    robustModel::RobustModel,
                    ω::Float64)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------
    Name of constraint
    """
    name = typeof_constraint(i, robustModel)
    subg::Array{Float64,1} = subgradient(xn, data, robustModel.regressionModel, i, name)
    pin::Array{Float64,1} = proj_pi(x->fconstraint(x, data, robustModel.regressionModel, i, name), subg, xn, i, data, robustModel.regressionModel)
    aux::Float64 = max(0, fconstraint(xn, data, robustModel.regressionModel, i, name))
    return vcat(ω*pin, ω*dot((pin-xn),(pin-xn)), aux)
end

function algo_proj_bis(x0::Array{Float64,1},
                    data::Array{Float64,2},
                    robustModel::RobustModel,
                    projParams::ProjParams)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------
    Name of constraint
    """
    dist_mem::Float64 = projParams.precision+1
    xn::Array{Float64,1} = copy(x0)
    iter::Int = 0
    while (iter < Int(round(length(robustModel.I0)/projParams.sample))+1) || ((dist_mem > projParams.precision) && (iter < projParams.ITER_MAX))
        # println("iter = ", iter, " dm = ", dist_mem)
        iter = iter +1
        I = StatsBase.sample(robustModel.I0, projParams.sample, replace=false, ordered=true)
        push!(I, length(robustModel.I0)+1)
        ω = ones(length(I))/(length(I))
        res::Array{Float64,1} = @distributed (+) for k in 1:length(I)
            aux_proj(I[k], xn, data, robustModel, ω[k])
        end
        dist_mem = res[end]/length(I)
        un::Array{Float64,1} = xn-res[1:end-2]
        Ln::Float64 = res[end-1]/(un'*un)
        if (dist_mem <= 0.0) || (dot(un,un) == 0)
            Ln = 1
            # dist_mem = (xn1-xn)'*(xn1-xn)
            I = StatsBase.sample(robustModel.I0, projParams.sample, replace=false, ordered=true)
            dist_mem = x_in_inter(xn, data, robustModel, I, projParams.para_inter)
        end
        xn = QR(x0,xn,xn - Ln*un)
    end
    return xn, dist_mem, iter
end

function run_algo_bis(xn::Array{Float64,1},
                    data::Array{Float64,2},
                    robustModel::RobustModel,
                    optParams::OptParams,
                    projParams::ProjParams)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------
    Name of constraint
    """
    # dist_mem = optParams.stability+1
    dm::Array{Float64,1} = []
    mini::Float64 = Inf
    mem::Array{Float64,1} = zeros(size(xn))
    y::Array{Float64,1} = zeros(size(xn))
    xn1::Array{Float64,1} = copy(xn)
    iter::Int = 0
    while (iter < optParams.itmax)
        iter = iter +1
        #print(iter," ")
        # xn = xn1 - optParams.learning_rate*robustModel.descent_direction
        y = xn + ((iter-2)/(iter+1))*(xn - xn1)
        xn, DM, it = algo_proj_bis(y-optParams.learning_rate*robustModel.descent_direction,
                                data,
                                robustModel,
                                projParams)
        (robustModel.descent_direction'*xn < mini) && (mini = min(mini,robustModel.descent_direction'*xn); mem = xn )
        push!(dm, robustModel.descent_direction'*xn)
        xn1 = xn
        ((iter % optParams.verbosity) == 0) && print(iter," ",DM," ")
    end
    # println("Exit after ", iter, " iterations")
    return xn, y, dm, mem, mini
end

# optimization algorithm
function run_algo(xn::Array{Float64,1},
                    data::Array{Float64,2},
                    robustModel::RobustModel,
                    optParams::OptParams,
                    projParams::ProjParams)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------
    Name of constraint
    """
    # dist_mem = optParams.stability+1
    dm::Array{Float64,1} = []
    t_iter::Array{Float64,1} = []
    mini::Float64 = Inf
    mem::Array{Float64,1} = zeros(size(xn))
    y::Array{Float64,1} = zeros(size(xn))
    xn1::Array{Float64,1} = copy(xn)
    iter::Int = 0
    # while (iter < optParams.itmax) || (iter > 2 &&  (dm[end]-dm[end-1]<optParams.stability))
    while (iter < optParams.itmax)
        iter = iter +1
        #print(iter," ")
        # xn = xn1 - optParams.learning_rate*robustModel.descent_direction
        y = xn + ((iter-2)/(iter+1))*(xn - xn1)
        start = time()
        xn, DM, it = algo_proj(y-optParams.learning_rate*robustModel.descent_direction,
                                data,
                                robustModel,
                                projParams)
        time_iter = time() - start
        (robustModel.descent_direction'*xn < mini) && (mini = min(mini,robustModel.descent_direction'*xn); mem = xn )
        push!(dm, robustModel.descent_direction'*xn)
        push!(t_iter, time_iter)
        xn1 = xn
        ((iter % optParams.verbosity) == 0) && print(iter," ",DM," ")
    end
    # println("Exit after ", iter, " iterations")
    return xn, y, dm, mem, mini, t_iter
end

max_stars_repo_name: Henri-Gerard/robox

max_issues_repo_name: Henri-Gerard/robox

max_forks_repo_name: Henri-Gerard/robox


text:
raw_data = OrderedDict{Symbol, DataFrame}()
low = DataFrame(
    class = [:negcontrol, :not, :not, :not],
    counts=[10, 100, 1000, 0],
    barcodeid=[1,2,3,4],
    gene=[1,2,2,3],
    behavior=[:linear,:linear,:linear,:linear]
)
high = DataFrame(
    class = [:not, :negcontrol, :not, :not],
    counts=[10, 100, 1000, 10],
    barcodeid=[2,1,3,4],
    gene=[2,1,2,3],
    behavior=[:linear,:linear,:linear,:linear]
)
raw_data[:low] = low
raw_data[:high] = high

guide_data, gene_data = differences_between_bins(raw_data)

total = sum([10, 100, 1000, 0] .+ 0.5)
@test guide_data[:freqs_low][1] == 10.5/total
@test guide_data[:freqs_low][4] == 0.5/total
@test guide_data[:rel_freqs_low][1] == 1.0 # only a single negative control should it should be divided by itself
@test guide_data[:rel_freqs_low][2] == 100.5/total/(10.5/total)

# make sure the log2 fold changes are correct
@test all(guide_data[:log2fc_high_div_low] .== log2.(guide_data[:rel_freqs_high]./guide_data[:rel_freqs_low]))

# > 2 bins
raw_data = OrderedDict{Symbol, DataFrame}()
raw_data[:low] = low

raw_data[:mid] = DataFrame(
    class = [:not, :negcontrol, :not, :not],
    counts=[10, 100, 1000, 100],
    barcodeid=[2,1,3,4],
    gene=[2,1,2,3],
    behavior=[:linear,:linear,:linear,:linear]
)

raw_data[:high] = high
guide_data, gene_data = differences_between_bins(raw_data)

# make sure all the log2fc comparisons are there
@test all(map(x->contains(==, names(guide_data), x), [:log2fc_mid_div_low, :log2fc_high_div_low, :log2fc_high_div_mid]))


text:
#!/usr/bin/julia --color=yes

# Usage: call it with .net files as ARGS.

function parse_file(file)

    io = open(file,"r")
    println(file,":")

    header = read(io,Int8,8)
    @assert join(map(Char,header),"") == "ANNEALER"
    IS_BIG_ENDIAN = ENDIAN_BOM == 0x01020304

    is_big_endian = read(io,Int64)
    @assert is_big_endian == IS_BIG_ENDIAN

    L = read(io,Int64)
    V = L^3

    nmeas = read(io,Int64)
    tws = read(io,Int64,nmeas)

    nts = read(io,Int64)
    ts = read(io,Int64,nts)

    total_meas = read(io, Int64)

    T = read(io, Float64)

    Jup    = read(io,UInt64,V)
    Jright = read(io,UInt64,V)
    Jfront = read(io,UInt64,V)

    println("    L             : ", L)
    print("    tw's          : ")
    @printf("∈ [%.4g, %.4g] in %d steps\n",
            minimum(tws), maximum(tws), length(tws))
    print("    t's           : ")
    @printf("∈ [%.4g, %.4g] in %d steps\n",
            minimum(ts), maximum(ts), length(ts))
    println("    Endianness    : ", is_big_endian==1? "Big Endian":"Little Endian")
    println("    T             : ", T)
    println("    J hash        : ",
            hex(hash(Jup) + hash(Jright) + hash(Jfront)))

    mcs = Int64[]

    for i in 1:total_meas
        try
            mc    = read(io,Int64)
            spins = read(io,UInt64,V)
            push!(mcs,mc)
        catch e
            if e == EOFError()
                warn("EOF was reached.")
                print()
                print_with_color(:red, "Latest mc's: ")
                print_with_color(:red, "$(mcs[end-0]), ")
                print_with_color(:red, "$(mcs[end-1]), ")
                print_with_color(:red, "$(mcs[end-2]), ")
                print_with_color(:red, "$(mcs[end-3])\n")
                break
            else
                throw(e)
            end
        end
    end

    mcdigest =
    @sprintf("mc steps ∈ (%.4g,%.4g), %.4g values of %.4g\n",
             minimum(mcs), maximum(mcs), length(mcs), total_meas)

    print_with_color(:blue, mcdigest)

    close(io)
end

for file in ARGS
    parse_file(file)
end

max_stars_repo_name: RedPointyJackson/tfg

max_issues_repo_name: RedPointyJackson/tfg

max_forks_repo_name: RedPointyJackson/tfg


text:
#
# Truncate to Support Sampler
#
# ensures that the sampled value is within defined support
# only DistributionSamplers are supported
#
# We use the Truncated functionality of the Distributions package, but this currently does not
# handle all underlying distributions.
# We leave it to the Distributions to raise an error when it is not compatible: the advantage is that
# additional distributions will be automatically handled should the Distributions package handle them
# in future versions.
#

type TruncateToSupportSampler <: ModifyingSampler
	subsampler::DistributionSampler
	function TruncateToSupportSampler(subsampler::DistributionSampler)
		new(subsampler)
	end
end

function sample(s::TruncateToSupportSampler, support, cc::ChoiceContext)
	# TODO: check whether it is better to store truncated distribution, and only
	# recreate when support changes.  If so, will need to trap setparams call so
	# as to clear any stored truncated distribution
	if support[1] == support[2]
		# can give an error when support[1] == support[2]
		x = support[1]
	else
		truncateddistribution = Truncated(s.subsampler.distribution, support[1], support[2])
		x = rand(truncateddistribution)
	end
	x, Dict{Symbol, Any}(:sub=>Dict{Symbol, Any}(:rnd=>x))
end


text:
using OpenCL
cl = OpenCL
function cl_platforms()
  for pf in cl.platforms()
    @printf "\n%s\n\n" pf
    @printf "Platform name:\t\t%s\n"  pf[:name]
    if pf[:name] == "Portable Computing Language"
      warn("PCL platform is not yet supported")
      continue
    else
      @printf "Platform profile\t\t:%s\n"  pf[:profile]
      @printf "Platform vendor:\t\t%s\n"  pf[:vendor]
      @printf "Platform version:\t\t%s\n\n" pf[:version]

      for dv in cl.available_devices(pf)
        @printf "Device name:\t\t%s\n" dv[:name]
        @printf "Device type:\t\t%s\n" dv[:device_type]
        gms = dv[:global_mem_size] / (1024*1024)
        @printf "Device memory:\t%i MB\n"  gms
        mma = dv[:max_mem_alloc_size] / (1024*1024)
        @printf "Device max memory alloc:\t%i MB\n" mma
        mcu = device[:max_compute_units]
        @printf "Device max compute units:\t%i\n"  mcu
        mwgs = device[:max_work_group_size]
        @printf "Device max work group size:\t%i\n"mwgs
        mwis = device[:max_work_item_size]
        @printf "Device max work item size:\t%s\n"  mwis
      end
    end
  end
end

max_stars_repo_path: Module 3/Chapter10/cl-platforms.jl

max_issues_repo_path: Module 3/Chapter10/cl-platforms.jl

max_forks_repo_path: Module 3/Chapter10/cl-platforms.jl


text:
using Base: splat

test = """
    2199943210
    3987894921
    9856789892
    8767896789
    9899965678
    """

function parse_input(input::AbstractString)::Matrix
    return [[parse(Int, c) for c in row]
            for row in split(input, r"\s+"; keepempty = false)] |>
                splat(hcat) |> permutedims
end

function q1(seafloor::Matrix)::Int
    sum(getindex.(Ref(seafloor), low_points(seafloor)) .+ 1)
end

function low_points(seafloor::Matrix)::Vector
    return findall(keys(seafloor)) do i
        return minimum(neighbors(seafloor, i)) == seafloor[i]
    end
end

neighbors(seafloor::AbstractMatrix, ci) = neighbors(seafloor, ci.I...)
function neighbors(seafloor::M, r::Int, c::Int)::M where {M <: AbstractMatrix}
    return seafloor[(r - 1:r + 1) ∩ (1:end),
                    (c - 1:c + 1) ∩ (1:end)]
end

function q2(seafloor::Matrix)::Int
    basins = basin_size.(Ref(seafloor), low_points(seafloor))
    return sort!(basins; rev = true)[1:3] |> prod
end

function basin_size(seafloor::Matrix, ci::CartesianIndex)::Int
    basin = [ci]
    for ci in basin
        for cj in map(f -> f(ci), [up, down, left, right])
            peak(seafloor, cj) || union!(basin, [cj])
        end
    end
    return length(basin)
end

up(ci::CartesianIndex) = CartesianIndex(ci.I .- (1, 0))
down(ci::CartesianIndex) = CartesianIndex(ci.I .+ (1, 0))
left(ci::CartesianIndex) = CartesianIndex(ci.I .- (0, 1))
right(ci::CartesianIndex) = CartesianIndex(ci.I .+ (0, 1))
peak(m::Matrix, ci::CartesianIndex) = ci ∉ keys(m) || m[ci] == 9

let seafloor = parse_input(test)
    @assert q1(seafloor) == 15
    @assert q2(seafloor) == 1134
end

seafloor = parse_input(read("day09.in", String))
println("Q1: ", q1(seafloor))
println("Q2: ", q2(seafloor))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
using HTTP.MessageRequest: body_was_streamed

struct DownloadsBackend <: AWS.AbstractBackend
    downloader::Union{Nothing,Downloads.Downloader}
end

DownloadsBackend() = DownloadsBackend(nothing)

const AWS_DOWNLOADER = Ref{Union{Nothing,Downloader}}(nothing)
const AWS_DOWNLOAD_LOCK = ReentrantLock()

# Here we mimic Download.jl's own setup for using a global downloader.
# We do this to have our own downloader (separate from Downloads.jl's global downloader)
# because we add a hook to avoid redirects in order to try to match the HTTPBackend's
# implementation, and we don't want to mutate the global downloader from Downloads.jl.
# https://github.com/JuliaLang/Downloads.jl/blob/84e948c02b8a0625552a764bf90f7d2ee97c949c/src/Downloads.jl#L293-L301
function get_downloader(downloader=nothing)
    lock(AWS_DOWNLOAD_LOCK) do
        yield() # let other downloads finish
        downloader isa Downloader && return nothing
        while true
            downloader = AWS_DOWNLOADER[]
            downloader isa Downloader && return nothing
            D = Downloader()
            D.easy_hook =
                (easy, info) -> Curl.setopt(easy, Curl.CURLOPT_FOLLOWLOCATION, false)
            AWS_DOWNLOADER[] = D
        end
    end
    return downloader
end

# https://github.com/JuliaWeb/HTTP.jl/blob/2a03ca76376162ffc3423ba7f15bd6d966edff9b/src/MessageRequest.jl#L84-L85
body_length(x::AbstractVector{UInt8}) = length(x)
body_length(x::AbstractString) = sizeof(x)

read_body(x::IOBuffer) = take!(x)
function read_body(x::IO)
    close(x)
    return read(x)
end

function _http_request(backend::DownloadsBackend, request::Request, response_stream::IO)
    # HTTP.jl sets this header automatically.
    request.headers["Content-Length"] = string(body_length(request.content))

    # We pass an `input` only when we have content we wish to send.
    input = !isempty(request.content) ? IOBuffer(request.content) : nothing

    downloader = @something(backend.downloader, get_downloader())

    # set the hook so that we don't follow redirects. Only
    # need to do this on per-request downloaders, because we
    # set our global one with this hook already.
    if backend.downloader !== nothing && downloader.easy_hook === nothing
        downloader.easy_hook =
            (easy, info) -> Curl.setopt(easy, Curl.CURLOPT_FOLLOWLOCATION, false)
    end

    local buffer
    local response
    try
        @repeat 4 try
            # Use a sacrificial I/O stream so that we only write the `response_stream` once
            # even with multiple attempts.
            buffer = Base.BufferStream()

            # Rewind the input on each attempt otherwise every subsequent attempt will send an
            # empty payload.
            input !== nothing && seekstart(input)

            r = @mock Downloads.request(
                request.url;
                input=input,
                # Compatibility with Downloads.jl versions below v1.5.2
                # See: https://github.com/JuliaLang/Downloads.jl/issues/131
                output=request.request_method != "HEAD" ? buffer : nothing,
                method=request.request_method,
                headers=request.headers,
                verbose=false,
                throw=true,
                downloader=downloader,
            )

            response = _http_response(request, r; throw=true)
        catch e
            @delay_retry if (
                (isa(e, HTTP.StatusError) && AWS._http_status(e) >= 500) ||
                isa(e, Downloads.RequestError)
            )
            end
        end
    finally
        close(buffer)

        # Transfer the contents of the `BufferStream` into `response_stream` variable.
        write(response_stream, buffer)
    end

    return AWS.Response(response, response_stream)
end

function _http_response(req::Request, res::Downloads.Response; throw::Bool=true)
    response = HTTP.Response(
        res.status, res.headers; body=body_was_streamed, request=nothing
    )

    if throw && HTTP.iserror(response)
        target = HTTP.resource(HTTP.URI(req.url))
        e = HTTP.StatusError(res.status, req.request_method, target, response)
        Base.throw(e)
    end

    return response
end