text:
module TexasData

using CSV, JSON
using Dates
using DataFrames

projectbasepath = "../../"
include(joinpath(projectbasepath, "src/processing/LocalDataCommon.jl"))
include(joinpath(projectbasepath, "src/processing/GeographicData.jl"))

function load_data_tx(
		date_range,
		los_dist, pct_beds_available,
		focus_locations, locations_limit;
		use_rounding=false,
	)
	@assert focus_locations isa Array || focus_locations isa Int

	rounddata = use_rounding ? x -> round.(Int, x) : identity

	rawdata = DataFrame(CSV.File(joinpath(projectbasepath, "data/local/texas/tx_tsa_hospitalizations.csv")))

	locations = sort(unique(rawdata.tsa_name))
	locations_limit = isnothing(locations_limit) ? length(locations) : locations_limit
	if focus_locations isa Int
		focus_locations = locations[1:focus_locations]
	end
	locations = vcat(focus_locations, setdiff(locations, focus_locations))
	locations = locations[1:min(length(locations), max(locations_limit, length(focus_locations)))]
	tsa_ids = [filter(row -> row.tsa_name == l, rawdata).tsa_id[1] for l in locations]

	data_common = (
		date_range = date_range,
		start_date = minimum(date_range),
		end_date = maximum(date_range),
		T = length(date_range),
		locations = locations,
		tsa_ids = tsa_ids,
		N = length(locations),
		focus_locations = focus_locations,
	)

	adj = compute_adj_tx(data_common)
	data_common = merge(data_common, (adj = adj,))

	data_dict = Dict((row.tsa_name, row.date) => (
		active_total          = row.patients_covid_total,
		active_icu            = row.patients_covid_icu,
		active_ward           = row.patients_covid_ward,
		capacity_allpat_total = row.beds_avail_total,
		capacity_allpat_icu   = row.beds_avail_icu,
		capacity_allpat_ward  = row.beds_avail_total - row.beds_avail_icu,
		allpat_total          = row.patients_all_total,
		allpat_icu            = row.patients_all_icu,
		allpat_ward           = row.patients_all_total - row.patients_all_icu,
		capacity_covid_total  = row.beds_avail_total - (row.patients_all_total - row.patients_covid_total),
		capacity_covid_icu    = row.beds_avail_icu - (row.patients_all_icu - row.patients_covid_icu),
		capacity_covid_ward   = (row.beds_avail_total - row.beds_avail_icu) - ((row.patients_all_total - row.patients_all_icu) - row.patients_covid_ward),
	) for row in eachrow(rawdata))
	start_date = minimum(date_range)

	data_allpat = merge(data_common, (
		active   = rounddata([data_dict[(c,d)].active_total for c in locations, d in date_range]),
		initial  = rounddata([data_dict[(c,start_date-Day(1))].active_total for c in locations]),
		capacity = rounddata([data_dict[(c,d)].capacity_covid_total for c in locations, d in date_range]),
		los_dist = los_dist.allpat,
		pct_beds_available = pct_beds_available.allpat,
	))
	adm_dis_all = LocalDataCommon.estimate_admitted_discharged(data_allpat, use_rounding=use_rounding)
	admitted_uncertainty_all = LocalDataCommon.estimate_admitted_uncertainty(data_allpat, adm_dis_all.admitted)
	beds_all = load_beds_tx(data_allpat, bedtype=:all, use_rounding=use_rounding)
	data_allpat = merge(data_allpat, adm_dis_all, (beds = beds_all, admitted_uncertainty = admitted_uncertainty_all))

	data_icu = merge(data_common, (
		active   = rounddata([data_dict[(c,d)].active_icu for c in locations, d in date_range]),
		initial  = rounddata([data_dict[(c,start_date-Day(1))].active_icu for c in locations]),
		capacity = rounddata([data_dict[(c,d)].capacity_covid_icu for c in locations, d in date_range]),
		los_dist = los_dist.icu,
		pct_beds_available = pct_beds_available.icu,
	))
	adm_dis_icu = LocalDataCommon.estimate_admitted_discharged(data_icu, use_rounding=use_rounding)
	admitted_uncertainty_icu = LocalDataCommon.estimate_admitted_uncertainty(data_icu, adm_dis_icu.admitted)
	beds_icu = load_beds_tx(data_icu, bedtype=:icu, use_rounding=use_rounding)
	data_icu = merge(data_icu, adm_dis_icu, (beds = beds_icu, admitted_uncertainty = admitted_uncertainty_icu))

	data_ward = merge(data_common, (
		active   = rounddata([data_dict[(c,d)].active_ward for c in locations, d in date_range]),
		initial  = rounddata([data_dict[(c,start_date-Day(1))].active_ward for c in locations]),
		capacity = rounddata([data_dict[(c,d)].capacity_covid_ward for c in locations, d in date_range]),
		los_dist = los_dist.ward,
		pct_beds_available = pct_beds_available.ward,
	))
	adm_dis_ward = LocalDataCommon.estimate_admitted_discharged(data_ward, use_rounding=use_rounding)
	admitted_uncertainty_ward = LocalDataCommon.estimate_admitted_uncertainty(data_ward, adm_dis_ward.admitted)
	beds_ward = load_beds_tx(data_ward, bedtype=:ward, use_rounding=use_rounding)
	data_ward = merge(data_ward, adm_dis_ward, (beds = beds_ward, admitted_uncertainty = admitted_uncertainty_ward))

	data = merge(data_common, (
		total = data_allpat,
		icu = data_icu,
		ward = data_ward,
	))

	data_carepathmodel = LocalDataCommon.build_carepaths_data(data)
	data = merge(data, (
		carepaths = data_carepathmodel,
	))

	return data
end;

function load_beds_tx(data; bedtype=:all, use_rounding=true)
	counties_list = JSON.parsefile(normpath(projectbasepath, "rawdata/texas/tsa_counties.json"));
	county_fips = [c["county_fips"] for c in counties_list]
	county_fips = parse.(Int, county_fips)

	beds_data = DataFrame(CSV.File(joinpath(projectbasepath, "data/hospitals/definitivehc.csv")))
	allowed_hospital_types = ["Critical Access Hospital", "Long Term Acute Care Hospital", "Short Term Acute Care Hospital"]
	filter!(row -> row.HQ_STATE == "TX", beds_data)
	filter!(row -> row.HOSPITAL_TYPE in allowed_hospital_types, beds_data)
	dropmissing!(beds_data, :FIPS)

	if bedtype == :all
		beds_col = :NUM_STAFFED_BEDS
	elseif bedtype == :icu
		beds_col = :NUM_ICU_BEDS
	elseif bedtype == :ward
		beds_col = :NUM_NON_ICU_BEDS
		beds_data.NUM_NON_ICU_BEDS = beds_data.NUM_STAFFED_BEDS - coalesce.(beds_data.NUM_ICU_BEDS, 0)
	end

	beds_data_agg = combine(groupby(beds_data, :FIPS), beds_col => (x -> sum(skipmissing(x))) => :beds)
	beds_data_dict = Dict(row.FIPS => row.beds for row in eachrow(beds_data_agg))
	beds_by_county = [(haskey(beds_data_dict, c) && !ismissing(beds_data_dict[c])) ? beds_data_dict[c] : 0 for c in county_fips]

	beds_by_tsa = Dict(zip(data.tsa_ids, zeros(Float64, data.N)))
	for (i, c) in enumerate(counties_list)
		tsa_id = c["tsa_id"]
		if haskey(beds_by_tsa, tsa_id)
			beds_by_tsa[tsa_id] += beds_by_county[i]
		end
	end

	beds = [haskey(beds_by_tsa, tsa_id) ? beds_by_tsa[tsa_id] : 0 for tsa_id in data.tsa_ids]
	beds = beds .* data.pct_beds_available

	if use_rounding
		beds = round.(Int, beds)
	end

	return beds
end;

function compute_adj_tx(data)
	# counties_list = JSON.parsefile(normpath(projectbasepath, "data/local/texas/tsa_counties.json"));
	# county_fips = [c["county_fips"] for c in counties_list]
	# county_fips = parse.(Int, county_fips)

	# adj = GeographicData.adjacencies(positions, level=:other, threshold=4.0)
	# return adj

	adj = GeographicData.adjacencies(data.locations, source=:fullyconnected, level=:other, threshold=4.0)
	return adj
end;

end;

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation


text:
@doc """
Match a set of electrodes to those provided

usage: lf, valid = match_sensors(electrodes, sensor_labels)
"""
function match_sensors(sens::Array{S}, lbls::Array{AS}) where {AS<:AbstractString,S<:Sensor}

    valid_idx = Int[]
    for label in lbls
        matched_idx = something(findfirst(isequal(label), labels(sens)), 0)
        if matched_idx != 0
            push!(valid_idx, matched_idx)
        end
        @debug(
            "Label $label matched to $( matched_idx == 0 ? "!! nothing !!" : sens[matched_idx].label)"
        )
    end

    sens = sens[valid_idx]

    return sens, valid_idx
end

function match_sensors(
    lf::Array,
    lf_labels::Array{S},
    labels::Array{S},
) where {S<:AbstractString}
    # Match the sensors in a leadfield array to those provided
    #
    # usage: lf, valid = match_sensors(leadfield, leadfield_labels, sensor_labels)

    valid_idx = Int[]
    for label in labels
        matched_idx = something(findfirst(isequal(label), lf_labels), 0)
        if matched_idx != 0
            push!(valid_idx, matched_idx)
        end
        @debug(
            "Label $label matched to $( matched_idx == 0 ? "!! nothing !!" : lf_labels[matched_idx])"
        )
    end

    @info(
        "Leadfield had $(length(lf_labels)) channels, now has $(length(valid_idx)) channels"
    )

    lf = lf[:, :, valid_idx]
    lf_labels = lf_labels[valid_idx]

    return lf, valid_idx
end

ext: jl


text:
#since the GLM package uses the Distribution types from the Distibutions package, we use these packages not to simulate from their existing functions but to
#use them as type dispatchers for the simulate_glm_trait function

# This super type of all response distribution types
abstract type ResponseDistribution end

#SIMULATE STUDENT-T distributed TRAITS needs to be edited
#T distribution with v degrees of freedom
struct TResponse <: ResponseDistribution
df::Int64
scale::Float64
end

function simulate_glm_trait(μ, dist::TResponse)
if dist.scale < 0
      error("Scale parameter cannot be negative for a T distribution!")
    end
 return t_deviate.(μ, dist.scale, dist.df)
end

#SIMULATE WEIBULL TRAITS
#weibull distributed reponse with shape
struct WeibullResponse <: ResponseDistribution
shape::Float64
end

function simulate_glm_trait(μ, dist::WeibullResponse)
if dist.shape < 0
      error("Shape parameter cannot be negative for a Weibull distribution!")
    end
 return weibull_deviate.(μ, dist.shape)
end

#SIMULATE POISSON traits
struct PoissonResponse <: ResponseDistribution
end

function simulate_glm_trait(μ, dist::PoissonResponse)
if any(μ_i -> μ_i < 0, μ)
      error("Location parameter cannot be negative for a Poisson distribution!")
    end
 return poisson_deviate.(μ)
end

#SIMULATE NORMAL TRAITS
# Normal response type with standard deviation σ
struct NormalResponse <: ResponseDistribution
scale::Float64
end

function simulate_glm_trait(μ, dist::NormalResponse)
if dist.scale <= 0
      error("Scale cannot be negative for a normal distribution!")
    end
 return normal_deviate.(μ, dist.scale) # the sigma is coming from dist.scale "ResponseType" and not additionally from ~N(0,1)
end

#SIMULATE Bernoulli TRAITS
struct BernoulliResponse <: ResponseDistribution
end

function simulate_glm_trait(μ, dist::BernoulliResponse)
 return bernoulli_deviate.(μ)
end

#SIMULATE BINOMIAL TRAITS
struct BinomialResponse <: ResponseDistribution
trials::Int64
end

function simulate_glm_trait(μ, dist::BinomialResponse)
if dist.trials < 0
      error("Trials cannot be negative for a binomial distribution!")
    end
 return binomial_deviate.(μ, dist.trials)
end

#SIMULATE GAMMA TRAITS
#simulate from gamma response with shape parameter
struct GammaResponse <: ResponseDistribution
shape::Float64
function GammaResponse(shape::Float64)
  if(shape <= 0)
    error("shape must be greater than zero in a Gamma distribution!")
  else
    new(shape)
  end
end
end

function simulate_glm_trait(μ, dist::GammaResponse)
if any(μ_i -> μ_i < 0, μ)
      error("Shape cannot be negative for a gamma distribution!")
    end
 return gamma_deviate.(dist.shape, μ)
end

# Inverse Gaussian with shape parameter
struct InverseGaussianResponse <: ResponseDistribution
shape::Float64
function InverseGaussianResponse(shape)
  if(shape <= 0)
    error("Shape parameter must be greater than zero in a inverse Gaussian distribution!")
  else
    return(new(shape))
  end
end
end

#********SIMULATE INVERSE GAUSSIAN TRAITS
function simulate_glm_trait(μ, dist::InverseGaussianResponse)
if dist.shape < 0
      error("Shape must be positive for an inverse Gaussian distribution!")
    end
 return inverse_gaussian_deviate.(dist.shape, μ)
end

#SIMULATE EXPONENTIAL TRAITS ## EDIT THIS GUY
#simulate from exponential response with shape parameter
struct ExponentialResponse <: ResponseDistribution
# scale::Float64
# function ExponentialResponse(scale::Float64)
#   if(scale <= 0)
#     error("scale must be greater than zero in a Exponential distribution!")
#   else
#     new(scale)
#   end
# end
end

function simulate_glm_trait(μ, dist::ExponentialResponse)
  if any(μ_i -> μ_i < 0, μ)
      error("Rate cannot be negative for an Exponential Distribution!")
    end
  return exponential_deviate.(μ)
end

#########BERNOULLI
"""Generates a Bernoulli random deviate with success probability p."""
#********SIMULATE BERNOULLI TRAITS

function bernoulli_deviate(p)
  if rand(eltype(p)) <= p
    return 1
  else
    return 0
  end
end

#########BINOMIAL
"""Generates a binomial random deviate with success probability p
and n trials."""

function binomial_deviate(p, n::Int)
  sucesses = 0
  for i = 1:n
    if rand(eltype(p)) <= p
      sucesses = sucesses + 1
    end
  end
  return sucesses
end

"""Generates an exponential random deviate with mean mu."""

function exponential_deviate(mu)
  return -log(rand()) / mu
end

"""Generates a gamma deviate with shape parameter alpha
and intensity lambda."""

function gamma_deviate(alpha, lambda)
  n = floor(Int, alpha)
  T = eltype(alpha)
  z = - log(prod(rand(T, n)))
  beta = alpha - n
  if beta <= one(T) / 10^6
    y = zero(T)
  elseif beta < one(T) / 10^3
    y = (beta / rand(T))^(one(T) / (one(T) - beta))
  else
    (r, s) = (one(T) / beta, beta - one(T))
    for i = 1:1000
      u = rand(T, 2)
      y = - log(one(T) - u[1]^r)
      if u[2] <= (y / (one(T) - exp(- y)))^s
        exit
      end
    end
  end
  return (z + y) / lambda
end

"""Generates an inverse Gaussian deviate with mean mu and
scale lambda."""

function inverse_gaussian_deviate(lambda, mu)
  w = mu * randn()^2
  c = mu / (2 * lambda)
  x = mu + c * (w - sqrt(w * (4 * lambda + w)))
  p = mu / (mu + x)
  if rand() < p
    return x
  else
    return mu^2 / x
  end
end

"""Generates a normal random deviate with mean mu and
standard deviation sigma."""

function normal_deviate(mu, sigma)
  return sigma * randn() + mu
end

"""Generates a Poisson random deviate with mean mu."""

function poisson_deviate(mu)
  (x, p, k) = (exp(-mu), one(eltype(mu)), 0)
  while p > x
    k = k + 1
    p = p * rand()
  end
  return k - 1
end

"""Generates a t random deviate with, location mu, scale sigma,
and degrees of freedom df."""

function t_deviate(mu, scale, df)
  x = randn(mu)
  y = 2 * gamma_deviate(df / 2, one(mu))
  return scale * (x / sqrt(y / df)) + mu
end

"""Generates a Weibull random deviate with scale lambda
and shape alpha."""

function weibull_deviate(lambda, alpha)
  return lambda * (- log(rand()))^(one(lambda) / alpha)
end

max_stars_repo_name: sarah-ji/Tutorials

max_issues_repo_name: sarah-ji/Tutorials

max_forks_repo_name: sarah-ji/Tutorials


text:
using Test
using Printf
using Logging
using Oceananigans
using LESbrary

Logging.global_logger(OceananigansLogger())

architectures = [CPU()]

max_stars_repo_name: ali-ramadhan/LESbrary

max_issues_repo_name: ali-ramadhan/LESbrary

max_forks_repo_name: ali-ramadhan/LESbrary


text:
using Distributions
using SDE
using Base.Test
srand(8)
# sum([ norm(mean([ sample(MvWiener(d),  linspace(0.,2.,5)).X[1,end] for i in 1:n])) for j in 1:1000] .< r*sqrt(2/n))
#tests with alpha 0.01

n = 1000
d= 2

#0.99 quantile vor mean of sum of squares of d-dimensional normals
#r = sqrt(quantile(Distributions.Chisq(d), 0.99))
r = sqrt(9.211) # for d = 2

#  call 'sample(t,MvWiener(d)).X[:,i]' n times
function Wn(d, t, n, i)
 wn = zeros(d, n)
 for j in 1:n
   wn[:, j] = sample(t, Wiener(d)).yy[:,i]
 end
 wn
end

@test size(sample(linspace(3.,4.,2), Wiener(d)).yy) == (d,2)

# the mean and variance of a Brownian motion at t=2 is 0 and 2
@test norm(mean([ sample(linspace(0.,2.,5),Wiener(d)).yy[:,end] for i in 1:n])) < r*sqrt(2/n) #scale with std(W_2)/sqrt(n)
#chiupper =  quantile(Distributions.Chisq(n), 0.995) #upper 0.005 percentile
#chilower = quantile(Distributions.Chisq(n), 0.005) #lower 0.005 percentile

chiupper = 1119 #n = 1000, 0.995
chilower = 888.5

@test chiupper > n*var( Wn(1, linspace(0.,2.,5), n, 5))/2 > chilower

# check that W(2) has the right covariance matrix
@test (d-> norm(cov( Wn(d, linspace(0.,2.,5), n, 5)') - diagm(2ones(d)))*sqrt(n)/sqrt(d))(50) < 5. # did not figure that out exactly, should fail less then 99 %
@test (d-> norm(cov( Wn(d,[0.,0.1, 0.3, 0.5, 2.], n, 5)') - diagm(2ones(d)))*sqrt(n)/sqrt(d))(50) < 5. # did not figure that out exactly, should fail less then 99 %

function Bn(d, u, v, t, n, i)
 wn = zeros(d, n)
 for j in 1:n
   wn[:, j] = samplebridge(t, u, t[end], v, Wiener(d)).yy[:,i]
 end
 wn
end

# a "deterministic" bridge with only start end endpoint
@test (global t1 = norm(samplebridge( linspace(1.,4.,2), [1.,2.], 4., [5.,7.], Wiener(2)).yy - [1. 5.; 2. 7.])) < eps(10.)

# Check that the bridge has the right mean and Variance

# Covariance of a Brownian bridge from t_1 to t_2 (t_2-t)(s-t_1)/(t_2-t_1)
# here (3-0.5)*(0.5-0)/(3) = 0.4166666666666667

#cq = quantile(Distributions.Chisq(2), 0.99)
qchisq = 9.21034037197618
@test norm(mean(Bn(2,[1.,2.], [5.,7.],[1.,1.1, 1.3, 1.5, 3.], n, 4),2) -  [1.,2.] - ([5.,7.] - [1.,2.])*.5/2) < sqrt(qchisq)*sqrt(0.416/n)

 # (3-0.5)*(0.5-0.1)/(2.9) = 0.3448275862068966
@test  chilower < n*var(Bn(1,[1.], [5.],[0.1, 0.3, 0.5, 3.], n, 3),2 )[1]/0.345  <  chiupper


text:
module musicXML
using EzXML

function generateMusicXML()

# ヘッダー
    doc = XMLDocument()
    root = ElementNode("score-partwise")
    root["version"] = "3.0"
    setroot!(doc, root)

# ノード定義
    partList = ElementNode("part-list")
    scorePart = ElementNode("score-part")
    partName = ElementNode("part-name")
    part = ElementNode("part")
    measure = ElementNode("measure")
    attributes = ElementNode("attributes")
    divisions = ElementNode("divisions")
    time = ElementNode("time")
    beats = ElementNode("beats")
    beatType = ElementNode("beat-type")
    clef = ElementNode("clef")
    sign = ElementNode("sign")
    line = ElementNode("line")
    note = ElementNode("note")
# タイトル
    addelement!(root, "movement-title", "trieMusics")
# パートリスト
    link!(root, partList)
    partName.content = "piano"
    link!(partList, scorePart)
    scorePart["id"] = "piano"
    prettyprint(doc)
    link!(scorePart, partName)
    link!(partName, TextNode("piano"))
    link!(root, part)
    part["id"] = "piano"

# パート開始
    link!(part, measure)
# 小節開始
    measure["number"] = "1"
    link!(measure, attributes)
# 属性開始
    link!(attributes, divisions)
    link!(divisions, TextNode("10080"))
    link!(attributes, time)
    link!(time, beats)
    link!(beats, TextNode("4"))
    link!(time, beatType)
    link!(beatType, TextNode("4"))
    link!(attributes, clef)
    link!(clef, sign)
    link!(sign, TextNode("G"))
    link!(clef, line)
    link!(line, TextNode("2"))

# 音符開始
    targetNote = makeNote("C")
    link!(measure, targetNote)

    for pitch in ["D", "E", "F"]
        nextNote = makeNote(pitch)
        println("nandesyo", nextNote)
        linknext!(targetNote, nextNote)
        targetNote = nextNote
    end

    prettyprint(doc)
    write("result.xml", doc)

end

function makeNote(pitch)
    noteNode =  ElementNode("note")
    pitchNode = ElementNode("pitch")
    stepNode = ElementNode("step")
    octaveNode = ElementNode("octave")
    durationNode = ElementNode("duration")
    typeNode = ElementNode("type")
    link!(noteNode, pitchNode)
    link!(pitchNode, stepNode)
    link!(stepNode, TextNode(pitch))
    link!(pitchNode, octaveNode)
    link!(octaveNode, TextNode("4"))
    link!(noteNode, durationNode)
    link!(durationNode, TextNode("10080"))
    link!(noteNode, typeNode)
    link!(typeNode, TextNode("quarter"))
    return noteNode
end

end


text:
using FactCheck
using SimpleHttpIO

include(joinpath(dirname(@__FILE__), "..", "src", "Chunk.jl"))

using Chunk
const C = Chunk

facts("Test chunk.jl") do
    context("types") do
        cio = IOChunked(IOBuffer())
        @fact isa(cio, AbstractIOSocket) --> true
    end

    data = String["For God so loved the world, that he gave his ",
            "only begotten Son, that whosoever ",
            "believeth in him should not perish, but ",
            "have everlasting life.",
            "(John 3:16)"]

    context("writer") do
        io = IOChunked(IOBuffer())

        for d in data
            C.write(io, d)
        end

        C.done(io)

        seekstart(io.io)

        @fact UTF8String(io.io.data) --> string("2d\r\nFor God so loved the world, that he gave his ",
            "\r\n22\r\nonly begotten Son, that whosoever \r\n28\r\nbelieveth in him should not ",
            "perish, but \r\n16\r\nhave everlasting life.\r\nb\r\n(John 3:16)\r\n0\r\n\r\n")
    end

    context("reader") do
        io = IOChunked(IOBuffer())

        for d in data
            C.write(io, d)
        end

        C.done(io)
        r = C.iterator(io)

        seekstart(io.io)
        d = String[UTF8String(d) for d in collect(r)]
        @fact d --> data
    end
end


text:
## Test discretize given edges
function discretize_test()
	@testset "Discretize" begin
        n = 55;
        inV = collect(range(1, 10, length = n));
        edgeV = [0.0, 6.0, 11.0];
        outV = discretize(inV, edgeV);
		@test all(outV[inV .< edgeV[1]] .== 0)
		@test all(outV[inV .> edgeV[end]] .== 0)

		for ix = 1 : n
			iCl = outV[ix];
			if iCl > 0
				@test inV[ix] >= edgeV[iCl]  &&  inV[ix] <= edgeV[iCl+1]
			else
				@test inV[ix] <= edgeV[1]  ||  inV[ix] >= edgeV[end]
			end
		end
	end
end

## Test bin edges (unweighted)
function bin_edges_test()
    @testset "Bin edges" begin
        n = 450;
        inV = collect(range(-4.0, 3.0, length = n));
        pctV = [0.25, 0.4, 0.85];
        edgeV = bin_edges_from_percentiles(inV, pctV);
        @test length(edgeV) == (length(pctV) + 1)

        # test that fractions match target
        cntV = count_bins(inV, edgeV);
        @test sum(cntV) == round(n * pctV[end])
        fracV = cntV ./ n;
        tgFracV = diff(vcat(0.0, pctV));
        @test isapprox(fracV, tgFracV, atol = 2/n)
    end
end

## Test bin edges (weighted)
function bin_edges_weighted_test()
    @testset "Bin edges weighted" begin
        n = 450;
        inV = collect(range(-4.0, 3.0, length = n));
        # FrequencyWeights must be integers
        wtV = collect(range(0.2, 0.8, length = n));
        pctV = [0.25, 0.4, 0.85];
        edgeV = bin_edges_from_percentiles(inV, wtV, pctV);
        @test length(edgeV) == (length(pctV) + 1)

        # test that fractions below edges match target
        massV = count_bins(inV, wtV, edgeV);
        @test sum(massV) ≈ sum(wtV[inV .<= edgeV[end]])
        fracV = massV ./ sum(wtV);
        tgFracV = diff(vcat(0.0, pctV));
        @test isapprox(fracV, tgFracV, atol = 2/n)
    end
end

## Discretize given percentiles
function discretize_given_percentiles_test()
	@testset "discretize_given_percentiles" begin
        n = 38;
        inV = collect(range(1.0, 5.0, length = n));
        pctV = [0.25, 0.4, 0.85, 1.0];
        classV = discretize_given_percentiles(inV, pctV);
		@test all(classV .>= 1);
		@test all(classV .<= length(pctV));
	end
end

## Discretize (weighted)
function discretize_given_pct_weighted_test()
    @testset "discretize pct weighted" begin
        n = 38;
        inV = collect(range(1.0, 5.0, length = n));
        wtV = collect(range(0.1, 0.5, length = n));
        pctV = [0.25, 0.4, 0.85, 1.0];
        classV = discretize_given_percentiles(inV, wtV, pctV);
	    @test all(classV .>= 1);
	    @test all(classV .<= length(pctV));
    end
end

@testset "Discretize" begin
   discretize_test()
   discretize_given_pct_weighted_test();
   bin_edges_test();
   bin_edges_weighted_test();
   discretize_given_percentiles_test()
end

# ----------------


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    isdefined(FilesystemDatastructures, Symbol("#13#14")) && precompile(Tuple{getfield(FilesystemDatastructures, Symbol("#13#14")),SizeConstrainedFileCache})
    isdefined(FilesystemDatastructures, Symbol("#15#17")) && precompile(Tuple{getfield(FilesystemDatastructures, Symbol("#15#17")),SizeConstrainedFileCache})
    isdefined(FilesystemDatastructures, Symbol("#19#21")) && precompile(Tuple{getfield(FilesystemDatastructures, Symbol("#19#21")),SizeConstrainedFileCache})
end


text:
# final scheme. adaptive kernel.

import Printf
import PyPlot
import Random
import Optim

using LinearAlgebra

using FFTW

import Statistics

import Distributions
#import HCubature
import Interpolations
import SpecialFunctions

import SignalTools

import Utilities # https://gitlab.com/RoyCCWang/utilities

import VisualizationTools # https://gitlab.com/RoyCCWang/visualizationtools
import ProbabilityUtilities

import Convex
import SCS

import Calculus
import ForwardDiff

import RiemannianOptim

include("../src/misc/declarations.jl")

include("../src/RKHS/RKHS.jl")
include("../src/RKHS/kernel.jl")
include("../src/RKHS/interpolators.jl")

include("../src/warp_map/Rieszanalysis.jl")
include("../src/misc/utilities.jl")

include("../src/optimization/sdp.jl")

include("../src/misc/front_end.jl")

PyPlot.close("all")
fig_num = 1

Random.seed!(25)

D = 2
# https://github.com/rdeits/RegionTrees.jl/blob/master/examples/demo/demo.ipynb

import Images, TestImages

rgb_image = TestImages.testimage("lighthouse")
gray_image = Images.Gray.(rgb_image)
#Images.mosaicview(rgb_image, gray_image; nrow = 1)

imA = convert(Matrix{Float64}, gray_image)

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.imshow(imA, cmap = "gray" )

PyPlot.title("oracle image")
PyPlot.legend()

down_factor = 5
imB = imA[1:down_factor:size(imA,1), 1:down_factor:size(imA,2)]

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.imshow(imB, cmap = "gray" )

PyPlot.title("input image")
PyPlot.legend()

function assemblepatchs(patch_inds::Matrix{Tuple{UnitRange{Int},UnitRange{Int}}},
    X::Matrix{Vector{T}},
    y::Matrix{T}) where T

    X_set = Matrix{Matrix{Vector{T}}}(undef, size(patch_inds))
    y_set = Matrix{Matrix{T}}(undef, size(patch_inds))

    for r = 1:size(patch_inds,1)
        for c = 1:size(patch_inds,2)

            inds = patch_inds[r,c]
            # println("inds = ", inds)
            # println("X[inds[1], inds[2]] = ", X[inds[1], inds[2]])

            X_set[r,c] = X[inds[1], inds[2]]
            y_set[r,c] = y[inds[1], inds[2]]
        end
    end

    return X_set, y_set
end

function getpatches(Y::Matrix{T}, patch_sz::Tuple{Int,Int}; half_overlap::Int = 2) where T

    N_rows, N_cols = size(Y)
    #V = collect( Iterators.product(1:1:N_rows, 1:1:N_cols) )

    M_rows = ceil(Int, size(Y,1)/patch_sz[1])
    M_cols = ceil(Int, size(Y,2)/patch_sz[2])

    patch_inds = Matrix{Tuple{UnitRange{Int},UnitRange{Int}}}(undef, M_rows, M_cols)
    for r = 1:M_rows
        for c = 1:M_cols

            r_st = max(1, 1+(r-1)*patch_sz[1] - half_overlap)
            r_fin = min(N_rows, r*patch_sz[1] + half_overlap)

            c_st = max(1, 1+(c-1)*patch_sz[1] - half_overlap)
            c_fin = min(N_cols, c*patch_sz[1] + half_overlap)

            if r == M_rows
                r_fin = N_rows
            end

            if c == M_cols
                c_fin = N_cols
            end

            #patch_inds[r,c] = V[r_st:r_fin, c_st:c_fin]
            patch_inds[r,c] = (r_st:r_fin, c_st:c_fin)
        end
    end

    return patch_inds
end

function fitRKHSpathces(imB::Matrix{T}, θ, σ², patch_inds) where T

    N_rows, N_cols = size(imB)

    # tall.
    x1 = 1#one(T)
    x2 = N_cols/N_rows

    # wide.
    if N_cols > N_rows
        x1 = N_rows/N_cols
        x2 = 1#one(T)
    end

    X0 = collect( Iterators.product(LinRange(0, x1, N_rows), LinRange(0, x2, N_cols)) )
    X = reshape(collect( collect(X0[n]) for n = 1:length(X0) ), size(X0))

    X_set, y_set = assemblepatchs(patch_inds, X, imB)

    η_set = Matrix{RKHSProblemType}(undef, size(patch_inds))

    for r = 1:size(patch_inds,1)
        for c = 1:size(patch_inds,2)
            η_set[r,c] = RKHSProblemType( zeros(T, length(X_set[r,c])),
                vec(X_set[r,c]), θ, σ²)

            fitRKHS!(η_set[r,c], vec(y_set[r,c]))
        end
    end

    return η_set, X_set, y_set
end

patch_sz = (20,20)
patch_inds = getpatches(imB, patch_sz)

σ² = 1e-5
θ = Spline34KernelType(0.2)

println("Training RKHS.")
@time η_set, X_set, y_set = fitRKHSpathces(imB, θ, σ², patch_inds)

function querypatches(N_rows::Int, N_cols::Int, upfactor::Int,
    η_set, patch_inds) where T

    for r = 1:size(η_set,1)
        for c = 1:size(η_set,2)

            x1_st = X_set[r,c][1,1]
            x1_fin = X_set[r,c][end,end]
            N_x1 = size(X_set[r,c],1)

            x2_st = X_set[r,c][1,1]
            x2_fin = X_set[r,c][end,end]
            N_x2 = size(X_set[r,c],2)

            Xq0 = collect( Iterators.product(LinRange(x1_st, x1_fin, N_x1*upfactor),
            LinRange(x2_st, x2_fin, N_x2*upfactor)) )

            Xq = collect( collect(Xq0[n]) for n = 1:length(Xq0) )

            yq = Vector{T}(undef, length(Xq))
            query!(yq, Xq, η_set[r,c])

            out[r,c] = reshape(yq, N_x1*upfactor, N_x2*upfactor)
        end
    end

    return η_set, X_set, y_set
end

# query.
upfactor = 2
Xq0 = collect( Iterators.product(LinRange(0, x1, upfactor*N_rows), LinRange(0, x2, upfactor*N_cols)) )
Xq = collect( collect(Xq0[n]) for n = 1:length(Xq0) )

yq = Vector{Float64}(undef, length(Xq))
println("Querying RKHS.")
@time query!(yq, Xq, η)

imC = reshape(yq, (upfactor*N_rows,upfactor*N_cols))

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.imshow(imC, cmap = "gray" )

PyPlot.title("output image")
PyPlot.legend()

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
function collision_condition(u, t, integrator)
    θ₁ = u[1]; θ₂ = u[3]
    if θ₁>-0.05
        return 1.0
    else
        return θ₂ + 2*θ₁

    end
end

function collision_affect!(integrator)
    terminate!(integrator)
end

collision_cb = ContinuousCallback(collision_condition, collision_affect!)


text:
export DirectlySpecifiedMotion, BasicDirectMotion

abstract type DirectlySpecifiedMotion end

struct BasicDirectMotion <: DirectlySpecifiedMotion

end


text:
###########################################################
#
#   Maximal cliques of undirected graph
#
##########################################################

function maximal_cliques{V}(g::AbstractGraph{V})
    """
    Finds all maximal cliques of an undirected graph.

    julia> using Graphs
    julia> g = simple_adjlist(3, is_directed=false)
    julia> add_edge!(g, 1, 2)
    julia> add_edge!(g, 2, 3)
    julia> maximal_cliques(g)
    2-element Array{Array{Int64,N},1}:
     [2,3]
     [2,1]

    Adapted from the NetworkX find_cliques function: http://networkx.lanl.gov/reference/generated/networkx.algorithms.clique.find_cliques.html#networkx.algorithms.clique.find_cliques
    """

    @graph_requires g vertex_list adjacency_list

    if is_directed(g)
        throw(ArgumentError("graph must be undirected."))
    end

    # Cache nbrs and find first pivot (highest degree)
    maxconn = -1
    nnbrs = Dict{V,Set{V}}()
    pivotnbrs = Set{V}() # handle empty graph
    pivotdonenbrs = Set{V}()  # initialize

    for n in vertices(g)
        nbrs = Set{V}()
        union!(nbrs, out_neighbors(n, g))
        delete!(nbrs, n) # ignore edges between n and itself
        conn = length(nbrs)
        if conn > maxconn
            pivotnbrs = nbrs
            nnbrs[n] = pivotnbrs
            maxconn = conn
        else
            nnbrs[n] = nbrs
        end
    end

    # Initial setup
    cand = Set{V}()
    union!(cand, keys(nnbrs))
    smallcand = setdiff(cand, pivotnbrs)
    done = Set{V}()
    stack = @compat(Tuple{Set{V}, Set{V}, Set{V}})[]
    clique_so_far = V[]
    cliques = Array{V}[]

    # Start main loop
    while !isempty(smallcand) || !isempty(stack)
        if !isempty(smallcand) # Any nodes left to check?
            n = pop!(smallcand)
        else
            # back out clique_so_far
            cand, done, smallcand = pop!(stack)
            pop!(clique_so_far)
            continue
        end
        # Add next node to clique
        push!(clique_so_far, n)
        delete!(cand, n)
        push!(done, n)
        nn = nnbrs[n]
        new_cand = intersect(cand, nn)
        new_done = intersect(done, nn)
        # check if we have more to search
        if isempty(new_cand)
            if isempty(new_done)
                # Found a clique!
                push!(cliques, collect(clique_so_far))
            end
            pop!(clique_so_far)
            continue
        end
        # Shortcut--only one node left!
        if isempty(new_done) && length(new_cand) == 1
            push!(cliques, cat(1, clique_so_far, collect(new_cand)))
            pop!(clique_so_far)
            continue
        end
        # find pivot node (max connected in cand)
        # look in done nodes first
        numb_cand = length(new_cand)
        maxconndone = -1
        for n in new_done
            cn = intersect(new_cand, nnbrs[n])
            conn = length(cn)
            if conn > maxconndone
                pivotdonenbrs = cn
                maxconndone = conn
                if maxconndone == numb_cand
                    break
                end
            end
        end
        # Shortcut--this part of tree already searched
        if maxconndone == numb_cand
            pop!(clique_so_far)
            continue
        end
        # still finding pivot node
        # look in cand nodes second
        maxconn = -1
        for n in new_cand
            cn = intersect(new_cand, nnbrs[n])
            conn = length(cn)
            if conn > maxconn
                pivotnbrs = cn
                maxconn = conn
                if maxconn == numb_cand - 1
                    break
                end
            end
        end
        # pivot node is max connected in cand from done or cand
        if maxconndone > maxconn
            pivotnbrs = pivotdonenbrs
        end
        # save search status for later backout
        push!(stack, (cand, done, smallcand))
        cand = new_cand
        done = new_done
        smallcand = setdiff(cand, pivotnbrs)
    end
    cliques
end


text:
"""
`module JAC.Nuclear`
    ... a submodel of JAC that contains procedures for defining a Nuclear.Model and for calculating various nuclear
        potentials.
"""
module Nuclear

    using  QuadGK, ..Basics,  ..Defaults, ..Radial, ..Math
    export Model

    """
    `struct  Nuclear.Model`  ... defines a type for the nuclear model, i.e. for its form and parameters.

        + Z        ::Float64         ... nuclear charge
        + model    ::String          ... identifier of the nuclear model: {"Fermi", "Point", "Uniform"}
        + mass     ::Float64         ... atomic mass
        + radius   ::Float64         ... (root-mean square) radius of a uniform or Fermi-distributed nucleus
        + spinI    ::AngularJ64      ... nuclear spin I, must be >= 0
        + mu       ::Float64         ... magnetic dipole moment in Bohr magnetons
        + Q        ::Float64         ... electric qadrupole moment
    """
    struct  Model
        Z          ::Float64
        model      ::String
        mass       ::Float64
        radius     ::Float64
        spinI      ::AngularJ64
        mu         ::Float64
        Q          ::Float64
    end

    """
    `Nuclear.Model(Z::Real)`
        ... to specify a Fermi-type nucleus with charge Z, and where the nuclear spin and nuclear moments are all set to zero.
    """
    function Model(Z::Real)
        Z < 0.1  &&  error("Z must be >= 0.1")
        model    = "Fermi"
        mass     = 2*Z + 0.005*Z^2
        radius   = Nuclear.Rrms(mass)
        spinI    = AngularJ64(0)
        mu       = 0.
        Q        = 0.

        Model(Z, model, mass, radius, spinI, mu, Q)
    end

    """
    `Nuclear.Model(Z::Real, model::String)`
        ... to specify a nucleus with charge Z, model = {"Fermi", "point", "uniform"}, and where the nuclear spin and
            nuclear moments are all set to zero.
    """
    function Model(Z::Real, model::String)
        Z < 0.1                                    &&  error("Z must be >= 0.1")
        !(model in ["Fermi", "point", "uniform"])  && error("Inappropriate model = $model")
        if  model == "point"
            mass     = 1000*Z
            radius   = 0.
        else
            mass     = 2*Z + 0.005*Z^2
            radius   = Nuclear.Rrms(mass)
        end
        spinI    = AngularJ64(0)
        mu       = 0.
        Q        = 0.

        Model(Z, model, mass, radius, spinI, mu, Q)
    end

    """
    `Nuclear.Model(nm::Nuclear.Model;`

                Z=..,         model=..,         mass=..,        radius=..,
                spinI=..,     mu=..,            Q=..)
        ... constructor for re-defining a nuclear model nm::Nuclear.Model.
    """
    function Model(nm::Nuclear.Model;            Z::Union{Nothing,Float64}=nothing,          model::Union{Nothing,String}=nothing,
        mass::Union{Nothing,Float64}=nothing,    radius::Union{Nothing,Float64}=nothing,     spinI::Union{Nothing,AngularJ64}=nothing,
        mu::Union{Nothing,Float64}=nothing,      Q::Union{Nothing,Float64}=nothing)

        if  Z         == nothing   Zx          = nm.Z           else   Zx          = Z          end
        if  model     == nothing   modelx      = nm.model       else   modelx      = model      end
        if  mass      == nothing   massx       = nm.mass        else   massx       = mass       end
        if  radius    == nothing   radiusx     = nm.radius      else   radiusx     = radius     end
        if  spinI     == nothing   spinIx      = nm.spinI       else   spinIx      = spinI      end
        if  mu        == nothing   mux         = nm.mu          else   mux         = mu         end
        if  Q         == nothing   Qx          = nm.Q           else   Qx          = Q          end

    	Model(Zx, modelx, massx, radiusx, spinIx, mux, Qx)
    end

    # `Base.show(io::IO, m::Model)`  ... prepares a proper printout of the variable  m::Model.
    function Base.show(io::IO, m::Model)
        if      m.model == "Fermi"
            print(io, "Fermi nuclear model for Z = $(m.Z) with mass = $(m.mass), radius R = $(m.radius) fm and ")
        elseif  m.model == "point"
            print(io, "Point nuclear model for Z = $(m.Z) with mass = $(m.mass), radius R = $(m.radius) fm and ")
        elseif  m.model == "uniform"
            print(io, "Uniform nuclear model for Z = $(m.Z) with mass = $(m.mass), radius R = $(m.radius) fm and ")
        else
            error("stop a")
        end

        print(io, "nuclear spin I = $(m.spinI), dipole moment mu = $(m.mu) and quadrupole moment Q = $(m.Q).")
    end

    """
   `Nuclear.fermiA`  ... provides a value::Float64 for the fermi_a parameter.
    """
    fermiA   = 2.3/(4 * log(3))

    """
    `Nuclear.Rrms(A)`
        ... provides a value::Float64 for the root-mean-squared radius R of a uniformly-distributed charge
            density of a nuclues with mass A.
    """
    function Rrms(A)
        return( 0.836 * A^(1/3) + 0.57 )
    end

    """
    `Nuclear.fermiDensity(r, b)`
        ... provides a value::Float64 for a fermi-distributed charge density at r and for given fermi_a
            and fermi_b parameters.
    """
    function fermiDensity(r, b)
        return( 1/(1 + exp((r-b)/fermiA)) )
    end

    """
    `Nuclear.fermiRrms(b::Float64)`
        ... provides a value::Float64 for the root-mean-squared radius of a fermi-distributed charge
            density for given fermi_a and fermi_b parameters.
    """
    function fermiRrms(b::Float64)

        return( sqrt(12 * fermiA^2 * Math.polylogExp(b/fermiA, 5) / Math.polylogExp(b/fermiA, 3)) )
    end

    """
    `Nuclear.computeFermiBParameter(R::Float64)`
        ... computes a value::Float64 for the fermi_b parameter for a fermi-distributed nuclear with
            root-mean square radius R.
    """
    function computeFermiBParameter(R::Float64)

        function eq(b :: Float64)
            return fermiRrms(b) - R
        end

        b = R
        diff = eq(b)

        eps = 1E-14
        count = 0
        max = 100

        while  abs(diff) > eps   &&   count < max
            b = b -  diff

            diff = eq(b)
            count = count + 1
        end
        # println("Found b = $b, Δ = $(eq(b)) after $count iterations")

        return( b )
    end

    """
    `Nuclear.fermiDistributedNucleus(Rrms::Float64, Z::Float64, grid::Radial.Grid)`
        ... computes the effective, radial-dependent charge Z(r) for a Fermi-distributed nucleus with rms
            radius R and nuclear charge Z. The full nuclear potential is then given by V_nuc = - Z(r)/r;
            a potential::Radial.Potential is returned.
    """
    function fermiDistributedNucleus(Rrms::Float64, Z::Float64, grid::Radial.Grid)

        zz = zeros(Float64, grid.NoPoints);   zznew = zeros(Float64, grid.NoPoints);   dx = zeros(Float64, grid.NoPoints)
        ##x function  rho(r::Float64)  1.0 / (1.0 + exp( (r-fermiC_au)/fermiA_au ) )   end
        function  r_rho(r::Float64)  r / (1.0 + exp( (r-fermiC_au)/fermiA_au ) )   end
        function  rr_rho(r::Float64)  r^2 / (1.0 + exp( (r-fermiC_au)/fermiA_au ) )  end

        b = computeFermiBParameter(Rrms);    b < 0  &&  error("Inappropriate R_rms radius.")

        fermiA_au = Defaults.convertUnits("length: from fm to atomic", fermiA)
        fermiC_au = Defaults.convertUnits("length: from fm to atomic", b)
        N = 1. / QuadGK.quadgk(rr_rho, 0., 1.3, rtol=1.0e-10)[1]
        #
        #
        for  i = 2:length(zznew)
            zznew[i] = 1 / grid.r[i] * quadgk(rr_rho, 0.0, 1.0e-3, grid.r[i], rtol=1.0e-8)[1]
            zznew[i] = zznew[i] + quadgk(r_rho, grid.r[i], 1.0, 1.0e+1, 1.0e+3, rtol=1.0e-8)[1]
            zznew[i] = Z * N * zznew[i] * grid.r[i]
        end

        potential = Radial.Potential("nuclear-potential: Fermi-distributed", zznew, deepcopy(grid))
        return( potential )
    end

    """
    `Nuclear.pointNucleus(Z::Float64, grid::Radial.Grid)`
        ... computes the effective, radial-dependent charge Z(r) for a point-like nucleus with nuclear charge Z.
            The full nuclear potential is then given by V_nuc = - Z(r)/r; a potential::Radial.Potential is returned.
    """
    function pointNucleus(Z::Float64, grid::Radial.Grid)

        zz = Z * ones(Float64, grid.NoPoints)

        potential = Radial.Potential("nuclear-potential: point-like", zz, deepcopy(grid))
        return( potential )
    end

    """
    `Nuclear.uniformNucleus(R::Float64, Z::Float64, grid::Radial.Grid)`
        ... computes the effective, radial-dependent charge Z(r) for a uniformly-distributed nucleus with radius
            R [in fm] and nuclear charge Z. The full nuclear potential is then given by V_nuc = - Z(r)/r outside
            the nucleus; a potential::Radial.Potential is returned.
    """
    function uniformNucleus(R::Float64, Z::Float64, grid::Radial.Grid)

        zz = zeros(Float64, grid.NoPoints);   R_au = Defaults.convertUnits("length: from fm to atomic", R)
        ##x println("uniformNucleus()::  R_au = $R_au")

        for i = 1:grid.NoPoints
            if     grid.r[i] <= R_au    zz[i] = Z / (2 * R_au) * (grid.r[i]^2/R_au^2 - 3.) * grid.r[i] * grid.r[i]
            else   zz[i] = Z
            end
        end

        potential = Radial.Potential("nuclear-potential: uniformly-distributed", zz, deepcopy(grid))
        return( potential )
    end

    """
    `Nuclear.nuclearPotential(model::Nuclear.Model, grid::Radial.Grid)`
        ... computes the effective, radial-dependent charge Z(r) for a nucleus with radius R and nuclear charge Z.
            The full nuclear potential is then given by V_nuc = - Z(r)/r; a potential::Radial.Potential
            is returned.
    """
    function nuclearPotential(nm::Nuclear.Model, grid::Radial.Grid)
        zz = zeros(Float64, grid.NoPoints)

        if      nm.model == "point"      potential = Nuclear.pointNucleus(nm.Z, grid)
        elseif  nm.model == "uniform"    potential = Nuclear.uniformNucleus(nm.radius, nm.Z, grid)
        elseif  nm.model == "Fermi"      potential = Nuclear.fermiDistributedNucleus(nm.radius, nm.Z, grid)
        else    error("stop a")
        end

        return( potential )
     end

    """
    `Nuclear.nuclearPotentialDH(model::Nuclear.Model, grid::Radial.Grid, lambda::Float64)`
        ... computes the effective, radial-dependent charge Z(r) for a nucleus with radius R, nuclear charge Z
            and for a Debye-Hueckel screening exp(-lambda r). The full Debye-Hueckel nuclear potential is then
            given by V_nuc = - Z(r)/r; a potential::Radial.Potential is returned.
    """
    function nuclearPotentialDH(nm::Nuclear.Model, grid::Radial.Grid, lambda::Float64)
        pot = nuclearPotential(nm, grid)
        Zr  = pot.Zr
        for  i = 1:length(Zr)   Zr[i] = Zr[i] *exp(-lambda * pot.grid.r[i])    end

        potential = Radial.Potential(pot.name* "+ Debey-Hueckel screening", Zr, deepcopy(pot.grid))
        return( potential )
     end

end # module


text:
export
    capability, devcount, devname, totalmem

function capability(dev::Int)
    major, minor = Cint[0], Cint[0]
    cuDeviceComputeCapability(major, minor, dev)
    Int(major[1]), Int(minor[1])
end

function devcount()
    ref = Cint[0]
    cuDeviceGetCount(ref)
    Int(ref[1])
end

function devname(dev::Int)
    ref = zeros(UInt8, 256)
    cuDeviceGetName(ref, length(ref), dev)
    unsafe_string(pointer(ref))
end

function totalmem(dev::Int)
    ref = Csize_t[0]
    cuDeviceTotalMem(ref, dev)
    Int(ref[1])
end

function attribute(attrib::Int, dev::Int)
    ref = Cint[0]
    cuDeviceGetAttribute(ref, attrib, dev)
    Int(ref[1])
end

function attributes(dev::Int)
    Dict(
    "MAX_THREADS_PER_BLOCK" => attribute(1, dev),
    "MAX_BLOCK_DIM_X" => attribute(2, dev),
    "MAX_BLOCK_DIM_Y" => attribute(3, dev),
    "MAX_BLOCK_DIM_Z" => attribute(4, dev),
    "MAX_GRID_DIM_X" => attribute(5, dev),
    "MAX_GRID_DIM_Y" => attribute(6, dev),
    "MAX_GRID_DIM_Z" => attribute(7, dev),
    "MAX_SHARED_MEMORY_PER_BLOCK" => attribute(8, dev),
    "TOTAL_CONSTANT_MEMORY" => attribute(9, dev),
    "WARP_SIZE" => attribute(10, dev),
    "MAX_PITCH" => attribute(11, dev),
    "MAX_REGISTERS_PER_BLOCK" => attribute(12, dev)
    )
end

function infodevices()
    count = devcount()
    if count == 0
        println("No CUDA-capable device found.")
        return
    end
    for dev = 0:count-1
        name = devname(dev)
        cap = capability(dev)
        mem = round(Int, totalmem(dev) / (1024^2))
        info("device[$dev]: $(name), capability $(cap[1]).$(cap[2]), totalmem = $(mem) MB")
    end
end


text:
abstract type HydroCascade <: PSY.HydroGen end

get_value = PSY.get_value


text:
function get_median(survey, type = "level", outfile = "temp.xlsx")
  url = url_median(survey, type)

  urldownload(url, parser = identity, save_raw = outfile)

  sn = XLSX.readxlsx(outfile)
  sheet = XLSX.sheetnames(sn)[1]

  df = DataFrame(XLSX.readtable(outfile, sheet)...)
  df = ifelse.(df .== "#N/A", missing, df)

  outfile == "temp.xlsx" ? rm(outfile) : nothing

  return df
end


text:
module RoboticsToolkit

# sub-modules
export SpatialMaths
export SerialManipulators
export MobilePlatforms
export Navigation
export Localization

# includes
include("SpatialMaths.jl")
include("SerialManipulators.jl")
include("MobilePlatforms.jl")
include("Navigation.jl")
include("Localization.jl")

end # module RoboticsToolkit


text:
# ---
# title: 892. Surface Area of 3D Shapes
# id: problem892
# author: Indigo
# date: 2021-02-01
# difficulty: Easy
# categories: Math, Geometry
# link: <https://leetcode.com/problems/surface-area-of-3d-shapes/description/>
# hidden: true
# ---
#
# On a `N * N` grid, we place some `1 * 1 * 1 `cubes.
#
# Each value `v = grid[i][j]` represents a tower of `v` cubes placed on top of
# grid cell `(i, j)`.
#
# Return the total surface area of the resulting shapes.
#
#
#
# **Example 1:**
#
#
#
#     Input: [[2]]
#     Output: 10
#
#
# **Example 2:**
#
#
#
#     Input: [[1,2],[3,4]]
#     Output: 34
#
#
# **Example 3:**
#
#
#
#     Input: [[1,0],[0,2]]
#     Output: 16
#
#
# **Example 4:**
#
#
#
#     Input: [[1,1,1],[1,0,1],[1,1,1]]
#     Output: 32
#
#
# **Example 5:**
#
#
#
#     Input: [[2,2,2],[2,1,2],[2,2,2]]
#     Output: 46
#
#
#
#
# **Note:**
#
#   * `1 <= N <= 50`
#   * `0 <= grid[i][j] <= 50`
#
#
## @lc code=start
using LeetCode

function surface_area(grid::Vector{Vector{Int}})
    res = 2 * sum(sum(1 for i in row if i != 0) for row in grid)
    res += sum(row[1] + row[end] + sum(abs, diff(row)) for row in grid)
    return res += sum(grid[1][i] +
                      grid[end][i] +
                      ((length(grid) > 1) ?
                       (sum(abs(grid[row][i] - grid[row - 1][i]) for row in 2:length(grid))) :
                       0) for i in 1:length(grid[1]))
end
function surface_area(grid::Matrix{Int})
    res = 2 * sum(sum(1 for i in row if i != 0) for row in grid)
    res += sum(row[1] + row[end] + sum(abs, diff(row)) for row in eachrow(grid))
    res += sum(col[1] + col[end] + sum(abs, diff(col)) for col in eachcol(grid))
    res
end
## @lc code=end


text:
# This file was generated, do not modify it. # hide
using HTTP, CSV, MLJ, StatsBase, PyPlot
MLJ.color_off() # hide
req = HTTP.get("https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data")
data = CSV.read(req.body,
                header=["Class", "Alcool", "Malic acid",
                        "Ash", "Alcalinity of ash", "Magnesium",
                        "Total phenols", "Flavanoids",
                        "Nonflavanoid phenols", "Proanthcyanins",
                        "Color intensity", "Hue",
                        "OD280/OD315 of diluted wines", "Proline"])
# the target is the Class column, everything else is a feature
y, X = unpack(data, ==(:Class), colname->true);

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
@doc Markdown.doc"""
  lift(a::padic) -> fmpz

  Returns the positive canonical representative in Z. a needs
  to be integral.
"""
function lift(a::padic)
  b = fmpz()
  R = parent(a)
  ccall((:padic_get_fmpz, :libflint), Nothing, (Ref{fmpz}, Ref{padic}, Ref{FlintPadicField}), b, a, R)
  return b
end

function Base.setprecision(f::Generic.Poly{padic}, N::Int)
  f = deepcopy(f)
  for i=1:length(f)
    f.coeffs[i].N = N
  end
  return f
end

function setprecision!(f::Generic.Poly{padic}, N::Int)
  for i=1:length(f)
    f.coeffs[i].N = N
  end
  return f
end


text:
@testset "Datasets Kernel Tests" begin
    df1 = createdummydatasetone()
    dataset1 = Persa.Dataset(df1)

    df2 = createdummydatasettwo()
    dataset2 = Persa.Dataset(df2, 10, 10)

    @testset "Dummy Tests" begin
        @test dataset1.users == 7
        @test dataset1.items == 6
        @test Persa.users(dataset1) == 7
        @test Persa.items(dataset1) == 6

        @test dataset2.users == 10
        @test dataset2.items == 10
        @test Persa.users(dataset2) == 10
        @test Persa.items(dataset2) == 10

        @test Persa.size(dataset1) == (7, 6)
        @test Persa.size(dataset2) == (10, 10)

        @test Persa.length(dataset1) == 35
        @test Persa.length(dataset2) == 35
    end

    @testset "Index Tests" begin
        @testset "Cartesian Index Tests" begin
            for i = 1:size(df1)[1]
                user = df1[:user][i]
                item = df1[:item][i]
                rating = df1[:rating][i]
                @test Persa.value(dataset1[user, item]) == rating
            end

            for i = 1:size(df2)[1]
                user = df2[:user][i]
                item = df2[:item][i]
                rating = df2[:rating][i]
                @test Persa.value(dataset2[user, item]) == rating
            end

            for i = 1:length(dataset1)
                (user, item, rating) = dataset1[i]
                @test dataset1[user, item] == rating
            end

            for i = 1:length(dataset2)
                (user, item, rating) = dataset2[i]
                @test dataset2[user, item] == rating
            end
        end

        @testset "Colun Index Tests" begin
            for user in 1:Persa.users(dataset1)
                for (item, rating) in dataset1[user, :]
                    @test dataset1[user, item] == rating
                end
            end

            for item in 1:Persa.items(dataset1)
                for (user, rating) in dataset1[:, item]
                    @test dataset1[user, item] == rating
                end
            end

            for user in 1:Persa.users(dataset2)
                for (item, rating) in dataset2[user, :]
                    @test dataset2[user, item] == rating
                end
            end

            for item in 1:Persa.items(dataset2)
                for (user, rating) in dataset2[:, item]
                    @test dataset2[user, item] == rating
                end
            end
        end
    end

    @testset "Matrix Conversion Tests" begin
        matrix = Array(dataset1)

        for user in 1:Persa.users(dataset1)
            for item in 1:Persa.items(dataset1)
                @test Persa.value(dataset1[user, item]) === matrix[user, item]
            end
        end

        matrix = Array(dataset2)

        for user in 1:Persa.users(dataset2)
            for item in 1:Persa.items(dataset2)
                @test Persa.value(dataset2[user, item]) === matrix[user, item]
            end
        end
    end

    @testset "Interator Tests" begin
        for (user, item, rating) in dataset1
            @test dataset1[user, item] == rating
        end

        for (user, item, rating) in dataset2
            @test dataset2[user, item] == rating
        end
    end
end


text:
"""
```
gensys2(m::AbstractDSGEModel, Γ0::Matrix{Float64}, Γ1::Matrix{Float64}, C::Vector{Float64},
        Ψ::Matrix{Float64}, Π::Matrix{Float64}, TTT::Matrix{Float64}, RRR::Matrix{Float64},
        CCC::Vector{Float64}, T_switch::Int)
```

calculates the state space transition matrices when a temporary alternative policy applies.
This function is **not** the same as the gensys2 written by Chris Sims, which computes
a second-order perturbation.
"""
function gensys2(m::AbstractDSGEModel, Γ0::Matrix{Float64}, Γ1::Matrix{Float64}, C::Vector{Float64},
                 Ψ::Matrix{Float64}, Π::Matrix{Float64}, TTT::Matrix{Float64}, RRR::Matrix{Float64},
                 CCC::Vector{Float64}, T_switch::Int)

    Γ0_til, Γ1_til, Γ2_til, C_til, Ψ_til =
        gensys_to_predictable_form(Γ0, Γ1, C, Ψ, Π; use_sparse = (haskey(get_settings(m), :gensys2_sparse_matrices) &&
                                                                  get_setting(m, :gensys2_sparse_matrices)))

    Tcal = Vector{Matrix{Float64}}(undef, T_switch)
    Rcal = Vector{Matrix{Float64}}(undef, T_switch)
    Ccal = Vector{Vector{Float64}}(undef, T_switch)

    Tcal[end] = TTT
    Rcal[end] = RRR
    Ccal[end] = CCC

    for t = 1:(T_switch-1)
        tmp = Γ2_til * TTT + Γ0_til
        Tcal[end-t] = tmp \ Γ1_til
        Rcal[end-t] = tmp \ Ψ_til
        Ccal[end-t] = tmp \ (C_til - Γ2_til * CCC)

        TTT = Tcal[end-t]
        CCC = Ccal[end-t]
        RRR = Rcal[end-t]
    end
    return Tcal, Rcal, Ccal
end

function gensys2(m::AbstractDSGEModel, Γ0s::Vector{Matrix{Float64}}, Γ1s::Vector{Matrix{Float64}},
                 Cs::Vector{Vector{Float64}}, Ψs::Vector{Matrix{Float64}}, Πs::Vector{Matrix{Float64}},
                 TTT::Matrix{Float64}, RRR::Matrix{Float64},
                 CCC::Vector{Float64}, T_switch::Int)

    ntil = length(Γ0s)
    use_sparse = haskey(get_settings(m), :gensys2_sparse_matrices) &&
        get_setting(m, :gensys2_sparse_matrices)
    Γ0_tils = use_sparse ? Vector{SparseMatrixCSC{Float64, Int64}}(undef, ntil) : Vector{Matrix{Float64}}(undef, ntil)
    Γ1_tils = Vector{Matrix{Float64}}(undef, ntil)
    Γ2_tils = use_sparse ? Vector{SparseMatrixCSC{Float64, Int64}}(undef, ntil) : Vector{Matrix{Float64}}(undef, ntil)
    C_tils = Vector{Vector{Float64}}(undef, ntil)
    Ψ_tils = Vector{Matrix{Float64}}(undef, ntil)

    for i in 1:length(Γ0s)
        Γ0_tils[i], Γ1_tils[i], Γ2_tils[i], C_tils[i], Ψ_tils[i] =
            gensys_to_predictable_form(Γ0s[i], Γ1s[i], Cs[i], Ψs[i], Πs[i];
                                       use_sparse = use_sparse)
    end

    Tcal = Vector{Matrix{Float64}}(undef, T_switch)
    Rcal = Vector{Matrix{Float64}}(undef, T_switch)
    Ccal = Vector{Vector{Float64}}(undef, T_switch)

    # maybe allocate sparse version of TTT b/c going to do a lot of calculations with it
    # and also make sure to check for essentially zero terms.
    # leave Gamma1_til as a full array. Only make Gamma2, Gamma0 sparse.
    Tcal[end] = TTT
    Rcal[end] = RRR
    Ccal[end] = CCC

    for t = 1:(T_switch-1)
        tmp = Γ2_tils[end-t] * TTT + Γ0_tils[end-t]
        Tcal[end-t] = tmp \ Γ1_tils[end-t]
        Rcal[end-t] = tmp \ Ψ_tils[end-t]
        Ccal[end-t] = tmp \ (C_tils[end-t] - Γ2_tils[end-t] * CCC)

        TTT = Tcal[end-t]
        CCC = Ccal[end-t]
        RRR = Rcal[end-t]
    end

    return Tcal, Rcal, Ccal
end

function gensys_to_predictable_form(Γ0::AbstractMatrix{S}, Γ1::AbstractMatrix{S}, C::AbstractVector{S},
                                    Ψ::AbstractMatrix{S}, Π::AbstractMatrix{S};
                                    use_sparse::Bool = false) where {S <: Real}

    nonexp_eq_ind = map(i -> all(Π[i, :] .== 0.), 1:size(Π, 1))
    if use_sparse
        Γ0_til = spzeros(S, size(Γ0)...)
        Γ1_til = zeros(S, size(Γ1))   # this one has to be a dense matrix for left-divide
        Γ2_til = spzeros(S, size(Γ0)...)
    else
        Γ0_til = zeros(S, size(Γ0))
        Γ1_til = zeros(S, size(Γ1))
        Γ2_til = zeros(S, size(Γ0))
    end
    C_til = C
    Ψ_til = Ψ

    # In this loop, note that rudimentary tests in REPL indicate that, for matrices of the size
    # which we typically encounter in DSGEs, using a view for rows is still faster
    for row in 1:length(nonexp_eq_ind)
        if nonexp_eq_ind[row]
            # Not expectational equation
            Γ1_til[row, :]  = view(Γ1, row, :)
            Γ0_til[row, :]  = view(Γ0, row, :)
        else
            # expectational equation, assumed to take the form zₜ = Eₜ₋₁[zₜ] + ηₜ,
            # hence why we may use findfirst (no need to do findall or some other inference approch)
            Γ0_til[row, findfirst(x -> x > 0., view(Γ1, row, :))] = -1.
            Γ2_til[row, findfirst(x -> x > 0., view(Γ0, row, :))] = 1.
        end
    end

    return Γ0_til, Γ1_til, Γ2_til, C_til, Ψ_til
end


text:
insts = OrderedDict()
insts["minimal"] = [
    ((2, 4, 1.0),),
    ((2, 4, 0.5),),
    ]
insts["fast"] = [
    ((15, 20, 1.0),),
    ((15, 20, 0.25),),
    ((100, 100, 1.0),),
    ((100, 100, 0.15),),
    ((500, 100, 1.0),),
    ((500, 100, 0.15),),
    ]
insts["various"] = [
    ((250, 1000, 0.1),),
    ((500, 5000, 0.05),),
    ((500, 2000, 1.0),),
    ]
return (LinearOptNative, insts)


text:
server = LanguageServerInstance(IOBuffer(), IOBuffer(), dirname(Pkg.Types.Context().env.project_file), first(Base.DEPOT_PATH))
server.runlinter = true

LanguageServer.process(LanguageServer.parse(LanguageServer.JSONRPC.Request, Dict("jsonrpc"=>"2.0","id"=>0,"method"=>"initialize","params"=>init_request)), server)

testtext = """
decode_overlong
"""
LanguageServer.process(LanguageServer.JSONRPC.Request{Val{Symbol("textDocument/didOpen")},LanguageServer.DidOpenTextDocumentParams}(0, LanguageServer.DidOpenTextDocumentParams(LanguageServer.TextDocumentItem("testdoc", "julia", 0, testtext))), server)

doc = LanguageServer.getdocument(server, LanguageServer.URI2("testdoc"))

@test !isempty(LanguageServer.process(LanguageServer.JSONRPC.Request{Val{Symbol("textDocument/codeAction")},LanguageServer.CodeActionParams}(0, LanguageServer.CodeActionParams(LanguageServer.TextDocumentIdentifier("testdoc"), LanguageServer.Range(LanguageServer.Position(0,1),LanguageServer.Position(0,1)), LanguageServer.CodeActionContext(LanguageServer.Diagnostic[doc.diagnostics[1]], missing))), server))


text:
# This file is a part of RadiationDetectorDSP.jl, licensed under the MIT License (MIT).

using RadiationDetectorDSP
using Test

using Random, LinearAlgebra, Statistics
using DSP

@testset "Filters" begin
    # Runs Biquad filter in Direct Form 1, for comparison (DSP.jl uses
    # Transposed Direct Form 2):
    df1_filt(f::Biquad, input::Vector{<:Real}) =
        filt([f.b0, f.b1, f.b2], [one(f.a1), f.a1, f.a2], input)

    plot(args...) = nothing
    plot!(args...) = nothing
    hline!(args...) = nothing

    @testset "test_csa_response" begin
        t_drift, τ_rise, τ_decay = 20, 20, 500
        current_signal = vcat(fill(0.0, 10), fill(1.0 / t_drift, t_drift), fill(0.0, 200))
        plot(cumsum(current_signal))
        output = filt(simple_csa_response_filter(τ_rise, τ_decay, 1), current_signal)
        plot!(output)
        output_deconv = filt(inv_cr_filter(τ_decay), output)
        plot!(output_deconv)
        tail = output_deconv[150:end]
        # Tail of reco should be flat:
        @test var(tail) < 1e-5
    end

    @testset "rc_filter" begin
        RC = 20
        X = vcat(fill(0.0, 10), fill(1.0, 40))
        plot(X)
        output = filt(rc_filter(RC), X)
        # output = df1_filt(rc_filter(RC), X)
        plot!(output)
        hline!([1 - exp(-1)])
        #output_deconv = filt(inv_rc_filter(RC), output)
        #plot!(output_deconv)
    end

    @testset "rc_filter" begin
        RC = 10
        X = vcat(fill(0.0, 10), fill(1.0, 30))
        plot(X)
        output = filt(cr_filter(RC), X)
        # output = df1_filt(cr_filter(RC), X)
        plot!(output)
        output_deconv = filt(inv_cr_filter(RC), output)
        # output_deconv = df1_filt(inv_cr_filter(RC), output)
        plot!(output_deconv)
        hline!([exp(-1)])
        @test X ≈ output_deconv
    end

    @testset "crmod_filter" begin
        RC = 10
        X = vcat(fill(0.0, 10), fill(1.0, 30))
        plot(X)
        output = filt(crmod_filter(RC), X)
        # output = df1_filt(crmod_filter(RC), X)
        plot!(output)
        output_deconv = filt(inv_crmod_filter(RC), output)
        # output_deconv = df1_filt(inv_crmod_filter(RC), output)
        plot!(output_deconv)
        hline!([exp(-1)])
        @test X ≈ output_deconv
    end
end


text:
using Base.Test
using StateEstimation

println("1. Testing State constructors and methods...")
include("state_tests.jl")

println("2. Testing System constructors and methods...")
include("system_tests.jl")

println("3. Testing Observer constructors and methods...")
include("observer_tests.jl")

println("4. Testing UnscentedTransform constructors and methods...")
include("unscented_transform_tests.jl")

println("5. Testing KalmanFilter constructors and methods...")
include("kalman_filter_tests.jl")

println("6. Testing ExtendedKalmanFilter constructors and methods...")
include("extended_kalman_filter_tests.jl")

println("7. Testing MultiTargetFilter constructors and methods...")
include("multi_target_filter_tests.jl")

println("8. Testing LeastSquaresEstimator constructors and methods...")
include("least_squares_tests.jl")

println("9. Testing EstimatorHistory constructors and methods...")
include("archive_tests.jl")


text:
using ..CUDASupport
import ..TopOpt: @init_cuda, whichdevice, CPU, GPU
@init_cuda()

export  CPU,
        GPU,
        whichdevice,
        @define_cu,
        @thread_local_index,
        @total_threads_per_block,
        @block_index,
        @total_blocks,
        @thread_global_index,
        @total_threads,
        callkernel,
        @mapreduce_block

whichdevice(s::AbstractArray) = s isa GPUArrays.GPUArray ? GPU() : CPU()

function getfieldnames end
function cufieldnames end

macro define_cu(T, fields...)
    _define_cu(T, fields...)
end

function _define_cu(T, fields...)
    all_fields = gensym()
    args = gensym()
    if eltype(fields) <: QuoteNode
        field_syms = Tuple(field.value for field in fields)
    elseif eltype(fields) <: Symbol || length(fields) == 0
        field_syms = fields
    else
        throw("Unsupported fields.")
    end
    esc(quote
        $all_fields = Tuple(fieldnames($T))
        @eval @inline GPUUtils.getfieldnames(::Type{<:$T}) = $(Expr(:$, all_fields))
        @inline GPUUtils.cufieldnames(::Type{<:$T}) = $field_syms
        $args = Expr[]
        for fn in $all_fields
            push!($args, :(GPUUtils._cu(s, s.$fn, $(Val(fn)))))
        end
        @eval begin
            function CuArrays.cu(s::$T)
                $T($(Expr(:$, Expr(:..., args))))
            end
        end
    end)
end

function _cu(s::T, f::F, ::Val{fn}) where {T, F, fn}
    if fn ∈ GPUUtils.cufieldnames(T)
        if F <: AbstractArray
            CuArrays.CuArray(f)
        else
            CuArrays.cu(f)
        end
    else
        f
    end
end

macro thread_local_index()
    :((threadIdx().z - 1) * blockDim().y * blockDim().x + (threadIdx().y - 1) * blockDim().x + threadIdx().x)
end

macro total_threads_per_block()
    :(blockDim().z * blockDim().y * blockDim().x)
end

macro block_index()
    :(blockIdx().x + (blockIdx().y - 1) * gridDim().x + (blockIdx().z - 1) * gridDim().x * gridDim().y)
end

macro total_blocks()
    :(gridDim().z * gridDim().x * gridDim().y)
end

macro thread_global_index()
    :((@block_index() - 1) * (blockDim().x * blockDim().y * blockDim().z) + @thread_local_index())
end

macro total_threads()
    :(@total_blocks() * @total_threads_per_block())
end

macro mapreduce_block(indvar, limit, op, T, LMEM, result, mapexpr)
    _mapreduce_block(indvar, limit, op, T, LMEM, result, mapexpr)
end
function _mapreduce_block(indvar, limit, op, T, LMEM, result, mapexpr)
    offset = gensym()
    out = gensym()
    tmp_local = gensym()
    local_index = gensym()
    esc(quote
        $indvar = @thread_global_index()
        $offset = @total_threads()
        $out = zero($T)
        # # Loop sequentially over chunks of input vector
        while $indvar <= $limit
            $out = $op($out, $mapexpr)
            $indvar += $offset
        end

        # Perform parallel reduction
        $tmp_local = @cuStaticSharedMem($T, $LMEM)
        $local_index = @thread_local_index()
        $tmp_local[$local_index] = $out
        sync_threads()

        $offset = @total_threads_per_block() ÷ 2
        while $offset > 0
            if ($local_index <= $offset)
                $tmp_local[$local_index] = $op($tmp_local[$local_index], $tmp_local[$local_index + $offset])
            end
            sync_threads()
            $offset = $offset ÷ 2
        end
        if $local_index == 1
            $result[@block_index()] = $tmp_local[1]
        end
    end)
end

function callkernel(dev, kernel, args)
    blocks, threads = getvalidconfig(dev, kernel, args)
    #@show blocks, threads
    @cuda blocks=blocks threads=threads kernel(args...)

    return
end

# From CuArrays.jl
function getvalidconfig(dev, kernel, parallel_args)
    R = parallel_args[1]
    Rlength = length(R)
    Ssize = size(R)
    Slength = prod(Ssize)
    GC.@preserve parallel_args begin
        parallel_kargs = cudaconvert.(parallel_args)
        parallel_tt = Tuple{Core.Typeof.(parallel_kargs)...}
        parallel_kernel = cufunction(kernel, parallel_tt)

        # we are limited in how many threads we can launch...
        ## by the kernel
        kernel_threads = CUDAnative.maxthreads(parallel_kernel)
        ## by the device
        block_threads = (x=CUDAdrv.attribute(dev, CUDAdrv.MAX_BLOCK_DIM_X),
                         y=CUDAdrv.attribute(dev, CUDAdrv.MAX_BLOCK_DIM_Y),
                         total=CUDAdrv.attribute(dev, CUDAdrv.MAX_THREADS_PER_BLOCK))

        # figure out a legal launch configuration
        y_thr = min(nextpow(2, Rlength ÷ 512 + 1), 512, block_threads.y, kernel_threads)
        x_thr = min(512 ÷ y_thr, Slength, block_threads.x,
                    ceil(Int, block_threads.total/y_thr),
                    ceil(Int, kernel_threads/y_thr))

        blk, thr = (Rlength - 1) ÷ y_thr + 1, (x_thr, y_thr, 1)
        blk = min(blk, ceil(Int, Rlength / prod(thr)))
    end
    return blk, thr
end


text:
"""
    AbstractSimSettings

Abstract supertype for [`SimSettings`](@ref) object and variants.
"""
abstract type AbstractSimSettings end

"""
    SimSettings <: AbstractSimSettings

Holds settings for the simulation, in a `Ruleset` or `SimData` object.
"""
Base.@kwdef struct SimSettings{B,P,Op,C,T} <: AbstractSimSettings
    boundary::B = Remove()
    proc::P = SingleCPU()
    opt::Op = NoOpt()
    cellsize::C = 1
    timestep::T = nothing
end

boundary(s::AbstractSimSettings) = s.boundary
proc(s::AbstractSimSettings) = s.proc
opt(s::AbstractSimSettings) = s.opt
cellsize(s::AbstractSimSettings) = s.cellsize
timestep(s::AbstractSimSettings) = s.timestep


text:
"""
    particle_in_image(x::Real, y::Real, z::Real, hsml::Real,
                                halfXsize::Real, halfYsize::Real, halfZsize::Real)

Checks if a periodically mapped particle is still in the image frame.
"""
@inline function particle_in_image( x::Real, y::Real, z::Real, hsml::Real,
                                    halfsize::Array{T}) where T

    if (( x - hsml ) > halfsize[1] || ( x + hsml ) < -halfsize[1] ||
        ( y - hsml ) > halfsize[2] || ( y + hsml ) < -halfsize[2] ||
        ( z - hsml ) > halfsize[3] || ( z + hsml ) < -halfsize[3] )

        return false
    else
        return true
    end
end

"""
    particle_in_image( pos::Array{T}, hsml::T,
                       halfsize::Array{Float64}) where T

Checks if a periodically mapped particle is still in the image frame.
"""
@inline function particle_in_image( pos::Vector{T}, hsml::T,
                                    halfsize::Vector{Float64}) where T

    if ( ( pos[1] - hsml ) > halfsize[1] || ( pos[1] + hsml ) < -halfsize[1]  ||
         ( pos[2] - hsml ) > halfsize[2] || ( pos[2] + hsml ) < -halfsize[2]  ||
         ( pos[3] - hsml ) > halfsize[3] || ( pos[3] + hsml ) < -halfsize[3]  )

        return false
    else
        return true
    end
end

"""
    add_subtr_box(pos::T, boxsize::T) where T

Helper function to add or subtract the boxsize from the given position
"""
@inline function add_subtr_box(pos::T, boxsize::T) where T
    if pos > 0.0
        return pos - boxsize
    else
        return pos + boxsize
    end
end

"""
    find_position_periodic( pos, k, boxsize )

Performs a periodic mapping of the particle position.
"""
@inline function find_position_periodic( pos, k, boxsize )

    # re-define here to avoid type-instability
    box   = eltype(pos[1])(boxsize)

    x = (k & 0x1) == 0 ? pos[1] : add_subtr_box(pos[1], box)

	y = (k & 0x2) == 0 ? pos[2] : add_subtr_box(pos[2], box)

	z = (k & 0x4) == 0 ? pos[3] : add_subtr_box(pos[3], box)

    return x, y, z
end

@inline function find_position_periodic!( _pos::Array{T}, pos::Array{T}, k::Integer, boxsize::Float64) where T

    # re-define here to avoid type-instability
    box   = eltype(pos[1])(boxsize)

    _pos[1] = (k & 0x1) == 0 ? pos[1] : add_subtr_box(pos[1], box)

	_pos[2] = (k & 0x2) == 0 ? pos[2] : add_subtr_box(pos[2], box)

	_pos[3] = (k & 0x4) == 0 ? pos[3] : add_subtr_box(pos[3], box)

    return _pos
end


text:
export animate_evolution

# This internal function creates a closure for easy plotting for each particle.
# arguments are mostly similar to animate_evolution
function setup_animation(p::AbstractParticle, bd:: Billiard, t::AbstractFloat,
                         ax::PyPlot.PyCall.PyObject, raysplitters = nothing;
                         dt = 0.01, taillength::Int, tailcolor = "C0",
                         particle_kwargs = NamedTuple(),
                         tail_kwargs = NamedTuple())

    # get data to plot
    x,y,vx,vy = timeseries(p, bd, t, raysplitters, dt = dt)

    # initial plot
    if particle_kwargs != nothing
        point, arrow = plot_particle(x[1], y[1], vx[1], vy[1]; ax = ax, zorder = 20,
                                     particle_kwargs...)
    end

    tail, = ax.plot(x[1:2], y[1:2], zorder = 1, color = tailcolor; tail_kwargs...)

    # frame counter
    count = 2

    function plot_frame()
        if particle_kwargs != nothing
            # replot arrow
            arrow.remove()
            # replot point
            point.remove()
            point, arrow = plot_particle(x[count], y[count], vx[count], vy[count];
                                         ax = ax, zorder = 20, particle_kwargs...)
        end
        # set tail data
        @views tail.set_xdata(x[clamp(count-taillength, 1, count):count])
        @views tail.set_ydata(y[clamp(count-taillength, 1, count):count])

        count += 1 # increment frame counter
    end

    function skip_frame()
        count += 1
    end

    return plot_frame, skip_frame, length(x)

end

"""
    animate_evolution(ps, bd, colnumber [, raysplitters]; kwargs...)
Animate the evolution of a vector of particles `ps` in billiard `bd` for a
total time `t` (always considered float). Optionally enable ray-splitting.

### Evolution kwargs
  * `dt = 0.01` : Time resolution used for production of time series (see
    [`timeseries!`](@ref). It is not recommended to significantly increase this value,
    to preserve the smoothness of the orbits.
  * `frameskip = 5` : The amount of `dt`-steps performed each frame.
    Increasing either `frameskip` and `dt` makes the animation progress faster.
  * `tailtime = 1.0` : The length of the "tail" trailing the particle in time
    units.
  * `resetting = reset_billiard!` : function called after evolving each individual
    particle in the billiard (so that ray-splitting doesn't break).
### Colors & plotting kwargs
  * `colors` : An array of valid Matplotlib colors for the "tails". If `colors`
    is shorter than `ps`, colors are reused. Defaults to the standard
    Matplotlib color sequence.
  * `particle_kwargs::NamedTuple` : Additional keyword arguments passed to the
    `plot` function for particles.
  * `tail_kwargs::NamedTuple`: Additional keyword arguments passed to the `plot`
    function for "tails" (line plot).
### Exporting and axis kwargs
  * `figsize = (7.2, 7.2))` : Size for new figure (if one is created).
    Must be divisible by 2 if you want to save the animation.
  * `dpi = 100` : DPI for saving the figures (and thus for also the resulting video).
  * `ax = (figure(figsize = figsize); plot(bd); gca())` : axis to plot on.
  * `savename = nothing` : If given the animation is exported to
    an `savetype` file (requires ffmpeg). The name can include path.
  * `savetype = "mp4"` or `"gif"`.
  * `disable_axis = false` : Remove the axis splines.
  * `deletefigs = true` : To create the animation a lot of figures are saved in
    the save directory and are deleted after the animation is done. You can choose
    to keep them.
  * `framerate = 20` : Animation framerate.
"""
function animate_evolution(ps::AbstractVector{<:AbstractParticle{T}},
                           bd::Billiard, t, raysplitters = nothing;
                           dt  = 0.01, frameskip = 5, tailtime = 1.0,
                           colors = ["C$i" for i ∈ 0:9],
                           particle_kwargs = NamedTuple(),
                           tail_kwargs = NamedTuple(),
                           figsize = (7.2, 7.2),
                           dpi = 100,
                           ax = (PyPlot.figure(figsize = figsize, dpi = dpi); ax = PyPlot.gca(); plot(bd; ax = ax); ax),
                           savename = nothing, deletefigs = true,
                           disable_axis = false, framerate = 20,
                           savetype = "mp4",
                           resetting = reset_billiard!
                           ) where {T}

    disable_axis && ax.axis("off")
    nps = length(ps)
    taillength = round(Int, tailtime/dt)
    savename != nothing && (colnumbers = Int[])
    if savename != nothing
        fig = ax.get_figure()
        width, height = fig.get_size_inches()
        if 100width % 2 != 0 || 100height % 2 != 0
            @warn "100*width or 100*height of figure are not divisible by 2. "*
            "FFMPEG won't be able to produce the animation... "*
            "Use `fig.set_size_inches` to ensure correct size."
        end
    end

    plotframes = Vector{Function}(undef, nps)
    skipframes = Vector{Function}(undef, nps)
    tslengths = Vector{Int}(undef, nps)

    for (i, p) ∈ enumerate(ps)

        resetting(bd)

        pf, sf, tl = setup_animation(
            p, bd, T(t), ax, raysplitters, dt = dt, taillength = taillength,
            tailcolor = colors[(i-1)%length(colors) + 1],
            tail_kwargs=tail_kwargs, particle_kwargs=particle_kwargs
        )

        plotframes[i] = pf
        skipframes[i] = sf
        tslengths[i] = tl
    end

    maxframe = minimum(tslengths)
    j = 1

    for i ∈ 2:maxframe
        if (i-1)%frameskip == 0
            for pf in plotframes; pf(); end
            sleep(0.000001) # this updates PyPlot live
            if savename != nothing
                s = savename*"_$(j).png"
                PyPlot.savefig(s, dpi = dpi)
                push!(colnumbers, j)
                j += 1
            end
        else
            for sf in skipframes; sf(); end
        end

    end

    if savename != nothing
        try
            if savetype == "mp4"
                anim = `ffmpeg -y -framerate $(framerate) -start_number 1 -i $(savename)_%d.png
                -c:v libx264 -pix_fmt yuv420p -preset veryslow -profile:v high -level 5.2 $(savename).mp4`
                run(anim)
            elseif savetype == "gif"
                # Code from https://github.com/JuliaPlots/Plots.jl/commit/7f7b543e1802aa85a4bc9c1069d9b9f0a00b5340
                bas = basename(savename)
                bas == savename && (bas = "")
                pal = joinpath(bas, "palette.png")
                # create a pallete
                run(`ffmpeg -v 0 -start_number 1 -i $(savename)_%d.png -vf palettegen -y $pal`)
                # then apply the palette to get better results
                run(`ffmpeg -v 0 -framerate $framerate -start_number 1 -i $(savename)_%d.png -i $pal -lavfi paletteuse -y $(savename).gif`)
                rm(pal)
                println("Done creating $(savename).gif")
            end
        catch e
            println("When attempting to create video through ffmpeg, got error: ")
            showerror(stdout, e)
        end

        if deletefigs
            for i in colnumbers
                rm(savename*"_$(i).png")
            end
        end
    end
end

animate_evolution(p::AbstractParticle, args...; kwargs...) =
animate_evolution([p], args...; kwargs...)


text:
using WebSockets

function client_one_message(ws)
    printstyled(stdout, "\nws|client input >  ", color=:green)
    msg = readline(stdin)
    if writeguarded(ws, msg)
        msg, stillopen = readguarded(ws)
        println("Received:", String(msg))
        if stillopen
            println("The connection is active, but we leave. WebSockets.jl will close properly.")
        else
            println("Disconnect during reading.")
        end
    else
        println("Disconnect during writing.")
    end
end
function main()
    while true
        println("\nSuggestion: Run 'minimal_server.jl' in another REPL")
        println("\nWebSocket client side. WebSocket URI format:")
        println("ws:// host [ \":\" port ] path [ \"?\" query ]")
        println("Example:\nws://127.0.0.1:8080")
        println("Where do you want to connect? Empty line to exit")
        printstyled(stdout, color = :green,  "\nclient_repl_input >  ")
        wsuri = readline(stdin)
        wsuri == "" && break
        res = WebSockets.open(client_one_message, wsuri)
    end
    println("Have a nice day")
end

main()


text:
module NonlinearOptics
    using Polynomials
    using OrdinaryDiffEq, DiffEqBase, StochasticDiffEq, Tensors
    using RecursiveArrayTools, Parameters

    import Unitful
    import Unitful: ps, m, J, THz, °, THz2π, W
    import Unitful: c0, ħ, k

    using RecipesBase

    import DiffEqCallbacks: PeriodicCallback
    import Base: getindex, setindex!
    import DiffEqBase: solve

    abstract type AbstractNLSEProblem{uType, tType, zType} <: DEProblem end

    abstract type AbstractNLSESolution{T,N} <: AbstractTimeseriesSolution{T,N} end

    abstract type AbstractNLSEAlgorithm <: DEAlgorithm end

    const ħ_over_kBT = ħ/(k*300(Unitful.K)) # T = 300 K, ps

    Unitful.preferunits(ps)

    include("interface/fit_util.jl")
    include("interface/types.jl")
    include("waveguides/waveguides.jl")
    include("interface/models.jl")
    include("interface/solvertypes.jl")
    include("interface/interface.jl")
    include("optics/utils.jl")
    include("optics/dispersion.jl")
    include("optics/pulses.jl")
    include("interface/build.jl")
    include("interface/solve.jl")
    include("materials/materials.jl")
    include("plot_recipes.jl")

    export NLSEProblem, NLSESolution, SymmetrizedSplitStep
    export prob_bright_soliton, prob_linear_pulse, prob_dispersive_pulse,
            prob_LL
    export solve
    export build_problem

    export wavelength, frequency, c0
    export OpticalAttr
    export Mode, Waveguide
    export PulsedLaser, CWLaser
    export frequency, wavelength, getω, get_beta, get_FSR, get_groupindex,
           get_property, get_label
    export add_mode!, add_interaction!
    export add_straight!, add_turn!
    export currentmedium, lastmedium, newmedium
    export Model, ToyModel
    export DynamicNLSE, DynamicLL, DynamicIkeda,
           SteadyStateNLSE, SteadyStateLL, SteadyStateIkeda,
           StochasticNLSE
    export FT
    export Glass, Silicon, Diamond

end # module


text:
using Base.Filesystem

function writehead(tgrfileObj )
    for i=1:512
        Base.write(tgrfileObj,0x00::UInt8 )
    end
end
function writeDataSize(tgrfileObj )
    Base.write(tgrfileObj,0x00::UInt8 )
    Base.write(tgrfileObj,0x02::UInt8 )
end

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

function write(tgrfileObj, cfgdata::cfgStruct)
        writehead(tgrfileObj )
        cfgdatas = Array{UInt16,1}(undef,7)
        cfgdatas.=0x0000
        cfgdatas[1] =floor(UInt16,(cfgdata.meterperSeg *5000+0.001) )
        cfgdatas[2] =cfgdata.chn2SegBegIdx
        cfgdatas[3] =floor(UInt16,(cfgdata.reflectorFactor *10000+0.001) )
        cfgdatas[4] =floor(UInt16,(cfgdata.attenutionFactor *10000+0.001) )
        cfgdatas[5] =cfgdata.scanrate
        # cfgdatas[1] =floor(Uint16,(cfgdata.meterperSeg *5000) )
        cfgdatas[6] =cfgdata.calSamples
        cfgdatas[7] =cfgdata.SegNumber
        cfgString= reinterpret(UInt8, cfgdatas)

        Base.write(tgrfileObj,cfgString   )
        pos2=position( tgrfileObj )
        println( "Position  :",pos2," ,",cfgdata.meterperSeg )
        # for i =1:7
        #     println( cfgdatas[i] )
        # end
        # println()
        # exit()

end

function write(  tgrfileObj,dtdata::DatetimeData )
    # dt
    pos1=position( tgrfileObj )
    Base.write(tgrfileObj,dtdata.dt )
    writeDataSize(tgrfileObj )
    pos2=position( tgrfileObj )
    # println( "write Length :",pos2-pos1)

end

function write(  tgrfileObj,row::Array{UInt16,1} )
    pos1=position( tgrfileObj )
    data8Array= reinterpret( UInt8,  row)
    Base.write(tgrfileObj, data8Array );
    pos2=position( tgrfileObj )
    # println( "write Length :",pos2-pos1)
end

function Trasfer1File( dataFilename ,cfgfilename,tgrfileObj)
    chnid=1
    startMeters=0
    endmeters=1000000
    maxTimeCnt1=-1
    ProtectMeters=5000
    cfgObj = cfgStruct()
    # aRNG = RNG()
    CalObj = spaceCumSumData( )
    read!(cfgObj,cfgfilename)
    write(tgrfileObj, cfgObj)
    show( cfgObj  )
    # getRNG!( aRNG,cfgObj, startMeters ,endmeters ,ProtectMeters,chnid )
    # show( aRNG)
    # CalObj.CfgInfoObj = cfgObj
    # InitObj(CalObj, aRNG.AlrmProcStartID,aRNG.AlrmProcEndtID1 , aRNG.startID1 ,aRNG.endID1,maxTimeCnt1 )
    readinDatainFile( dataFilename,CalObj,tgrfileObj,maxTimeCnt1 )
end

function readinDatainFile( filename,spaceCumSumDataObj ,tgrfileObj,Count=-1)
    global extensionName
    fn = basename( filename )
    extensionName  =uppercase( splitext( fn)[2]  )
    SegNumber=8192
    println(  abspath( filename ) ,"   ",fn,"   " ,extensionName ,"   ",splitext( fn)[1]   )
    if  extensionName == ".RAW1"  ||   extensionName  == ".RAW"
        SegNumber =8192*512
    elseif  extensionName  in (   ".RAW3"  ,".RAW4"  ,".LC")
        SegNumber =8192*1
    end

    rawdatafileObj= open( filename,"r")
    readin(rawdatafileObj,NULLDATA)
    cnt=1
    dtdata = DatetimeData("")
    row = Array{UInt16,1}(undef,SegNumber)

    while( !eof(rawdatafileObj ) )
        try
            readin!(rawdatafileObj,dtdata)
            readin!(rawdatafileObj,row)
            # setDTValue( spaceCumSumDataObj , dtdata.dt)
            write(  tgrfileObj,dtdata )
            #@time
            # saveLogInfo(spaceCumSumDataObj ,row  )
            write(  tgrfileObj,row )
            # postProcFeatureData(spaceCumSumDataObj ,row )
            cnt= cnt+1
            if cnt > Count && Count >= 0
                break
            end
        catch(e)
            println("File End :$e")
        end
        end
    close( rawdatafileObj )
end

function transferfiles( files,confgfilename,ext )
    # ext=".RAW4"
    for dataFilename in files
        println( dataFilename)
        # tgrfilename = dirname( dataFilename )
        # basefilename = basename(dataFilename)
        extensition = splitext(dataFilename)
        tgr =  extensition[1] * ext
        tgr =  gettgrFilename(dataFilename)
        tgrfileObj = open(tgr,"w+")
		#if ispath( dataFilename)
        		Trasfer1File( dataFilename ,confgfilename,tgrfileObj)
		# else
		# 		close(tgrfileObj )
		# 		break
		# end
        close(tgrfileObj )
        # break
    end
end

function getmachineID()
    return Main.machineID
    return "Test-001"
end

function renamefiles( files )
    # ext=".RAW4"
    for dataFilename in files
        # println( dataFilename)
        if occursin("^", dataFilename)
            #println( dataFilename)
            continue
        end
        folder = dirname( dataFilename )
        basefilename = basename(dataFilename)
        filename,ext = splitext(basefilename)
        #println( folder," ", basefilename, " ",filename," ",ext," "  )
        fn = replace(filename,Pair("T","") )
        tgrname= "Data^Ver00^$(getmachineID())^"* fn*"^ID00"*ext
        fullname= joinpath( folder,tgrname)
        println( fullname )
        mv(dataFilename ,fullname  )
        # tgr =  extensition[1] * ext
        # tgrfileObj = open(tgr,"w+")
        # Trasfer1File( dataFilename ,confgfilename,tgrfileObj)
        # close(tgrfileObj )
        # break
    end
end

function gettgrFilename(srcname)
        if occursin("^", srcname)
            #println( dataFilename)
            return ""
        end
        folder = dirname( srcname )
        basefilename = basename(srcname)
        filename,ext = splitext(basefilename)
        #println( folder," ", basefilename, " ",filename," ",ext," "  )
        fn = replace(filename,Pair("T","") )
        ext =uppercase( ext )
        ext = replace(ext,Pair(".RAW",".RAW1"))
        ext = replace(ext,Pair(".LC",".RAW3"))
        ext = replace(ext,Pair(".ENERGY",".RAW4"))
        tgrname= "Data^Ver00^$(getmachineID())^"* fn*"^ID00"*ext
        fullname= joinpath( folder,tgrname)
        println( fullname )
        return fullname
        # mv(srcname ,fullname  )
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
module SpecialWrapper

using JavaCall

const BesselJ = @jimport org.hipparchus.special.BesselJ
const BesselJ_BesselJResult = @jimport org.hipparchus.special.BesselJ$BesselJResult
const Beta = @jimport org.hipparchus.special.Beta
const Class = @jimport java.lang.Class
const Erf = @jimport org.hipparchus.special.Erf
const Gamma = @jimport org.hipparchus.special.Gamma
const JString = @jimport java.lang.String
const Object = @jimport java.lang.Object

export BesselJ
export Beta
export Erf
export Gamma
export digamma
export equals
export erf
export erf_inv
export erfc
export erfc_inv
export gamma
export get_class
export hash_code
export inv_gamma1pm1
export lanczos
export log_beta
export log_gamma
export log_gamma1p
export notify
export notify_all
export regularized_beta
export regularized_gamma_p
export regularized_gamma_q
export rj_besl
export to_string
export trigamma
export value
export wait

include("bessel_j.jl")
include("beta.jl")
include("erf.jl")
include("gamma.jl")

end


text:
## Exercise 1-2
## Start the Julia REPL and use it as a calculator.

## 1. How many seconds are there in 42 minutes 42 seconds?
println("Ans 1:")

## 2. How many miles are there in 10 kilometers?
## TIP
## There are 1.61 kilometers in a mile.
println("Ans 2:")

## If you run a 10-kilometer race in 37 minutes 48 seconds, what is your average pace (time per mile in minutes and seconds)? What is your average speed in miles per hour?
println("Ans 3:")

println("End.")


text:
using Gadfly
using RDatasets

p1 = plot(dataset("datasets", "iris"), x="SepalLength", y="SepalWidth", Stat.binmean, Geom.point)
# color works
p2 = plot(dataset("datasets", "iris"), x="SepalLength", y="SepalWidth", color = "Species", Stat.binmean, Geom.point)
# n changes number of bins
p3 = plot(dataset("datasets", "iris"), x="SepalLength", y="SepalWidth", color="Species", Stat.binmean(n=3), Geom.point)
# integer works
p4 = plot(x= rand(1:10, 100), y = rand(1:10, 100), Stat.binmean, Geom.point)
#  surimpose regression
p5 = plot(
		layer(dataset("datasets", "iris"), x="SepalLength", y="SepalWidth", color="Species", Stat.binmean(n=5), Geom.point),
		layer(dataset("datasets", "iris"), x="SepalLength", y="SepalWidth", color="Species", Geom.smooth(method=:lm))
	)
vstack(p1, p2, p3, p4, p5)


text:
include("solver_preamble.jl")
using MosekTools
# It cannot be used in Direct mode yet as it gets:
# MathOptInterface.DeleteNotAllowed{MathOptInterface.ConstraintIndex{MathOptInterface.VectorOfVariables,MathOptInterface.RotatedSecondOrderCone}}
# for `sdsos_bivariate_quadratic`
factory = with_optimizer(Mosek.Optimizer, QUIET=true)
config = MOI.Test.TestConfig(atol=1e-5, rtol=1e-5)
@testset "Linear" begin
    Tests.linear_test(factory, config)
end
@testset "SOC" begin
    Tests.soc_test(factory, config)
end
@testset "SDP" begin
    Tests.sd_test(factory, config)
end


text:
function revers(side :: HorizonSide)
    side = HorizonSide((Int(side)+2)%4)
    return side
end

function movements_to_side!(r, side :: HorizonSide)

    while isborder(r,side) ==  false
        move!(r,side)
    end

end

function count_move!(r, side :: HorizonSide, count)
    for i = 1:count
        move!(r,side)
    end
end

function go_home(r, x, y, side_part :: HorizonSide, side :: HorizonSide)
    count_move!(r, side_part, y)
    count_move!(r, side, x)
end

function chess_move!(r, side :: HorizonSide, count)

    while isborder(r,side) ==  false
        if count % 2 == 0
            putmarker!(r)
        end

        move!(r,side)
        count = count + 1
    end

    if count % 2 == 0
        putmarker!(r)
    end

    return count
end

function part_chess(r, side_part :: HorizonSide, side :: HorizonSide)

    count = 0
    x = 0
    y = 0

    home = side
    flag = true

    while isborder(r,side_part) == false

        count = chess_move!(r, side, count)

        if flag
            x = count
            flag = false
        end

        move!(r, side_part)
        y = y + 1
        count = count + 1
        side = revers(side)

    end

    count = chess_move!(r, side, count)

    movements_to_side!(r, home)
    return x, y
end

function chess(r)

    x, y = part_chess(r, Sud, West)
    go_home(r, x, y, Nord, Ost)
    x, y = part_chess(r, Nord, Ost)
    go_home(r, x, y, Sud, West)

end

chess(r)
show(r)

max_stars_repo_name: AntonDudDm/mirea-progs

max_issues_repo_name: AntonDudDm/mirea-progs

max_forks_repo_name: AntonDudDm/mirea-progs


text:
include("mpctest.jl")
include("RunTests.jl")
include("../test/runtests.jl")
"""
test
"""
using Test

import Test.@test_skip, Test.@test_broken

# When testing the example solution, all tests must pass, even ones marked as skipped or broken.
# The track user will not be affected by this.
# Overwrite @test_skip, @test_broken with @test
macro test_skip(ex)
    @test eval(current_module(), ex)
end

macro test_broken(ex)
    @test eval(current_module(), ex)
end

# Store base path to cd back to it later
const path = pwd()

for exercise in readdir("exercises")
    # Allow only testing specified execises
    if !isempty(ARGS) && !(exercise in ARGS)
        continue
    end

    exercise_path = joinpath("exercises", exercise)
    isdir(exercise_path) || continue

    # Create temporary directory
    temp_path = mktempdir()

    # Copy tests & example to the temporary directory
    cp(joinpath(exercise_path, "example.jl"), joinpath(temp_path, "$exercise.jl"))
    cp(joinpath(exercise_path, "runtests.jl"), joinpath(temp_path, "runtests.jl"))

    try
        @info "Testing $exercise"
        # Run the tests within their own Julia process to prevent definitions from
        # one exercise leaking into another.
        cd(temp_path)
        run(`julia --color=yes runtests.jl`)
        cd(path)
        println()
    finally
        # Delete the temporary directory
        rm(temp_path, recursive=true)
    end
end

runtests()


text:
using StatsBase

abstract type AbstractStatistics end

abstract type StatsType end
struct GaussianStatsType <: StatsType end

mutable struct NonWeightedStatistics <: AbstractStatistics
  samples::Vector{Float64}
  isSorted::Bool
end

NonWeightedStatistics() = NonWeightedStatistics(Vector{Float64}(), false)

mutable struct GenericRiskStatistics{S <: StatsType} <: AbstractStatistics
  statsType::S
  samples::Vector{Float64}
  sampleWeights::StatsBase.Weights{Float64, Float64, Vector{Float64}}
  samplesMatrix::Matrix{Float64}
  isSorted::Bool
end

gen_RiskStatistics(dims::Int = 0) = GenericRiskStatistics(GaussianStatsType(), Vector{Float64}(undef, dims), weights(Vector{Float64}(undef, dims)), zeros(dims, 2), false)

const RiskStatistics = GenericRiskStatistics{GaussianStatsType}

mutable struct GenericSequenceStats{S <: AbstractStatistics} <: AbstractStatistics
  dimension::Int
  stats::Vector{S}
  results::Vector{Float64}
  quadraticSum::Matrix{Float64}
end

function GenericSequenceStats(dimension::Int, dim2::Int = dimension)
  # initial init
  gss = GenericSequenceStats(0, Vector{GenericRiskStatistics}(), Vector{Float64}(), Matrix{Float64}(undef, 0, 0))

  reset!(gss, dimension, dim2)

  return gss
end

GenericSequenceStats() = reset!(GenericSequenceStats(0, Vector{GenericRiskStatistics}(), Vector{Float64}(), Matrix{Float64}(undef, 0, 0)), 0)

function reset!(gss::GenericSequenceStats, dimension::Int, dim2::Int = dimension)
  # re-init
  if dimension > 0
    if dimension == gss.dimension
      for i in eachindex(gss.stats)
        reset!(gss.stats[i])
      end
    else
      gss.dimension = dimension
      gss.stats = GenericRiskStatistics[gen_RiskStatistics(dim2) for i = 1:dimension]
      gss.results = zeros(dimension)
    end
    gss.quadraticSum = zeros(dimension, dimension)
  else
    gss.dimension = dimension
  end

  return gss
end

function add_sample!(stat::GenericRiskStatistics, price::Float64, weight::Float64, idx::Int)
  # stat.samplesMatrix[idx, 1] = price
  # stat.samplesMatrix[idx, 2] = weight
  stat.samplesMatrix[idx, :] = [price weight]
  stat.isSorted = false
  return stat
end

add_sample!(stat::GenericRiskStatistics, price::Float64, idx::Int) = add_sample!(stat, price, 1.0, idx)

function add_sample!(stat::NonWeightedStatistics, price::Float64)
  push!(stat.samples, price)
  stat.isSorted = false
  return stat
end

function add_sample!(stat::GenericSequenceStats, vals::Vector, idx::Int, weight::Float64 = 1.0)
  if stat.dimension == 0
    # stat wasn't initialized
    dimension = length(vals)
    reset!(stat, dimension)
  end

  length(vals) == stat.dimension || error("sample size mismatch")

  BLAS.ger!(1.0, vals, vals, stat.quadraticSum)

  @simd for i in eachindex(stat.stats)
    @inbounds add_sample!(stat.stats[i], vals[i], weight, idx)
  end

  return stat
end

function adding_data!(stat::GenericRiskStatistics, sz::Int)
  # add to matrix
  stat.samplesMatrix = vcat(stat.samplesMatrix, zeros(sz, 2))
  # add to samples and sample weights
  append!(stat.samples, zeros(sz))
  stat.sampleWeights = weights(vcat(values(stat.sampleWeights), zeros(sz)))
  stat.isSorted = false
  return stat
end

function adding_data!(stat::GenericSequenceStats, sz::Int, sz2::Int)
  stat.dimension += sz
  append!(stat.results, zeros(sz))
  stat.quadraticSum = vcat(stat.quadraticSum, zeros(sz, 2))

  # now we have to add the additional stats
  newStats = GenericRiskStatistics[gen_RiskStatistics(sz2) for i = 1:sz]

  append!(stat.stats, newStats)

  return stat
end

function sort_samples!(stat::GenericRiskStatistics)
  if ~stat.isSorted
    stat.samplesMatrix = sortslices(stat.samplesMatrix, dims=1)
    stat.samples = stat.samplesMatrix[:, 1]
    stat.sampleWeights = weights(stat.samplesMatrix[:, 2])
    stat.isSorted = true
  end
  return stat
end

function sort_samples!(stat::NonWeightedStatistics)
  sort!(stat.samples)
  stat.isSorted = true
end

function error_estimate(stat::AbstractStatistics)
  sort_samples!(stat)
  return sqrt(var(stat.samples)/ length(stat.samples))
end

function weight_sum(stat::AbstractStatistics)
  sort_samples!(stat)
  return sum(stat.sampleWeights)
end

weight_sum(stat::GenericSequenceStats) = length(stat.stats) == 0 ? 0.0 : weight_sum(stat.stats[1])

function sample_num(stat::AbstractStatistics)
  sort_samples!(stat)

  return length(stat.samples)
end

sample_num(stat::GenericSequenceStats) = length(stat.stats) == 0 ? 0 : sample_num(stat.stats[1])

function stats_mean(stat::AbstractStatistics)
  sort_samples!(stat)
  return mean(stat.samples)
end

function stats_mean(stat::GenericSequenceStats)
  return Float64[stats_mean(x) for x in stat.stats]
end

function stats_std_deviation(stat::AbstractStatistics)
  sort_samples!(stat)
  return std(stat.samples)
end

function stats_skewness(stat::AbstractStatistics)
  sort_samples!(stat)
  return skewness(stat.samples)
end

function stats_kurtosis(stat::AbstractStatistics)
  sort_samples!(stat)
  return kurtosis(stat.samples)
end

function stats_covariance(stat::GenericSequenceStats)
  sampleWeight = weight_sum(stat)
  sampleWeight > 0.0 || error("sample weight of 0 insufficient")

  sampleNumber = float(sample_num(stat))
  sampleNumber > 1.0 || error("sample number <= 1.0 is insufficient")

  m = stats_mean(stat)
  inv = 1.0 / sampleWeight

  result = inv * stat.quadraticSum

  BLAS.ger!(-1.0, m, m, result)

  result *= (sampleNumber / (sampleNumber - 1.0))

  return result
end