text:
### BasicDiscParamIOStream

mutable struct BasicDiscParamIOStream <: ParameterIOStream{Discrete}
  value::Union{IOStream, Void}
  loglikelihood::Union{IOStream, Void}
  logprior::Union{IOStream, Void}
  logtarget::Union{IOStream, Void}
  diagnosticvalues::Union{IOStream, Void}
  names::Vector{AbstractString}
  size::Tuple
  n::Integer
  diagnostickeys::Vector{Symbol}

  function BasicDiscParamIOStream(
    size::Tuple,
    n::Integer,
    streams::Vector{Union{IOStream, Void}},
    diagnostickeys::Vector{Symbol}=Symbol[],
    filenames::Vector{AbstractString}=[(streams[i] == nothing) ? "" : streams[i].name[7:end-1] for i in 1:5]
  )
    instance = new()

    fnames = fieldnames(BasicDiscParamIOStream)
    for i in 1:5
      setfield!(instance, fnames[i], streams[i])
    end

    instance.names = filenames

    instance.size = size
    instance.n = n
    instance.diagnostickeys = diagnostickeys

    instance
  end
end

function BasicDiscParamIOStream(
  size::Tuple,
  n::Integer,
  filenames::Vector{AbstractString},
  diagnostickeys::Vector{Symbol}=Symbol[],
  mode::AbstractString="w"
)
  fnames = fieldnames(BasicDiscParamIOStream)
  BasicDiscParamIOStream(
    size,
    n,
    [isempty(filenames[i]) ? nothing : open(filenames[i], mode) for i in 1:5],
    diagnostickeys,
    filenames
  )
end

function BasicDiscParamIOStream(
  size::Tuple,
  n::Integer;
  monitor::Vector{Bool}=[true; fill(false, 3)],
  filepath::AbstractString="",
  filesuffix::AbstractString="csv",
  diagnostickeys::Vector{Symbol}=Symbol[],
  mode::AbstractString="w"
)
  fnames = fieldnames(BasicDiscParamIOStream)

  filenames = Array{AbstractString}(5)
  for i in 1:4
    filenames[i] = (monitor[i] == false ? "" : joinpath(filepath, string(fnames[i])*"."*filesuffix))
  end
  filenames[5] = (isempty(diagnostickeys) ? "" : joinpath(filepath, "diagnosticvalues."*filesuffix))

  BasicDiscParamIOStream(size, n, filenames, diagnostickeys, mode)
end

function BasicDiscParamIOStream(
  size::Tuple,
  n::Integer,
  monitor::Vector{Symbol};
  filepath::AbstractString="",
  filesuffix::AbstractString="csv",
  diagnostickeys::Vector{Symbol}=Symbol[],
  mode::AbstractString="w"
)
  fnames = fieldnames(BasicDiscParamIOStream)
  BasicDiscParamIOStream(
    size,
    n,
    monitor=[fnames[i] in monitor ? true : false for i in 1:4],
    filepath=filepath,
    filesuffix=filesuffix,
    diagnostickeys=diagnostickeys,
    mode=mode
  )
end

function open(iostream::BasicDiscParamIOStream, filenames::Vector{S}, mode::AbstractString="w") where {S<:AbstractString}
  iostream.names = filenames
  fnames = fieldnames(BasicDiscParamIOStream)

  for i in 1:5
    if getfield(iostream, fnames[i]) != nothing
      setfield!(iostream, fnames[i], open(filenames[i], mode))
    end
  end
end

function close(iostream::BasicDiscParamIOStream)
  fnames = fieldnames(BasicDiscParamIOStream)
  for i in 1:5
    if getfield(iostream, fnames[i]) != nothing
      close(getfield(iostream, fnames[i]))
    end
  end
end

function mark(iostream::BasicDiscParamIOStream)
  fnames = fieldnames(BasicDiscParamIOStream)
  for i in 1:5
    if getfield(iostream, fnames[i]) != nothing
      mark(getfield(iostream, fnames[i]))
    end
  end
end

function reset(iostream::BasicDiscParamIOStream)
  fnames = fieldnames(BasicDiscParamIOStream)
  for i in 1:5
    if getfield(iostream, fnames[i]) != nothing
      reset(getfield(iostream, fnames[i]))
    end
  end
end

function flush(iostream::BasicDiscParamIOStream)
  fnames = fieldnames(BasicDiscParamIOStream)
  for i in 1:5
    if getfield(iostream, fnames[i]) != nothing
      flush(getfield(iostream, fnames[i]))
    end
  end
end

function write(iostream::BasicDiscParamIOStream, state::ParameterState{Discrete, F}) where {F<:VariateForm}
  fnames = fieldnames(BasicDiscParamIOStream)
  for i in 1:5
    if getfield(iostream, fnames[i]) != nothing
      write(getfield(iostream, fnames[i]), join(getfield(state, fnames[i]), ','), "\n")
    end
  end
end

function write(iostream::BasicDiscParamIOStream, nstate::BasicDiscUnvParameterNState)
  fnames = fieldnames(BasicDiscParamIOStream)
  for i in 1:4
    if getfield(iostream, fnames[i]) != nothing
      writedlm(getfield(iostream, fnames[i]), nstate.(fnames[i]))
    end
  end
  if iostream.diagnosticvalues != nothing
    writedlm(iostream.diagnosticvalues, nstate.diagnosticvalues', ',')
  end
end

function write(iostream::BasicDiscParamIOStream, nstate::BasicDiscMuvParameterNState)
  fnames = fieldnames(BasicDiscParamIOStream)
  for i in 2:4
    if getfield(iostream, fnames[i]) != nothing
      writedlm(getfield(iostream, fnames[i]), nstate.(fnames[i]))
    end
  end
  for i in (1, 5)
    if getfield(iostream, fnames[i]) != nothing
      writedlm(getfield(iostream, fnames[i]), nstate.(fnames[i])', ',')
    end
  end
end

function read!(iostream::BasicDiscParamIOStream, nstate::BasicDiscUnvParameterNState{NI, NR}) where {NI<:Integer, NR<:Real}
  t = [NI, fill(NR, 3)]
  fnames = fieldnames(BasicDiscParamIOStream)
  for i in 1:4
    if getfield(iostream, fnames[i]) != nothing
      setfield!(nstate, fnames[i], vec(readdlm(getfield(iostream, fnames[i]), ',', t[i])))
    end
  end
  if iostream.diagnosticvalues != nothing
    nstate.diagnosticvalues = readdlm(iostream.diagnosticvalues, ',', Any)'
  end
end

function read!(iostream::BasicDiscParamIOStream, nstate::BasicDiscMuvParameterNState{NI, NR}) where {NI<:Integer, NR<:Real}
  fnames = fieldnames(BasicDiscParamIOStream)
  if getfield(iostream, fnames[1]) != nothing
    setfield!(nstate, fnames[1], readdlm(getfield(iostream, fnames[1]), ',', NI)')
  end
  for i in 2:4
    if getfield(iostream, fnames[i]) != nothing
      setfield!(nstate, fnames[i], vec(readdlm(getfield(iostream, fnames[i]), ',', NR)))
    end
  end
  if iostream.diagnosticvalues != nothing
    nstate.diagnosticvalues = readdlm(iostream.diagnosticvalues, ',', Any)'
  end
end

function read(iostream::BasicDiscParamIOStream, TI::Type{NI}, TR::Type{NR}) where {NI<:Integer, NR<:Real}
  local nstate::DiscreteParameterNState
  fnames = fieldnames(BasicDiscParamIOStream)
  l = length(iostream.size)

  if l == 0
    nstate = BasicDiscUnvParameterNState(
      iostream.n,
      [getfield(iostream, fnames[i]) != nothing ? true : false for i in 1:13],
      iostream.diagnostickeys,
      TI,
      TR
    )
  elseif l == 1
    nstate = BasicDiscMuvParameterNState(
      iostream.size[1],
      iostream.n,
      [getfield(iostream, fnames[i]) != nothing ? true : false for i in 1:13],
      iostream.diagnostickeys,
      TI,
      TR
    )
  else
    error("BasicDiscParamIOStream.size must be a tuple of length 0 or 1, got $(iostream.size) length")
  end

  read!(iostream, nstate)

  nstate
end

function show(io::IO, iostream::BasicDiscParamIOStream)
  fnames = fieldnames(BasicDiscParamIOStream)
  fbool = map(n -> getfield(iostream, n) != nothing, fnames[1:4])
  indentation = "  "

  println(io, "BasicDiscParamIOStream:")

  println(io, indentation*"state size = $(iostream.size)")
  println(io, indentation*"number of states = $(iostream.n)")

  print(io, indentation*"monitored components:")
  if !any(fbool)
    println(io, " none")
  else
    print(io, "\n")
    for i in 1:4
      if fbool[i]
        println(io, string(indentation^2, fnames[i]))
      end
    end
  end

  print(io, indentation*"diagnostics:")
  if isempty(iostream.diagnostickeys)
    print(io, " none")
  else
    for k in iostream.diagnostickeys
      print(io, "\n")
      print(io, string(indentation^2, k))
    end
  end
end


text:
"""
    Find the price of capital if the firm is not selling everything it produces.
"""
function find_price_info_unc(par)
    f(p_I) = begin
        @set! par.p_I = p_I^2
        k_I_b_A, k_I_b_B = solve_firm_problems_unc(par)
        diff = (1 - par.γ) * k_I_b_B + (1 - par.τ) * par.γ * k_I_b_A
        return diff # What one firm sells needs to be equal to what the other firm buys
    end
    find_zero(f, 0.25, tol = 1e-16)^2
end

function solve_firm_problems_unc(par)
    k_I_P_A = k_I_prod_A(par)
    if par.ϕ > 0 && loss_A(-k_I_P_A, par) > 0
        k_I_b_A = -k_I_P_A
    else
        f_A(x) = loss_A(x, par)
        k_I_b_A = find_zero(f_A, [-k_I_P_A + 0.00001, 5])
    end

    # B is always unconstrained, because they are buying capital
    k_I_P_B = k_I_prod_B(par)
    f_B(x) = loss_B(x, par)
    k_I_b_B = find_zero(f_B, [-k_I_P_B + 0.0001, 5])

    return k_I_b_A, k_I_b_B
end

"""
    This is much more robust!
"""
function solve_agg_C_iterative_unc(par; tol = 1e-10)
    diff = 10
    C_guess = par.C
    w_guess = par.w
    while diff > tol
        @set! par.C = C_guess
        #@show par.C
        C_sol, w_sol = find_agg_C_unc(par)
        @set! par.w = w_sol
        diff = copy((C_sol - C_guess)^2 + (w_sol - w_guess)^2)
        if diff > tol
            C_guess = copy(C_sol)
            w_guess = copy(w_sol)
        else
            @set! par.w = w_sol
            @set! par.C = C_guess
            @set! par.p_I = find_price_info_unc(par)
            @set! par.k_I_P_A_ss = k_I_prod_A(par)
            @set! par.k_I_P_B_ss = k_I_prod_B(par)

            f_A(x) = loss_A(x, par)
            @set! par.k_I_b_A_ss = find_zero(f_A, [-par.k_I_P_A_ss + 0.00001, 5])

            f_B(x) = loss_B(x, par)
            @set! par.k_I_b_B_ss = find_zero(f_B, [-par.k_I_P_B_ss + 0.00001, 5])

            @set! par.k_I_A_ss = total_capital_A(par.k_I_b_A_ss, par.k_I_P_A_ss, par)
            @set! par.k_I_B_ss = total_capital_B(par.k_I_b_B_ss, par.k_I_P_B_ss)

            @set! par.l_A_ss = labor_demand_A(par.k_I_A_ss, par)
            @set! par.l_B_ss = labor_demand_B(par.k_I_B_ss, par)
        end
    end
    return par
end

function find_agg_C_unc(par)
    @set! par.p_I = find_price_info_unc(par)
    k_I_b_A, k_I_b_B = solve_firm_problems_unc(par)

    k_I_A, l_A = capital_and_labor_unc_A(k_I_b_A, par)
    k_I_B, l_B = capital_and_labor_B(k_I_b_B, par)

    C_A = firm_prod(l_A, k_I_A; par)
    C_B = firm_prod(l_B, k_I_B; par)

    return agg_prod(C_A, C_B, par), wage_rate(k_I_A, k_I_B; par)
end

"""
    Profit for both type of firms.
"""
function profit_unc(par)
    @unpack_Pars_Simple par

    profit_A = γ * v * C^α_C * l_A_ss^α_hat * k_I_A_ss^δ_hat - w * l_A_ss - p_I * k_I_b_A_ss - c_A(k_I_P_A_ss, par)
    profit_B = (1 - γ) * v * C^α_C * l_B_ss^α_hat * k_I_B_ss^δ_hat - w * l_B_ss - p_I / (1 - τ) * k_I_b_B_ss - c_B(k_I_P_B_ss, par)

    return profit_A, profit_B
end

#! Now, go to the constrained problem, but both firms can sell and buy!
"""
    Solve the firms problem if
"""
function solve_firm_problems_con(par)
    k_I_P_A = k_I_prod_A(par)
    if par.ϕ > 0 && loss_A(-k_I_P_A, par) > 0
        k_I_b_A = -k_I_P_A
    else
        f_A(x) = loss_A(x, par)
        k_I_b_A = find_zero(f_A, [-k_I_P_A + 0.00001, 5])
    end

    # B is always unconstrained, because they are buying capital
    k_I_P_B = k_I_prod_B(par)
    f_B(x) = loss_B(x, par)
    k_I_b_B = find_zero(f_B, [-k_I_P_B + 0.0001, 5])

    return k_I_b_A, k_I_b_B
end

ext: jl


text:
l=length(ARGS)

buildscript="BuildExecutable/src/build_executable.jl"

default_version=string(VERSION)

version=l<1?default_version:ARGS[1]
name=l<2?"app":ARGS[2]
dir=l<3?"build":ARGS[3]
script=l<4?"main.jl":ARGS[4]

debug = false
julia = replace(abspath(joinpath(JULIA_HOME, debug ? "julia-debug" : "julia")), default_version,version)

path = replace(dirname(Base.source_path()),"\\","/")

targetdir_build=joinpath(path,dir)
targetdir_version=joinpath(targetdir_build,version)
targetdir_project=joinpath(targetdir_version,name)
targetdir = targetdir_project

sourcescript=joinpath(path,script)
buildscriptpath=buildscript

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

if !isfile(julia * (is_windows() ? ".exe" : "")) error("'$julia' not found.") end

# find buildscript...
if !isfile(buildscriptpath)
  buildscriptpath=joinpath(path,buildscript)
  if !isfile(buildscriptpath)
    current = path
    prev = current
    found=false
    for i=1:10
      current=abspath(joinpath(current,"../"))
      buildscriptpath=joinpath(current,buildscript)
      if isfile(buildscriptpath) found=true; break end
      if current == prev break end #dublicate, sp break
      prev = current
    end
  end
  if !found error("'$buildscript' not found.") end
end

if !isfile(sourcescript) error("'$sourcescript' not found.") end

if !isdir(targetdir_build) mkdir(targetdir_build) end
if !isdir(targetdir_version) mkdir(targetdir_version) end
if !isdir(targetdir_project) mkdir(targetdir_project) end

if !isdir(targetdir_build) error("'$targetdir_build' failed to create.") end
if !isdir(targetdir_version) error("'$targetdir_version' failed to create.") end
if !isdir(targetdir_project) error("'$targetdir_project' failed to create.") end

cmd=`$julia $buildscriptpath --force "$name" $sourcescript $targetdir`

info(cmd)
run(cmd)

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
using Roots

root = fzero(sin, -0.5, 0.5)
println(sin(root))

f(x) = (x - 1)^3
root = fzero(f, 0.5, 2.0)
println(f(root))

f(x) = 2x*exp(-20) - 2*exp(-20x) + 1

roots = fzero(f, 0.0, 1.0; tolerance=1e-10, max_iter=100)
println(f(root))

f(x) = exp(x) - cos(x)
fp(x) = exp(x) + sin(x)

root = newton(f, fp, 3.0)
println(f(root))

f(x) = exp(x) - cos(x)
fp(x) = exp(x) + sin(x)
fpp(x) = exp(x) + cos(x)

root = halley(f, fp, fpp, 3.0)
println(f(root))

using Calculus

f(x) = exp(x) - cos(x)
root = newton(f, derivative(f), 3.0)
println(f(root))


text:
using RemoteS, GMT, Dates
using Test, Statistics, Printf

include("../src/sat_tracks.jl")		# Even without SatelliteToolbox loaded we can still test some functions

truecolor(mat2img(rand(UInt16,128,128)), mat2img(rand(UInt16,128,128)), mat2img(rand(UInt16,128,128)));
@test RemoteS.guess_increment_from_coordvecs([1., 1, 1, 1], [1., 1, 1, 1]) == [1.0, 1.0]
@test RemoteS.helper_find_sds("AA", "xxxxxxxx:AA", findall("\n", @sprintf("aA\nbbbbbnnn\n"))) == "xxx:AA"

truecolor(mat2img(rand(UInt16, 16, 16, 3), noconv=1), [1,2,3]);
truecolor(mat2img(rand(UInt16, 16, 16), noconv=1), mat2img(rand(UInt16, 16, 16), noconv=1), mat2img(rand(UInt16, 16, 16), noconv=1));

@test RemoteS.assign_description(["LC08_B2.TIF", "LC08_B5.TIF"], String[])[1] == "Band 2 - Blue [0.45-0.51]"
@test RemoteS.assign_description(["AC08_B2.TIF", "AC08_B5.TIF"], String[])  == ["AC08_B2", "AC08_B5"]
@test RemoteS.assign_description(["AC08_B2.TIF", "AC08_B5.TIF"], ["B2", "B5"]) == ["B2", "B5"]

@test reportbands(mat2img(rand(UInt16, 4,4,3), names=["Band 1", "Band 2", "Band 3"]), 3)[1] == "Band 3"
@test reportbands(mat2img(rand(UInt16, 4,4,3), names=["Band 1", "Band 2", "Band 3"]), layers=[1,3]) == ["Band 1", "Band 3"]

#sat_tracks(position=true);
#sat_tracks(position=true, tle=["1 27424U 02022A   21245.83760660  .00000135  00000-0  39999-4 0  9997"; "2 27424  98.2123 186.0654 0002229  67.6025 313.3829 14.57107527 28342"]);

@test_throws ErrorException("Bad input type Symbol. Must be a DateTime, a String or a Tuple(Int)") getitDTime(:val)
@test get_sat_name(Dict(:sat => "TERRA")) == "TERRA"
@test get_MODIS_scene_name(datetime2julian(DateTime("2020-09-20")), "A") == "AQUA_MODIS.20200920T000001.L2.SST.NRT.nc"
@test get_MODIS_scene_name(datetime2julian(DateTime("2020-09-20")), "A", false) == "A2020264000000.L2_LAC_OC.nc"

	println("	reading cube")
cube = gmtread("LC08_cube.tiff", layout="TRB");
truecolor("LC08_cube.tiff");	println("	truecolor 1")
truecolor(cube, [3,2,1]);		println("	truecolor 2")
truecolor(cube, layers=[3,2,1]); println("	truecolor 3")
	println("	reportbands 1")
@test reportbands("LC08_cube.tiff", 2)[1] == "Band 3 - Green [0.53-0.59]"
	println("	reportbands 2")
reportbands(cube);
	println("	subcube")
Irgb = subcube("LC08_cube.tiff", bandnames = ["red", "green", "blue"]);
RemoteS.parse_lsat8_file("LC08_cube.tiff", band=10);
dn2temperature("LC08_cube.tiff", band=10);	println("	dn2temperature")
dn2radiance("LC08_cube.tiff", band=5);		println("	dn2radiance")
dn2reflectance("LC08_cube.tiff", band=5);	println("	dn2reflectance")
reflectance_surf("LC08_cube.tiff", band=5); println("	reflectance_surf")

R = subcube("LC08_cube.tiff", bandnames=["red"]);	println("	subcube R")
G = subcube("LC08_cube.tiff", bandnames=["green"]);	println("	subcube G")
B = subcube("LC08_cube.tiff", bandnames=["blue"]);	println("	subcube B")
NIR = subcube("LC08_cube.tiff", bandnames=["nir"]);	println("	subcube NIR")
SW1 = subcube("LC08_cube.tiff", bandnames=["swir 1"]);	println("	subcube SW1")
SW2 = subcube("LC08_cube.tiff", bandnames=["swir 2"]);	println("	subcube SW2")
G1 = evi("LC08_cube.tiff"); println("	evi 1")
G2 = evi(B, R, NIR);		println("	evi 2")
G3 = evi(cube,[1,3,4]);		println("	evi 3")
@test G1.range[5:6] == G2.range[5:6]
G1 = evi2("LC08_cube.tiff"); println("	evi2 1")
G2 = evi2(R, NIR);			println("	evi2 2")
@test G1.range[5:6] == G2.range[5:6]
G1 = msavi("LC08_cube.tiff"); println("	msavi 1")
G2 = msavi(R, NIR);			println("	msavi 2")
@test G1.range[5:6] == G2.range[5:6]
G1 = ndvi("LC08_cube.tiff"); println("	ndvi 1")
G2 = ndvi(R, NIR);			println("	ndvi 2")
@test G1.range[5:6] == G2.range[5:6]
G1 = satvi("LC08_cube.tiff"); println("	satvi 1")
G2 = satvi(R, SW1, SW2);	println("	satvi 2")
@test G1.range[5:6] == G2.range[5:6]
G1 = savi("LC08_cube.tiff"); println("	savi 1")
G2 = savi(R, NIR);			println("	savi 2")
@test G1.range[5:6] == G2.range[5:6]
G1 = slavi("LC08_cube.tiff"); println("	slavi 1")
G2 = slavi(R, NIR, SW2);	println("	slavi 2")
@test G1.range[5:6] == G2.range[5:6]

max_stars_repo_name: GenericMappingTools/RemoteS

max_issues_repo_name: GenericMappingTools/RemoteS

max_forks_repo_name: GenericMappingTools/RemoteS


text:
const td = mktempdir()
J.FOLDER_PATH[] = td

J.def_GLOBAL_PAGE_VARS!()
J.def_GLOBAL_LXDEFS!()

@testset "Paths" begin
    P = J.set_paths!()

    @test J.PATHS[:folder]       == td
    @test J.PATHS[:src]      == joinpath(td, "src")
    @test J.PATHS[:src_css]  == joinpath(td, "src", "_css")
    @test J.PATHS[:src_html] == joinpath(td, "src", "_html_parts")
    @test J.PATHS[:libs]    == joinpath(td, "libs")
    @test J.PATHS[:pub]     == joinpath(td, "pub")
    @test J.PATHS[:css] == joinpath(td, "css")

    @test P == J.PATHS

    mkdir(J.PATHS[:src])
    mkdir(J.PATHS[:src_pages])
    mkdir(J.PATHS[:libs])
    mkdir(J.PATHS[:src_css])
    mkdir(J.PATHS[:src_html])
    mkdir(J.PATHS[:assets])
end

# copying _libs/katex in the J.PATHS[:libs] so that it can be used in testing
# the js_prerender_math
cp(joinpath(dirname(dirname(pathof(JuDoc))), "test", "_libs", "katex"), joinpath(J.PATHS[:libs], "katex"))

@testset "Set vars" begin
    d = J.PageVars(
    	"a" => 0.5 => (Real,),
    	"b" => "hello" => (String, Nothing))
    J.set_vars!(d, ["a"=>"5", "b"=>"nothing"])

    @test d["a"].first == 5
    @test d["b"].first === nothing

    @test_logs (:warn, "Doc var 'a' (type(s): (Real,)) can't be set to value 'blah' (type: String). Assignment ignored.") J.set_vars!(d, ["a"=>"\"blah\""])
    @test_logs (:error, "I got an error (of type 'DomainError') trying to evaluate '__tmp__ = sqrt(-1)', fix the assignment.") J.set_vars!(d, ["a"=> "sqrt(-1)"])

    # assigning new variables

    J.set_vars!(d, ["blah"=>"1"])
    @test d["blah"].first == 1
end

@testset "Def+coms" begin # see #78
    st = raw"""
        @def title = "blah" <!-- comment -->
        @def hasmath = false
        etc
        """ * J.EOS
    (m, jdv) = J.convert_md(st)
    @test jdv["title"].first == "blah"
    @test jdv["hasmath"].first == false
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ ebcb6ec6-61f9-11eb-05a3-817618350415
using PlutoUI

# ╔═╡ 2d1693e6-61fa-11eb-3cb6-9f5431449a41
begin
	using Plots
	pyplot()
end

# ╔═╡ f5d3f88c-620e-11eb-3bc8-e3e0e1598482
using Random

# ╔═╡ 6b69c792-62bd-11eb-0f1d-2d0081598dfa
using Distributions

# ╔═╡ b80b8e58-62bf-11eb-11fb-1397b06ae3e1
using QuadGK

# ╔═╡ ebd217ba-61f8-11eb-2337-9db63173ed34
md"# Module 2: _Probability Quick Refresher in Julia_
**CS6741: Statistical Foundations of Data Science**

In this module, we will quickly revise some aspects of Probability through exercises on Julia.

"

# ╔═╡ 4e0bacac-621c-11eb-17be-a968500dc6bc
md"### Dice Throws"

# ╔═╡ 48a0d2ac-6a1d-11eb-2441-897ca192e8ad
md"👉 How would you simulate a fair dice?"

# ╔═╡ 29fb09a2-61f9-11eb-0699-4dbf7daca0d7
faces = 1:6

# ╔═╡ bbae804a-61f9-11eb-03fe-89049682eae7
rand(faces)

# ╔═╡ e977219e-620d-11eb-2ca1-2141cc7c7d2a
md"What is the distribution of values obtained with throws?"

# ╔═╡ ee1bd118-61f9-11eb-1a18-df7d29a96c97
@bind throwdice Button("Throw dice")

# ╔═╡ 059c87a6-61fa-11eb-01b4-09c40af1be94
throws = []

# ╔═╡ 0a10f34e-61fa-11eb-365d-01ff16513530
begin
	throwdice
	for _ in 1:1000
		push!(throws, rand(faces))
	end
end

# ╔═╡ e00f4d02-61f9-11eb-2c5e-19ab03f4bdeb
begin
	throwdice
	histogram(throws, bins = [1, 2, 3, 4, 5, 6, 7], normalize=true, legend=false)
end

# ╔═╡ 20ea39f4-620e-11eb-05b3-716b88bcd205
md"👉 What about the distribution of throws of two dice?"

# ╔═╡ c0b6a868-61fb-11eb-2190-83d6d760ede4
begin
	N = 1000
	two_throws =  [(rand(faces), rand(faces)) for _ in 1:N]
end

# ╔═╡ 947c880a-6203-11eb-1bfe-95c8e630b395
begin
	freq_matrix = zeros(length(faces), length(faces))
	for i in faces
		for j in faces
			freq_matrix[i, j] = count(two_throws .== [(i, j)]) / N
		end
	end
end

# ╔═╡ fb898100-6204-11eb-1faf-cd75283f54d4
freq_matrix

# ╔═╡ b25c21be-6203-11eb-2f40-5db01b4beb89
heatmap(faces, faces, freq_matrix, clim=(0, 0.1))

# ╔═╡ 2b3a4f0c-620e-11eb-2021-bdaab9318fc1
md"👉 What is the probability that the sum of two dice is prime?"

# ╔═╡ c2ae2510-620a-11eb-07e0-91cf24657a55
function isPrime(x)
	if x == 1
		return false
	elseif x == 2
		return true
	end
	for i = 2:sqrt(x)
		if x%i == 0
			return false
		end
	end
	return true
end

# ╔═╡ 853a1ab6-6a1d-11eb-2e08-93e5bcc73c7c
isPrime.(2:12)

# ╔═╡ 33d3be8c-620e-11eb-3f7d-8f9a60cfdefc
md"We have 11 unique sums possible 2 through 12 out of which only 5 (2, 3, 5, 7, 11) are prime. Hence probability is:"

# ╔═╡ 6205d7aa-620d-11eb-3c75-83b491a84995
count(isPrime.(2:12)) / 11

# ╔═╡ 6f0eec76-6a1d-11eb-2585-17959367c09f
md"👉 Is the above correct? Why or why not?"

# ╔═╡ ba50d33e-6a1d-11eb-0687-47cb0b6f8242
md"👉 What visualiation would you use to help you with this?"

# ╔═╡ ec8b9cd0-6211-11eb-06fa-37880ae43e40
histogram([(rand(faces) + rand(faces)) for _ in 1:10000])

# ╔═╡ e37eefb6-6a1d-11eb-000b-89dbddabce5d
md"👉 How do we fix this?"

# ╔═╡ 5b2accf0-620e-11eb-06b9-ef354ebb3261
md"Correct theoretical approach assumes equal probability of the individual dice throws and not the sum"

# ╔═╡ 7dced6ec-6205-11eb-3e95-5f2fa0d5adf3
theorySol = sum([isPrime(i + j) for i in faces, j in faces]) / length(faces)^2

# ╔═╡ 7bbe53a6-620e-11eb-296a-dbf3eaab558d
md"Now we will compute the same thing by running lots of experiments with randomly sampled data. Usually called _Monte Carlo experiments_"

# ╔═╡ 014acff6-6206-11eb-0f3f-f17ef6edf8dd
@bind n_exp html"<input type=range min=1 max=6>"

# ╔═╡ 987d2bb2-620d-11eb-07b6-1736ac133029
n_exp

# ╔═╡ 0c205d04-6205-11eb-2d31-5f689e25a3f1
est = sum([isPrime(rand(faces) + rand(faces)) for i in 1:10^n_exp]) / 10^n_exp

# ╔═╡ beccfb34-620e-11eb-177a-711b597f5f6f
md"These random experiments give different results each time"

# ╔═╡ b0c7733e-620e-11eb-0fb6-f3a9f7631823
[sum([isPrime(rand(faces) + rand(faces)) for i in 1:10^5]) / 10^5 for j = 1: 10]

# ╔═╡ ce084cfc-620e-11eb-2881-8f30a7558bba
md"For reproducibility, we can set the seed of the random number generator"

# ╔═╡ d8b7f274-620e-11eb-13b1-3933af0f7157
begin
	Random.seed!(0)
	sum([isPrime(rand(faces) + rand(faces)) for i in 1:10^5]) / 10^5
end

# ╔═╡ 77fb84bc-621c-11eb-2cc3-71425a74f823
md"Notice that in the above we used the probability concepts of _independence_ of two random variables and _Monte Carlo_ simulations."

# ╔═╡ 59437fd2-621c-11eb-1e31-f1973639a0e2
md"### Picking balls"

# ╔═╡ 8fbbc5a6-6212-11eb-1020-7382ef578cea
md"Selection with and without replacement"

# ╔═╡ beeacf18-6212-11eb-1929-a5c830dddbcb
md"Let us suppose a bag contains 5 balls - 3 tennis balls, 2 cricket balls"

# ╔═╡ 6aa6b4ec-6212-11eb-39ea-db9d36b22cce
bag = ["🎾", "🥎", "🎾", "🎾", "🥎"]

# ╔═╡ bc4822e0-6212-11eb-3a49-e5307fe49156
rand(bag)

# ╔═╡ d5215c20-6212-11eb-32f5-c71343b60d22
md"👉 What is the probability that if I draw two balls both are tennis balls"

# ╔═╡ e6f2b97e-6212-11eb-0ef4-6593e840b93a
md"The easier case is with replacement. It is simply $P(🎾) \times P(🎾)$"

# ╔═╡ e49d22cc-6212-11eb-2732-27ad5e863ffc
3/5 * 3/5

# ╔═╡ f462afb0-6212-11eb-1a80-239d798e25f8
md"👉 What if we do not allow replacement?"

# ╔═╡ 383c0f2e-6213-11eb-355e-5b16663f1897
md"We are looking for $P(🎾|\phi) \times P(🎾|🎾)$"

# ╔═╡ 559fea86-6213-11eb-185e-b5fc917251fe
3/5 * 2/4

# ╔═╡ 73f6494e-6213-11eb-2a42-7b58fe47d597
md"👉 How do we implement this in Julia when running for more complex examples? Which one do you think is easier to implement - with or without replacements?"

# ╔═╡ 1ad550e8-6214-11eb-108b-0d6e08518018
md"With replacement:"

# ╔═╡ e455c886-6213-11eb-01b8-5df9e2bb7ea8
sum([(rand(bag) == "🎾" && rand(bag) == "🎾") for i in 1:100000])/100000

# ╔═╡ 33bc8ece-6a1f-11eb-20c9-fd32413cc2ac
md"👉 How would you implement this without replacement?"

# ╔═╡ 17ebe7ac-6214-11eb-38f1-3b757972e958
md"Without replacement:"

# ╔═╡ 3d281fc0-6214-11eb-16ae-89900d7f437f
begin
	num = 10000
	count_ = 0
	for i = 1:num
		bag1 = copy(bag)
		ball1 = rand(bag1)
		deleteat!(bag1, findfirst(x-> x == "🎾", bag1))
		ball2 = rand(bag1)
		if ball1 == "🎾" && ball2 == "🎾"
			count_ += 1
		end
	end
	count_ / num
end

# ╔═╡ 94d49374-621c-11eb-2985-b5ee5644c67e
md"Notice that in the above example we use the idea of _conditional probability_"

# ╔═╡ 5eff3d3a-621c-11eb-0824-f76604aed501
md"### Inclusion exclusion principle"

# ╔═╡ ba16f60e-623f-11eb-190a-af4758c32068
md"You have to inform $n$ people about some important news, but you have a fancy quantum phone that randomly dials a person. If you are allowed to make $r \geq n$ calls, what is the probability that all $n$ people are informed of the important news."

# ╔═╡ 51413a62-6a1f-11eb-1400-5918c088122c
md"👉 How would you do this? Think first of the mathematical approach, and then we will look at the computational approach."

# ╔═╡ eed974c0-623f-11eb-021e-ebfc3c85b934
md"Let $A_i$ denote the event that the $i$th person is not informed. Then the probability we seek is
```math
P(B) = 1 - P(A_1 \cup \ldots \cup A_n)
```"

# ╔═╡ 2cf52d58-6240-11eb-1e21-f3c0a71b4525
md"The probability of the union of those events is not as easy as it may seem. It requires application of the _inclusion exclusion_ principle repeatedly to obtain
```math
P(\large\cup_{i = 1}^n A_i) = \sum_{i = 1}^n P(A_i) - \sum_{\mbox{pairs}} P(A_i \cap A_j) + \sum_{\mbox{triplets}} P(A_i \cap A_j \cap A_k) + \ldots
```
"

# ╔═╡ e8c3415a-6240-11eb-328c-fd7e571ad7f9
md"So our probability becomes
```math
P(B) = 1 - \sum_{k = 1}^n (-1)^{k + 1} {n \choose k}p_k
```
where $p_k$ is the probability of at least $k$ people hearing about the news. From counting principles you can find out that
```math
p_k = \dfrac{(n - k)^r}{n^r} = \left(1 - \dfrac{k}{n}\right)^r
```
"

# ╔═╡ 74a09470-6241-11eb-00df-cb5032e755c3
function callprobabilty(n, r)
	return sum([(-1)^k * binomial(n, k) * (1 - k/n)^r for k in 0:n])
end

# ╔═╡ 9d189772-6241-11eb-16a3-25c54c8d4148
theorycurve = [callprobabilty(big(n), big(2*n)) for n in 1:100]

# ╔═╡ 024a8056-6242-11eb-23b6-a14e3fe68a74
plot(1:100, theorycurve)

# ╔═╡ 619109a4-6242-11eb-34eb-2527db9629bc
theorycurves = [ [callprobabilty(big(n), big(r_ratio*n)) for n in 1:200] for r_ratio in 1:3]

# ╔═╡ 85f02e2e-6242-11eb-2a27-7746c98beb47
plot(1:100, theorycurves, label=["r=n" "r=2n" "r=3n"])

# ╔═╡ ce1be85a-6242-11eb-026b-eb7215cc8c18
function simulatecalls(n, r, N)
	failed = float(0)
	for exp in 1:N
		counts = zeros(n)
		for k in 1:r
			person = rand(1:n)
			counts[person] += 1
		end
		if minimum(counts) == 0
			failed += 1
		end
	end
	return 1.0 - failed / N
end

# ╔═╡ 363c8308-6249-11eb-3580-8b721c00a424

# ╔═╡ e0efb4ea-6245-11eb-2423-b7d55a99bf0e
empiricalcurves = [[simulatecalls(n, r_ratio*n, 1000) for n in 1:100] for r_ratio in 1:3]

# ╔═╡ 1aaa0faa-6246-11eb-1f5c-69cbf277efea
begin
	scatter(1:100, empiricalcurves, label="")
	plot!(1:100, theorycurves, label=["r=n" "r=2n" "r=3n"])
end

# ╔═╡ 55b868f4-62bd-11eb-16a8-e14408e3e4bd
md"### Testing independence in random number generators"

# ╔═╡ a1c8ec8c-6a1f-11eb-0a85-fb5cac062c56
md"Randomness is a big requirement in statistical analysis and computational methods that depend on sampling."

# ╔═╡ 95bd11e8-6a1f-11eb-3287-0b32ba568406
md"👉 Do you know of any simple way to generate pseudo random numbers?"

# ╔═╡ abe1778a-6249-11eb-3b17-0de274e8e3fe
md"
A pseudo-random generator is given by the following sequence of outputs:
```math
X_{n + 1} = (a X_n + c) \mbox{ mod } m
```
"

# ╔═╡ 277b4f60-624a-11eb-2e11-9b6c5aeb3b4a
function prng(n)
	c = 7
	m = 13
	a = 5
	out = zeros(Int, n)
	out[1] = floor(rand() * m)
	for i in 2:n
		out[i] = (a * out[i - 1] + c) % m
	end
	return out
end

# ╔═╡ 9645b9a6-624a-11eb-09e9-c399cb136a5c
prng(5)

# ╔═╡ f6e2aff6-624c-11eb-16b5-318c5f230a56
n_vals = 1000

# ╔═╡ b97aa654-6a1f-11eb-11da-f91221057751
md"👉 Are two consecutive random numbers generated with `prng` independent? How would you check this computationally?"

# ╔═╡ ff30f85c-624b-11eb-3cf2-c170108d7092
vals1 = [prng(2) for i in 1:n_vals]

# ╔═╡ e59bc748-6a1f-11eb-2463-7dc3313f4b35
md"Remember that two random variables $X$ and $Y$ are indpendent if
```math
P(X, Y) = P(X) P(Y)
```"

# ╔═╡ b197bcf2-624c-11eb-3b99-efbedf0c6f3b
begin
	a = rand(1:13)
	b = rand(1:13)
	p_a = sum([vals1[i][1] == a for i in 1:n_vals]) / n_vals
	p_b = sum([vals1[i][2] == b for i in 1:n_vals]) / n_vals
	p_ap_b = p_a * p_b
	p_ab = sum([vals1[i][1] == a && vals1[i][2] == b for i in 1:n_vals]) / n_vals
	"a=$a b=$b p(x = a)=$p_a p(y = b)=$p_b p(x = a, y = b)=$p_ab p(x = a)p(y = b)=$p_ap_b"
end

# ╔═╡ 6b0245fe-624f-11eb-38f1-e702814d57c3
vals2 = [[rand(1:13), rand(1:13)] for i in 1:n_vals]

# ╔═╡ 5193bb62-62bc-11eb-1df7-1b7a7e9983ed
begin
	p_a2 = sum([vals2[i][1] == a for i in 1:n_vals]) / n_vals
	p_b2 = sum([vals2[i][2] == b for i in 1:n_vals]) / n_vals
	p_a2p_b2 = p_a2 * p_b2
	p_ab2 = sum([vals2[i][1] == a && vals2[i][2] == b for i in 1:n_vals]) / n_vals
	"a=$a b=$b p(x = a)=$p_a2 p(y = b)=$p_b2 p(x = a, y = b)=$p_ab2 p(x = a)p(y = b)=$p_a2p_b2"
end

# ╔═╡ 76fffb4e-62bd-11eb-153b-690eeaed1160
md"### Distributions and their plots"

# ╔═╡ fae8040a-62be-11eb-18a4-9fe8dee70bed
D_n = Normal(0, 1)

# ╔═╡ 0ea764ea-62bf-11eb-0e98-437744ee9fa9
plot(-5:0.1:5, [pdf(D_n, x) for x in -5:0.1:5], label="Normal PDF")

# ╔═╡ 3771b0ec-62bf-11eb-0300-6d9503c5716b
plot(-5:0.1:5, [cdf(D_n, x) for x in -5:0.1:5], label="CDF")

# ╔═╡ 66abd1d6-62c3-11eb-2e4e-1f67f89c3526
begin
	plot(-5:0.1:5, [cdf(D_n, x) for x in -5:0.1:5], label="CDF")
	plot!(-5:0.1:5, [quadgk(x -> pdf(D_n, x), (-5, y)...)[1] for y in -5:0.1:5], line = (4, :dash), label="From integration")
end

# ╔═╡ c420bab4-62c0-11eb-13c5-51ef6be81102
D_b = Binomial(5, 0.5)

# ╔═╡ ee11743a-62c0-11eb-04ef-8f30696c46ee
plot(0:5, [pdf(D_b, x) for x in 0:5], line=:stem, marker=:circle, label="Binomial")

# ╔═╡ 66a1bfe2-62bf-11eb-36d4-cfcaff920522
mean(D_n)

# ╔═╡ 93a82e9c-62cc-11eb-1b9c-cbc0e30b623e
md"### Moments of distributions"

# ╔═╡ f5b41e14-62bf-11eb-3257-6f5f5d6ceee6
md"What is the mean?
For a discrete distribution:
```math
\mathbb{E}[x]  = \sum_x xp(x)
```
For a continuous distribution:
```math
\mathbb{E}[x]  = \int_{-\infty}^{\infty} x f(x) dx,
```
where $f(x)$ is the pdf of the distribution.
"

# ╔═╡ 2be88b12-62c2-11eb-14de-6912598bc6bb
area = quadgk(x -> pdf(D_n, x), (-5, 5)...)

# ╔═╡ ec35489c-62bf-11eb-0311-2fed6b29d0c6
computed_mean = quadgk(x -> x * pdf(D_n, x), (-5, 5)...)

# ╔═╡ 5dffb66c-62c1-11eb-32f6-4b0b17d2aed1
mean(Binomial(5, 0.5))

# ╔═╡ 35406a36-62c2-11eb-347e-9b03c368119a
tot_length = sum([pdf(D_b, x) for x in 0:5])

# ╔═╡ 653515e4-62c1-11eb-25b0-930f782be996
computed_mean_binominal = sum([x * pdf(D_b, x) for x in 0:5])

# ╔═╡ 6a1d4ba4-62c2-11eb-3b71-431ee0c1fdf5
md"
```math
\mbox{Var}(X) = \mathbb{E}[X^2] - (\mathbb{E}[X])^2
```
"

# ╔═╡ 565657a8-62c2-11eb-22ad-9fbf7e7b44e0
var(D_n)

# ╔═╡ 641557fe-62c2-11eb-1ffd-0b1486231fa0
quadgk(x -> x^2 * pdf(D_n, x), (-5, 5)...)[1] - mean(D_n)^2

# ╔═╡ a46d7728-62cc-11eb-3486-d799a3d29298
md"### Sampling and sums of distributions"

# ╔═╡ 5811e2d8-62cc-11eb-266a-3d1ed4d9270b
n_samples = 1000

# ╔═╡ 2d0191b4-62c3-11eb-1479-21830ac879eb
x1 = rand(D_n, n_samples)

# ╔═╡ c845fae0-62c6-11eb-2a29-bdbc0c07d6d4
histogram(x1, normalize=true)

# ╔═╡ bd903520-62c6-11eb-1f71-49b662373de7
x2 = rand(D_b, n_samples)

# ╔═╡ 485386b2-62cc-11eb-1dd4-31857e128d7f
histogram(x2, normalize=true)

# ╔═╡ c4f908fa-62c6-11eb-259d-576dee85efb9
x = x1 + x2

# ╔═╡ 505c1658-62cc-11eb-1068-659d4286a4da
histogram(x, normalize=true)

# ╔═╡ 6e365fee-62cc-11eb-28d4-cd498c21895d
begin
	histogram(x, normalize=true)
	plot!(-5:0.1:5, [pdf(D_n, x)+pdf(D_b, x) for x in -5:0.1:5])
end

# ╔═╡ ad175b5a-62cc-11eb-2913-8741eca3d89c
md"The random variable which is the sum of two independent random variable does not necessarily have a pdf as the sum of the pdfs of the two random variables."

# ╔═╡ c19c5abc-62cc-11eb-3bcd-d9ee89e4619d
begin
	conv(x) = sum(pdf(D_n,x-k)*pdf(D_b,k) for k=-5:0.1:8)
	plot(-5:0.1:8, conv.(-5:0.1:8))
end

# ╔═╡ e6d9958c-62cd-11eb-164e-4b52aedd1584
md"The convolution operation on the individual pdfs creates the pdf of the sum of two random variables. It is defined as:
```math
h(z) = (p1 \otimes p2)(z) = \int_{-\infty}^{\infty} f(t) g(z - t) dt
```
This needs to be changed to discrete versions if distribution/s are discrete.
"

# ╔═╡ b68f9fdc-62cd-11eb-1852-0116fe475f36
begin
	histogram(x, normalize=true)
	plot!(-5:0.1:8, conv.(-5:0.1:8))
end

# ╔═╡ 7b02ff44-62ce-11eb-2c40-b5a734f869cb
md"### Inverse"

# ╔═╡ 7e5765f6-62ce-11eb-222e-9786c5fc0911
md"When we say that the $\mbox{CDF}(x) = u$, we mean that the probability of the random number being less than or equal to $x$ is $u$. We can now ask an opposite question: For what value of $x$ is the probability of the random number being less than $x$ is a given $u$. This is usually referred to as the _inverse function_."

# ╔═╡ beea6cda-62ce-11eb-0ef1-253f1a40504b
md"Though not mathematically tight, we can define the inverse function $\mbox{CDF}^{-1}$ as:
```math
\mbox{CDF}^{-1}(u) = \inf\{x\ :\ \mbox{CDF}(x) \geq u\}
```
"

# ╔═╡ 0047b980-62cf-11eb-3ca2-ad2f18f0e091
cdf_n(x) = cdf(Normal(0, 1), x)

# ╔═╡ 1676307e-62cf-11eb-39e5-bd732850ab2f
cdf_n_inv(u) = minimum(filter((x) -> (cdf_n(x) >= u), -5:0.01:5))

# ╔═╡ 343f9c62-62cf-11eb-2d09-cff6967a1d65
plot(-5:0.1:5, [cdf_n(x) for x in -5:0.1:5])

# ╔═╡ 5fab01a2-62cf-11eb-1c61-f9da64f9a75f
plot(0.1:0.01:0.9, [cdf_n_inv(u) for u in 0.1:0.01:0.9])

# ╔═╡ 9a083902-62d2-11eb-2570-bdb7ea1e1f5e
md"Let us relate the inverse function to area under the PDF curve"

# ╔═╡ c483fc1c-62d1-11eb-14b6-f50474a58fb2
cdf_n_inv(0.8)

# ╔═╡ a9a39152-62d2-11eb-3b4a-5d37f7195f3a
cdf_n_inv(0.95)

# ╔═╡ f437e702-62d1-11eb-1ac8-63c7f6ec16c4
xslider = @bind x_slider html"<input type=range min=-5 max=5 step=0.01>"

# ╔═╡ 8439e2ec-62d2-11eb-3772-d17e099debb5
x_slider

# ╔═╡ 794c86e6-62d2-11eb-3e27-ebd8a1b3d4e3
u_cdf = cdf_n(x_slider)

# ╔═╡ cb5e2858-62d1-11eb-1489-037703d411ec
begin
	plot(x->x, x->pdf(Normal(0, 1), x), -10, 10, label="")
	plot!(x->x, x->pdf(Normal(0, 1), x), x_slider, 10, fill=(0, :orange), label="")
end

# ╔═╡ 2a8682e2-62d3-11eb-119e-cbd804e1b0b2
md"### Fitting given data to distributions"

# ╔═╡ a089b976-62d3-11eb-2702-8945a47bf858
n_slider = @bind n_samples_exp html"<input type=range min=1 max=5 step=1>"

# ╔═╡ 3151c634-62d3-11eb-2698-95e4e14aab2d
samples = rand(Normal(10, 2.5), 10^n_samples_exp) + rand(Normal(0, 1), 10^n_samples_exp)

# ╔═╡ 62e55e4a-62d3-11eb-21e3-7b02354ddd02
histogram(samples)

# ╔═╡ 786f3d9c-62d3-11eb-28a7-b387242e3760
fit(Normal, samples)

# ╔═╡ 9531fffe-62d3-11eb-124c-0d2cccea92ae
md"Notice how the estimate of μ is reasonably unbiased even for small sizes of the samples, but how σ seems to be higher."

# ╔═╡ Cell order:
# ╟─ebd217ba-61f8-11eb-2337-9db63173ed34
# ╠═ebcb6ec6-61f9-11eb-05a3-817618350415
# ╠═2d1693e6-61fa-11eb-3cb6-9f5431449a41
# ╟─4e0bacac-621c-11eb-17be-a968500dc6bc
# ╟─48a0d2ac-6a1d-11eb-2441-897ca192e8ad
# ╠═29fb09a2-61f9-11eb-0699-4dbf7daca0d7
# ╠═bbae804a-61f9-11eb-03fe-89049682eae7
# ╟─e977219e-620d-11eb-2ca1-2141cc7c7d2a
# ╠═ee1bd118-61f9-11eb-1a18-df7d29a96c97
# ╠═059c87a6-61fa-11eb-01b4-09c40af1be94
# ╠═0a10f34e-61fa-11eb-365d-01ff16513530
# ╠═e00f4d02-61f9-11eb-2c5e-19ab03f4bdeb
# ╟─20ea39f4-620e-11eb-05b3-716b88bcd205
# ╠═c0b6a868-61fb-11eb-2190-83d6d760ede4
# ╠═947c880a-6203-11eb-1bfe-95c8e630b395
# ╠═fb898100-6204-11eb-1faf-cd75283f54d4
# ╠═b25c21be-6203-11eb-2f40-5db01b4beb89
# ╟─2b3a4f0c-620e-11eb-2021-bdaab9318fc1
# ╟─c2ae2510-620a-11eb-07e0-91cf24657a55
# ╠═853a1ab6-6a1d-11eb-2e08-93e5bcc73c7c
# ╟─33d3be8c-620e-11eb-3f7d-8f9a60cfdefc
# ╠═6205d7aa-620d-11eb-3c75-83b491a84995
# ╟─6f0eec76-6a1d-11eb-2585-17959367c09f
# ╟─ba50d33e-6a1d-11eb-0687-47cb0b6f8242
# ╠═ec8b9cd0-6211-11eb-06fa-37880ae43e40
# ╟─e37eefb6-6a1d-11eb-000b-89dbddabce5d
# ╟─5b2accf0-620e-11eb-06b9-ef354ebb3261
# ╠═7dced6ec-6205-11eb-3e95-5f2fa0d5adf3
# ╟─7bbe53a6-620e-11eb-296a-dbf3eaab558d
# ╠═014acff6-6206-11eb-0f3f-f17ef6edf8dd
# ╠═987d2bb2-620d-11eb-07b6-1736ac133029
# ╠═0c205d04-6205-11eb-2d31-5f689e25a3f1
# ╟─beccfb34-620e-11eb-177a-711b597f5f6f
# ╠═b0c7733e-620e-11eb-0fb6-f3a9f7631823
# ╟─ce084cfc-620e-11eb-2881-8f30a7558bba
# ╠═f5d3f88c-620e-11eb-3bc8-e3e0e1598482
# ╠═d8b7f274-620e-11eb-13b1-3933af0f7157
# ╟─77fb84bc-621c-11eb-2cc3-71425a74f823
# ╟─59437fd2-621c-11eb-1e31-f1973639a0e2
# ╟─8fbbc5a6-6212-11eb-1020-7382ef578cea
# ╟─beeacf18-6212-11eb-1929-a5c830dddbcb
# ╠═6aa6b4ec-6212-11eb-39ea-db9d36b22cce
# ╠═bc4822e0-6212-11eb-3a49-e5307fe49156
# ╟─d5215c20-6212-11eb-32f5-c71343b60d22
# ╟─e6f2b97e-6212-11eb-0ef4-6593e840b93a
# ╠═e49d22cc-6212-11eb-2732-27ad5e863ffc
# ╟─f462afb0-6212-11eb-1a80-239d798e25f8
# ╟─383c0f2e-6213-11eb-355e-5b16663f1897
# ╠═559fea86-6213-11eb-185e-b5fc917251fe
# ╟─73f6494e-6213-11eb-2a42-7b58fe47d597
# ╟─1ad550e8-6214-11eb-108b-0d6e08518018
# ╠═e455c886-6213-11eb-01b8-5df9e2bb7ea8
# ╟─33bc8ece-6a1f-11eb-20c9-fd32413cc2ac
# ╟─17ebe7ac-6214-11eb-38f1-3b757972e958
# ╠═3d281fc0-6214-11eb-16ae-89900d7f437f
# ╟─94d49374-621c-11eb-2985-b5ee5644c67e
# ╟─5eff3d3a-621c-11eb-0824-f76604aed501
# ╟─ba16f60e-623f-11eb-190a-af4758c32068
# ╟─51413a62-6a1f-11eb-1400-5918c088122c
# ╟─eed974c0-623f-11eb-021e-ebfc3c85b934
# ╟─2cf52d58-6240-11eb-1e21-f3c0a71b4525
# ╟─e8c3415a-6240-11eb-328c-fd7e571ad7f9
# ╠═74a09470-6241-11eb-00df-cb5032e755c3
# ╠═9d189772-6241-11eb-16a3-25c54c8d4148
# ╠═024a8056-6242-11eb-23b6-a14e3fe68a74
# ╠═619109a4-6242-11eb-34eb-2527db9629bc
# ╠═85f02e2e-6242-11eb-2a27-7746c98beb47
# ╠═ce1be85a-6242-11eb-026b-eb7215cc8c18
# ╠═363c8308-6249-11eb-3580-8b721c00a424
# ╠═e0efb4ea-6245-11eb-2423-b7d55a99bf0e
# ╠═1aaa0faa-6246-11eb-1f5c-69cbf277efea
# ╟─55b868f4-62bd-11eb-16a8-e14408e3e4bd
# ╟─a1c8ec8c-6a1f-11eb-0a85-fb5cac062c56
# ╟─95bd11e8-6a1f-11eb-3287-0b32ba568406
# ╟─abe1778a-6249-11eb-3b17-0de274e8e3fe
# ╠═277b4f60-624a-11eb-2e11-9b6c5aeb3b4a
# ╠═9645b9a6-624a-11eb-09e9-c399cb136a5c
# ╠═f6e2aff6-624c-11eb-16b5-318c5f230a56
# ╟─b97aa654-6a1f-11eb-11da-f91221057751
# ╠═ff30f85c-624b-11eb-3cf2-c170108d7092
# ╟─e59bc748-6a1f-11eb-2463-7dc3313f4b35
# ╠═b197bcf2-624c-11eb-3b99-efbedf0c6f3b
# ╠═6b0245fe-624f-11eb-38f1-e702814d57c3
# ╠═5193bb62-62bc-11eb-1df7-1b7a7e9983ed
# ╟─76fffb4e-62bd-11eb-153b-690eeaed1160
# ╠═6b69c792-62bd-11eb-0f1d-2d0081598dfa
# ╠═fae8040a-62be-11eb-18a4-9fe8dee70bed
# ╠═0ea764ea-62bf-11eb-0e98-437744ee9fa9
# ╠═3771b0ec-62bf-11eb-0300-6d9503c5716b
# ╠═66abd1d6-62c3-11eb-2e4e-1f67f89c3526
# ╠═c420bab4-62c0-11eb-13c5-51ef6be81102
# ╠═ee11743a-62c0-11eb-04ef-8f30696c46ee
# ╠═66a1bfe2-62bf-11eb-36d4-cfcaff920522
# ╠═b80b8e58-62bf-11eb-11fb-1397b06ae3e1
# ╟─93a82e9c-62cc-11eb-1b9c-cbc0e30b623e
# ╟─f5b41e14-62bf-11eb-3257-6f5f5d6ceee6
# ╠═2be88b12-62c2-11eb-14de-6912598bc6bb
# ╠═ec35489c-62bf-11eb-0311-2fed6b29d0c6
# ╠═5dffb66c-62c1-11eb-32f6-4b0b17d2aed1
# ╠═35406a36-62c2-11eb-347e-9b03c368119a
# ╠═653515e4-62c1-11eb-25b0-930f782be996
# ╟─6a1d4ba4-62c2-11eb-3b71-431ee0c1fdf5
# ╠═565657a8-62c2-11eb-22ad-9fbf7e7b44e0
# ╠═641557fe-62c2-11eb-1ffd-0b1486231fa0
# ╟─a46d7728-62cc-11eb-3486-d799a3d29298
# ╠═5811e2d8-62cc-11eb-266a-3d1ed4d9270b
# ╠═2d0191b4-62c3-11eb-1479-21830ac879eb
# ╠═c845fae0-62c6-11eb-2a29-bdbc0c07d6d4
# ╠═bd903520-62c6-11eb-1f71-49b662373de7
# ╠═485386b2-62cc-11eb-1dd4-31857e128d7f
# ╠═c4f908fa-62c6-11eb-259d-576dee85efb9
# ╠═505c1658-62cc-11eb-1068-659d4286a4da
# ╠═6e365fee-62cc-11eb-28d4-cd498c21895d
# ╟─ad175b5a-62cc-11eb-2913-8741eca3d89c
# ╠═c19c5abc-62cc-11eb-3bcd-d9ee89e4619d
# ╟─e6d9958c-62cd-11eb-164e-4b52aedd1584
# ╠═b68f9fdc-62cd-11eb-1852-0116fe475f36
# ╟─7b02ff44-62ce-11eb-2c40-b5a734f869cb
# ╟─7e5765f6-62ce-11eb-222e-9786c5fc0911
# ╟─beea6cda-62ce-11eb-0ef1-253f1a40504b
# ╠═0047b980-62cf-11eb-3ca2-ad2f18f0e091
# ╠═1676307e-62cf-11eb-39e5-bd732850ab2f
# ╠═343f9c62-62cf-11eb-2d09-cff6967a1d65
# ╠═5fab01a2-62cf-11eb-1c61-f9da64f9a75f
# ╟─9a083902-62d2-11eb-2570-bdb7ea1e1f5e
# ╠═c483fc1c-62d1-11eb-14b6-f50474a58fb2
# ╠═a9a39152-62d2-11eb-3b4a-5d37f7195f3a
# ╠═f437e702-62d1-11eb-1ac8-63c7f6ec16c4
# ╠═8439e2ec-62d2-11eb-3772-d17e099debb5
# ╠═794c86e6-62d2-11eb-3e27-ebd8a1b3d4e3
# ╠═cb5e2858-62d1-11eb-1489-037703d411ec
# ╟─2a8682e2-62d3-11eb-119e-cbd804e1b0b2
# ╠═a089b976-62d3-11eb-2702-8945a47bf858
# ╠═3151c634-62d3-11eb-2698-95e4e14aab2d
# ╠═62e55e4a-62d3-11eb-21e3-7b02354ddd02
# ╠═786f3d9c-62d3-11eb-28a7-b387242e3760
# ╟─9531fffe-62d3-11eb-124c-0d2cccea92ae

max_stars_repo_name: divkakwani/pluto-on-binder


text:
using Test
using DFTK

include("./testcases.jl")

@testset "Diagonalization of a free-electron Hamiltonian" begin
    # Construct a free-electron Hamiltonian
    Ecut = 5
    fft_size = [15, 15, 15]
    model = Model(silicon.lattice, n_electrons=silicon.n_electrons, terms=[Kinetic()])  # free-electron model
    basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)
    ham = Hamiltonian(basis)

    tol = 1e-8
    nev_per_k = 10
    ref_λ = [
        [0.00000000000, 0.56219939834, 0.56219939834, 0.56219939834, 0.56219939834,
         0.56219939834, 0.56219939834, 0.56219939834, 0.56219939834, 0.74959919778],
        [0.06246659981, 0.24986639926, 0.49973279852, 0.49973279852, 0.49973279852,
         0.56219939834, 0.56219939834, 0.56219939834, 0.74959919778, 0.74959919778],
        [0.08328879975, 0.33315519901, 0.39562179883, 0.39562179883, 0.39562179883,
         0.39562179883, 0.83288799753, 0.83288799754, 0.83288799754, 0.83288799754],
        [0.16657759951, 0.22904419932, 0.22904419932, 0.41644399877, 0.41644399877,
         0.66631039803, 0.72877699784, 0.72877699784, 0.72877699784, 0.72877699784],
    ]

    @test length(ref_λ) == length(silicon.kcoords)
    @testset "without Preconditioner" begin
        res = diagonalize_all_kblocks(lobpcg_hyper, ham, nev_per_k, tol=tol,
                                      prec_type=nothing, interpolate_kpoints=false)

        @test res.converged
        for ik in 1:length(silicon.kcoords)
            @test ref_λ[ik] ≈ res.λ[ik]
            @test maximum(res.residual_norms[ik]) < 100tol  # TODO Why the 100?
            @test res.iterations[ik] < 200
        end
    end

    @testset "with Preconditioner" begin
        res = diagonalize_all_kblocks(lobpcg_hyper, ham, nev_per_k, tol=tol,
                                      prec_type=PreconditionerTPA, interpolate_kpoints=false)

        @test res.converged
        for ik in 1:length(silicon.kcoords)
            @test ref_λ[ik] ≈ res.λ[ik]
            @test maximum(res.residual_norms[ik]) < 100tol  # TODO Why the 100?
            @test res.iterations[ik] < 50
        end
    end
end

if !isdefined(Main, :FAST_TESTS) || !FAST_TESTS
    @testset "Diagonalization of kinetic + local PSP" begin
        Ecut = 25
        fft_size = [33, 33, 33]

        Si = ElementPsp(silicon.atnum, psp=load_psp("hgh/lda/si-q4"))
        model = Model(silicon.lattice; atoms=[Si => silicon.positions],
                      terms=[Kinetic(),AtomicLocal()])
        basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops;
                               fft_size=fft_size)
        ham = Hamiltonian(basis)

        res = diagonalize_all_kblocks(lobpcg_hyper, ham, 6, tol=1e-8)

        ref = [
            [-4.087198659513310, -4.085326314828677, -0.506869382308294,
             -0.506869382280876, -0.506869381798614],
            [-4.085824585443292, -4.085418874576503, -0.509716820984169,
             -0.509716820267449, -0.508545832298541],
            [-4.086645155119840, -4.085209948598607, -0.514320642233337,
             -0.514320641863231, -0.499373272772206],
            [-4.085991608422304, -4.085039856878318, -0.517299903754010,
             -0.513805498246478, -0.497036479690380]
        ]
        for ik in 1:length(silicon.kcoords)
            @test res.λ[ik][1:5] ≈ ref[ik] atol=5e-7
        end
    end
end

@testset "Diagonalization of a core Hamiltonian" begin
    Ecut = 10
    fft_size = [21, 21, 21]

    Si = ElementPsp(silicon.atnum, psp=load_psp("hgh/lda/si-q4"))
    model = Model(silicon.lattice; atoms=[Si => silicon.positions],
                  terms=[Kinetic(), AtomicLocal(), AtomicNonlocal()])

    basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)
    ham = Hamiltonian(basis)

    res = diagonalize_all_kblocks(lobpcg_hyper, ham, 5, tol=1e-8, interpolate_kpoints=false)
    ref = [
        [0.067955741977536, 0.470244204908046, 0.470244204920801,
         0.470244204998022, 0.578392222232969],
        [0.111089041747288, 0.304724122513625, 0.445322298067717,
         0.445322298101198, 0.584713217756577],
        [0.129419322499919, 0.293174377882115, 0.411932220567084,
         0.411932220611853, 0.594921264868345],
        [0.168662148987539, 0.238552367551507, 0.370743978236562,
         0.418387442903058, 0.619797227001203],
    ]
    for ik in 1:length(silicon.kcoords)
        @test res.λ[ik] ≈ ref[ik] atol=0.02
    end
end

@testset "Full diagonalization of a LDA Hamiltonian" begin
    Ecut = 2

    Si = ElementPsp(silicon.atnum, psp=load_psp("hgh/lda/si-q4"))
    model = model_DFT(silicon.lattice, [Si => silicon.positions], :lda_xc_teter93)
    basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops)
    ham = Hamiltonian(basis; ρ=guess_density(basis))

    res1 = diagonalize_all_kblocks(lobpcg_hyper, ham, 5, tol=1e-8, interpolate_kpoints=false)
    res2 = diagonalize_all_kblocks(diag_full, ham, 5, tol=1e-8, interpolate_kpoints=false)
    for ik in 1:length(silicon.kcoords)
        @test res1.λ[ik] ≈ res2.λ[ik] atol=1e-6
    end
end


text:
if length(ARGS) < 1
    println("Please provide file -- fasta format")
    exit(0)
else
    fastaf = ARGS[1]
end

function parse_F(fa_file)
    f_h = open(fa_file, "r")
    n = countlines(f_h)
    seekstart(f_h)

    seqs_Dict = Dict{String, String}()
    id_array = Array{String, 1}()

    for i = 1:n
        Line = readline(f_h)

        if (Line[1] == '>')
            Name= rstrip(Line)[2:end]
            seqs_Dict[Name] = ""
            push!(id_array, Name)

        else

            t_name = id_array[end]

            try
                seqs_Dict[t_name] *= rstrip(Line)
            catch UndefVarError       #catch KeyError
                seqs_Dict[t_name] = rstrip(Line)

            #throw OverflowError()
            #    continue
            end
        end
    end
    return seqs_Dict, id_array
end

seq_d, seq_names = parse_F(fastaf)

for (key, value) in seq_d
    println(key," => ", value)
end

for name in seq_names
    println(name)
end


text:
using JLD, GibbsTypePriors, DataFramesMeta, RCall, Optim, OptimTestProblems, StatsBase

include("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/scripts_for_figures/common_functions.jl")

R"
StickBreakingPY <- function(alpha,sigma, N) {
  a<- rep(1-sigma,N-1)
  b<- alpha+ sigma*(1:N-1)
  V <- rbeta(N-1 , a, b)
  p    <- vector('numeric',length=N)
  p[1] <- V[1]
  for(l in 2:(N - 1))p[l] <- prod(1 - V[1:(l - 1)])*V[l]
  p[N] <- prod(1 - V)
  p
  return(p)
}

Prior_on_K_SB_PY<- function(alpha, sigma, N_s=100,N_tr=100, runs=10^4 ){
  array_nc<-c()
  i=1
  while (i<=runs){
    weights_PY<- StickBreakingPY(alpha,sigma,N_tr)
    c <- sample(1:N_tr,size=N_s, replace=TRUE, prob=weights_PY )
    n_c<- length(unique(c))
    array_nc[i]<- n_c
    i=i+1
  }
  p_k = tibble(k=as.numeric(names(table(array_nc))),
               p_k=as.vector(table(array_nc))/sum(table(array_nc)))
  p_zeros= tibble(k=(1:N_s)[!1:N_s%in%p_k$k],
                  p_k=rep(0,length((1:N_s)[!1:N_s%in%p_k$k])))
  pk_padded= rbind(p_k, p_zeros)%>% arrange(k)
  E_k=sum((1:N_s) *(pk_padded$p_k))
  V_k = sum((((1:N_s) - E_k)^2) *(pk_padded$p_k))
  return(list(pk=pk_padded$p_k, Ek= E_k,Vk=V_k ))
}
"

function EV_of_number_of_clusters_PY(n,par_vec)
    #pk_py = Pkn_NGG.(1:n, n, par_vec[1], par_vec[2])
    E =  GibbsTypePriors.E_PY(n, par_vec[1], par_vec[2]) |> Float64
    V  = GibbsTypePriors.V_PY(n, par_vec[1], par_vec[2])  |> Float64
    return E,V
end

sigma = collect(range(0.05,0.99, length=10))
alpha = collect(exp.(range(log(1), log(200), length =10)))
grid = collect(Iterators.product(alpha, sigma))
n=100
grid_vec = vec(grid)
H= 250

save("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_7_5.jld", "grid_vec" ,grid_vec)

## PY exact

EV_PY = EV_of_number_of_clusters_PY.(n, grid_vec)

DF_PY= DataFrame(Exp = first.(EV_PY),Var =last.(EV_PY), beta = first.(grid_vec), sigma = last.(grid_vec))
DF_PY.std = sqrt.(DF_PY.Var)
DF_PY.beta_log = log.(DF_PY.beta)

DF_PY

save("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_7_1.jld", "DF_PY" ,DF_PY)

## PY multinomial process

function EV_of_number_of_clusters_PYM(n,H,par_vec)
    pk_pym = GibbsTypePriors.Pkn_PYM.(1:n, n, H,  par_vec[1], par_vec[2])
    E_pym =  pk_pym|> ar -> map(*, ar, 1:n) |> sum
    x = ((1:n).-E_pym).^2
    V_pym = pk_pym |> ar -> map(*, ar, x) |> sum
    return E_pym, V_pym
end

EV_PYM = EV_of_number_of_clusters_PYM.(n,H,grid_vec)

DF_PYM= DataFrame(Exp = first.(EV_PYM),Var =last.(EV_PYM), beta = first.(grid_vec), sigma = last.(grid_vec))
DF_PYM.std = sqrt.(DF_PYM.Var)
DF_PYM.beta_log = log.(DF_PYM.beta)

DF_PYM

save("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_7_2.jld", "DF_PYM" ,DF_PYM)

## PY predictive approximation

function EV_of_number_of_clusters_PY_approx(n,par_vec)
    pk_approx = GibbsTypePriors.Pkn_PY_pred_approx(n, par_vec[1], par_vec[2])
    E_approx =  pk_approx|> ar -> map(*, ar, 1:n) |> sum
    x = ((1:n).-E_approx).^2
    V_approx = pk_approx |> ar -> map(*, ar, x) |> sum
    return E_approx, V_approx
end

EV_PY_apr= EV_of_number_of_clusters_PY_approx.(n, grid_vec)

DF_PY_apr= DataFrame(Exp = first.(EV_PY_apr),Var =last.(EV_PY_apr), beta = first.(grid_vec), sigma = last.(grid_vec))
DF_PY_apr.std = sqrt.(DF_PY_apr.Var)
DF_PY_apr.beta_log = log.(DF_PY_apr.beta)

DF_PY_apr

save("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_7_3.jld", "DF_PY_apr" ,DF_PY_apr)

##PY stick-breaking approximation

function EV_of_number_of_clusters_PY_SB(n,H, par_vec)
    R"SB_PY = Prior_on_K_SB_PY($(par_vec[1]),$(par_vec[2]), $n, $H, runs= 2*10^2)
    E_py_sb = SB_PY$Ek
    V_py_sb = SB_PY$Vk"
    @rget E_py_sb
    @rget V_py_sb
    return E_py_sb, V_py_sb
end

EV_PY_sb= EV_of_number_of_clusters_PY_SB.(n,H,  grid_vec)

DF_PY_sb= DataFrame(Exp = first.(EV_PY_sb),Var =last.(EV_PY_sb), beta = first.(grid_vec), sigma = last.(grid_vec))
DF_PY_sb.std = sqrt.(DF_PY_sb.Var)
DF_PY_sb.beta_log = log.(DF_PY_sb.beta)

DF_PY_sb

save("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_7_4.jld", "DF_PY_sb" ,DF_PY_sb)

max_stars_repo_name: konkam/approximatingBNPpriors

max_issues_repo_name: konkam/approximatingBNPpriors

max_forks_repo_name: konkam/approximatingBNPpriors


text:
struct Point{T}
    x::T
    y::T
end

@show Point{Float64} <: Point
@show Float64 <: Real
@show Point{Float64} <: Point{Real}
@show supertype(Point)

function norm(p::Point{Real})
    sqrt(p.x^2+p.y^2)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using HDF5
using Merlin

include("parser.jl")
include("model.jl")

const wordembeds_file = ".data/glove.6B.100d.h5"

#traindoc = CoNLL.read(".data/wsj_00-18.conll")
#testdoc = CoNLL.read(".data/wsj_22-24.conll")
#info("# sentences of train doc: $(length(traindoc))")
#info("# sentences of test doc: $(length(testdoc))")

#traindata = setup_data(traindoc, worddict, chardict)
#testdata = setup_data(testdoc, worddict, chardict)
#info("# words: $(length(worddict))")
#info("# chars: $(length(chardict))")
#traindata, testdata, worddict, chardict

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
@testset "auto_detect_ci_service" begin
    @testset "GitHub: env var present" begin
        withenv("GITHUB_REPOSITORY" => "true") do
            expected = CompatHelper.GitHubActions()

            @test CompatHelper.auto_detect_ci_service() == expected
        end
    end

    @testset "GitLab: env var present" begin
        withenv("GITLAB_CI" => "true", "GITHUB_REPOSITORY" => "false") do
            delete!(ENV, "GITHUB_REPOSITORY")
            expected = CompatHelper.GitLabCI()

            @test CompatHelper.auto_detect_ci_service() == expected
        end
    end

    @testset "env var dne" begin
        withenv("GITHUB_REPOSITORY" => "foobar", "GITLAB_CI" => "foobar") do
            delete!(ENV, "GITHUB_REPOSITORY")
            delete!(ENV, "GITLAB_CI")
            @test_throws CompatHelper.UnableToDetectCIService CompatHelper.auto_detect_ci_service()
        end
    end
end

@testset "$(func)" for func in [CompatHelper.ci_repository, CompatHelper.ci_token]
    value = "value"
    ci = CompatHelper.GitHubActions()

    @testset "exists" begin
        withenv("GITHUB_REPOSITORY" => value, "GITHUB_TOKEN" => value) do
            @test func(ci) == value
        end
    end

    @testset "dne" begin
        withenv("GITHUB_REPOSITORY" => value, "GITHUB_TOKEN" => value) do
            delete!(ENV, "GITHUB_REPOSITORY")
            delete!(ENV, "GITHUB_TOKEN")

            @test_throws KeyError func(ci)
        end
    end
end

@testset "$(func)" for func in [CompatHelper.ci_repository, CompatHelper.ci_token]
    value = "value"
    ci = CompatHelper.GitLabCI()

    @testset "exists" begin
        withenv("CI_PROJECT_PATH" => value, "GITLAB_TOKEN" => value) do
            @test func(ci) == value
        end
    end

    @testset "dne" begin
        withenv("CI_PROJECT_PATH" => value, "GITLAB_TOKEN" => value) do
            delete!(ENV, "CI_PROJECT_PATH")
            delete!(ENV, "GITLAB_TOKEN")

            @test_throws KeyError func(ci)
        end
    end
end


text:
const SymbolCollection = Union{Symbol, Tuple{Symbol, Vararg{Symbol}}, AbstractVector{Symbol}}

function _group_var(vc::VarCorr, group::Symbol)
    return sum(abs2, getproperty(getproperty(vc.σρ, group), :σ))
end

function _group_var(vc::VarCorr, groups::SymbolCollection)
    return sum(_group_var.(Ref(vc), groups))
end

function _group_var(vc::VarCorr)
    groups = collect(propertynames(vc.σρ))
    return _group_var(vc, groups)
end

"""
    icc(model::MixedModel, [groups])

Compute the intra-class correlation coefficient (ICC) for a mixed model.

The ICC is defined as the variance attributable to the `groups` divided by
the total variance from all groups and the observation-level (residual) variance.
In other words, the ICC can be interpreted as the proportion of the variance explainable
by the grouping/nesting structure.

A single `group` can be specified as a Symbol, e.g. `:subj` or a number of groups
can be specified as an array: `[:subj, :item]`. If no `groups` are specified, then
all grouping variables are used.

!!! note
    The value returned here is sometimes called the "adjusted ICC" and does not take
    the variance of the fixed effects into account (the "conditional ICC").

!!! note
    The result returned aggregates across groups. If you require the ICC for
    each group separately, then you must call `icc` separately for each group.
"""
function icc(model::GeneralizedLinearMixedModel,
             groups::Union{Symbol, SymbolCollection})
    dispersion_parameter(model) && throw(ArgumentError("GLMMs with dispersion parameters are not currently supported."))

    if model.resp.d isa Union{Binomial, Bernoulli}
        σ²res = π^2 / 3
    elseif model.resp.d isa Poisson
        σ²res = 1.0
    else
        throw(ArgumentError("Family $(typeof(model.resp.d)) currently unsupported, please file an issue."))
    end

    return _icc(VarCorr(model), groups, σ²res)
end

icc(model::LinearMixedModel,
    groups::Union{Symbol, SymbolCollection}) = _icc(VarCorr(model), groups, varest(model))

icc(model::MixedModel) = icc(model, fnames(model))

function _icc(vc::VarCorr, groups::Union{Symbol, SymbolCollection}, σ²res)
    σ²_α  = _group_var(vc, groups) # random effect(s) in numerator
    σ² = σ²res + _group_var(vc)
    return σ²_α / σ²
end


text:
using Pkg
Pkg.add("Plots")
Pkg.add("DifferentialEquations")
Pkg.add("StaticArrays")
Pkg.add("Parameters")


text:
# simple discrete hmm with testable output
A = eye(2)
B_matrix = [0.0  1.0  0.0 ;
            1.0  0.0  0.0 ]

model = HMM(A,B_matrix)
model.p = vec([1 0])

# state sequence s, observation sequence o
n_obs = 10
s,o = rand(model,n_obs)
@test(length(s) == length(o) == n_obs)
@test(all(s .== 1))
@test(all(o .== 2))

# check forward-backward algorithm on perfect sequence
alpha,beta,log_p_obs = forward_backward(model,o; scaling=false)
@test(log_p_obs == 0.0)
@test(all(alpha[:,1] .== 1))
@test(all(alpha[:,2] .== 0))
@test(all(beta[end,:] .== 1))
@test(all(beta[:,1] .== 1))

# test forward-backward with scaling
alpha,beta,log_p_obs = forward_backward(model,o)
@test(log_p_obs == 0.0)
@test(all(alpha[:,1] .== 1))
@test(all(alpha[:,2] .== 0))
@test(all(beta[end,:] .== 1))
@test(all(beta[:,1] .== 1))

# check viterbi algorithm on perfect sequence
@test(all(viterbi(model,o) .== s))

# check forward-backward algorithm on impossible sequence
alpha,beta,log_p_obs = forward_backward(model,o+1; scaling=false)
@test(log_p_obs == -Inf)
@test(all(alpha .== 0))
@test(all(beta[1:end-1,:] .== 0))

######
# Compare output to Michael Hamilton's "dishonest casino" example
# see Hamilton's HMM implementation in python:
# http://www.cs.colostate.edu/~hamiltom/code.html#python-hidden-markov-model

A = [ 0.99 0.01 ;
      0.01 0.99 ]
B = [ 1/6   1/6   1/6   1/6   1/6   1/6 ;
      1/10  1/10  1/10  1/10  1/10  1/2 ]
p = [0.5, 0.5]
model = HMM(A,B,p)

@test(model.n == 2)
for i = 1:model.n
  @test(typeof(model.B[i]) == Categorical)
  @test(model.B[i].K == 6)
end
@test(length(model.p) == model.n)

o = [1, 2, 1, 6, 6] # observation sequence

# check forward/backward (without scaling)
true_alpha = [ 8.3333e-02   1.3833e-02   2.2909e-03   3.7885e-04   6.295e-05  ;
               5.0000e-02   5.0333e-03   5.1213e-04   2.6496e-04   1.3305e-04 ]

true_beta = [ 8.9630e-04   5.2841e-03   3.0533e-02   1.7000e-01   1.0000e+00 ;
              2.4262e-03   2.4418e-02   2.4613e-01   4.9667e-01   1.0000e+00 ]

alpha,beta,log_p_obs = forward_backward(model, o; scaling=false)

@test(all(round(alpha,6) .== round(true_alpha',6)))
@test(all(round(beta,5) .== round(true_beta',5)))

# check forward/backward (with scaling)
alpha_scaled = [ 0.625   0.7332  0.8173  0.5884  0.3212 ;
                 0.375   0.2668  0.1827  0.4116  0.6788 ]

beta_scaled = [  4.5729   3.5946   2.9391   2.4312  3.2847 ;
                 12.3785  16.6108  23.6922  7.1028  3.2847 ]

alpha,beta,log_p_obs = forward_backward(model,o; scaling=true)

@test(all(round(alpha,4) .== round(alpha_scaled',4)))
@test(all(round(beta,4) .== round(beta_scaled',4)))

# check viterbi
true_vit = [2,2,2,2,2]
@test(all(viterbi(model,o) .== true_vit))
@test(all(viterbi(model,[1,2,3,4,5]) .== ones(Int,5)))

# another check for viterbi
model.A = [ 0.7 0.3 ;
          0.3 0.7 ]
o = [ 1, 2, 1, 6, 6, 3, 3, 2, 1 ]
true_vit = [ 1, 1, 1, 2, 2, 1, 1, 1, 1 ]
@test(all(viterbi(model,o) .== true_vit))

# check Baum-Welch algorithm
o = [ 1, 2, 1, 6, 6, 3, 3, 2, 1 ]
expected_result_A = [ 0.75  0.25 ;
                      0.25  0.75 ]
expected_result_B = [ 0.6  0.4  0.0  0.0  0.0  0.0 ;
                      0.0  0.0  0.5  0.0  0.0  0.5 ]
expected_result_p = [1.0, 0.0]

model1 = deepcopy(model) # check without scaling alpha/beta
fit!(model1, o; max_iter=60, scaling=false)

@test(all(round(model1.A,4) .== expected_result_A))
@test(all(round(model1.B[1].p,4) .== expected_result_B[1,:]'))
@test(all(round(model1.B[2].p,4) .== expected_result_B[2,:]'))
@test(all(round(model1.p,4) .== expected_result_p))

model2 = deepcopy(model) # check with scaling alpha/beta
fit!(model2, o; max_iter=60, scaling=true)
@test(all(round(model2.A,4) .== expected_result_A))
@test(all(round(model2.B[1].p,4) .== expected_result_B[1,:]'))
@test(all(round(model2.B[2].p,4) .== expected_result_B[2,:]'))
@test(all(round(model2.p,4) .== expected_result_p))


text: q="print(\"q=\$(repr(q)); eval(Meta.parse(q))\")"; eval(Meta.parse(q))

max_stars_repo_name: calvang/various-quines

max_issues_repo_name: calvang/various-quines

max_forks_repo_name: calvang/various-quines


text:
###### This file is part of the MendelImpute.jl package.
###### It contains relevant code to impute and output a genotype matrix

"""
    write(outfile, X, compressed_haplotypes, X_sampleID, XtoH_idx)

Writes imputed `X` into `outfile`. All genotypes in `outfile` are non-missing.
If `X` is a matrix, output VCF will be unphased. Otherwise if `X` is a Tuple
of matrix (i.e. `X = X1 + X2`), then all output will be phased.

# Notes
Here the writing routine is emulating `write_dlm` in Base at
https://github.com/JuliaLang/julia/blob/3608c84e6093594fe86923339fc315231492484c/stdlib/DelimitedFiles/src/DelimitedFiles.jl#L736
"""
function Base.write(
    outfile::AbstractString,
    X::Union{AbstractMatrix, Tuple{AbstractMatrix, AbstractMatrix}},
    compressed_haplotypes::CompressedHaplotypes,
    X_sampleID::AbstractVector,
    XtoH_idx::Union{Nothing, AbstractVector} = nothing,
    )
    threads = Threads.nthreads()
    snps = typeof(X) <: AbstractMatrix ? size(X, 1) : size(X[1], 1)
    len = div(snps, threads)
    files = ["tmp$i.vcf.gz" for i in 1:threads]

    # retrieve reference file information
    chr = (isnothing(XtoH_idx) ? compressed_haplotypes.chr :
                                 compressed_haplotypes.chr[XtoH_idx])
    pos = (isnothing(XtoH_idx) ? compressed_haplotypes.pos :
                                 compressed_haplotypes.pos[XtoH_idx])
    ids = (isnothing(XtoH_idx) ? compressed_haplotypes.SNPid :
                                 compressed_haplotypes.SNPid[XtoH_idx])
    ref = (isnothing(XtoH_idx) ? compressed_haplotypes.refallele :
                                 compressed_haplotypes.refallele[XtoH_idx])
    alt = (isnothing(XtoH_idx) ? compressed_haplotypes.altallele :
                                 compressed_haplotypes.altallele[XtoH_idx])

    # write minimal meta information to outfile
    io = [openvcf(files[i], "w") for i in 1:threads]
    pb = [PipeBuffer() for _ in 1:threads]
    print(pb[1], "##fileformat=VCFv4.2\n")
    print(pb[1], "##source=MendelImpute\n")
    print(pb[1], "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n")

    # header line should match reffile (i.e. sample ID's should match)
    print(pb[1], "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT")
    for id in X_sampleID
        print(pb[1], "\t", id)
    end
    print(pb[1], "\n")
    bytesavailable(pb[1]) > 1048576 && write(io[1], take!(pb[1]))

    # each thread writes `len` SNPs
    pmeter = Progress(snps, 5, "Writing to file...")
    Threads.@threads for t in 1:threads
        id = Threads.threadid()
        cur_ranges = (id == threads ?
            ((threads-1)*len+1:snps) : (1:len) .+ (t-1)*len)

        @inbounds for i in cur_ranges
            # write meta info (chrom/pos/snpid/ref/alt)
            print(pb[id], chr[i], "\t", string(pos[i]), "\t", ids[i][1], "\t",
                ref[i], "\t", alt[i][1], "\t.\tPASS\t.\tGT")
            # print ith record
            write_snp!(pb[id], X, i)
            bytesavailable(pb[id]) > 1048576 && write(io[id], take!(pb[id]))
            next!(pmeter)
        end
        write(io[id], take!(pb[id]))
    end
    close.(io); close.(pb) # close io and buffer

    # concatenate all files into 1 VCF file
    run(pipeline(`cat $files`, stdout=outfile))

    # delete intermediate files
    for i in 1:threads
        rm("tmp$i.vcf.gz", force=true)
    end

    return nothing
end

"""
Helper function for saving a record (SNP), not tracking phase information.
"""
function write_snp!(pb::IOBuffer, X::AbstractMatrix, i::Int)
    x = @view(X[i, :]) # current record
    n = length(x)
    @inbounds for j in 1:n
        if x[j] == 0
            print(pb, "\t0/0")
        elseif x[j] == 1
            print(pb, "\t1/0")
        elseif x[j] == 2
            print(pb, "\t1/1")
        else
            error("imputed genotypes can only be 0, 1, 2 but got $(x[j])")
        end
    end
    print(pb, "\n")
    nothing
end

"""
Helper function for saving a record (SNP), tracking phase information.
Here `X = X1 + X2`.
"""
function write_snp!(pb::IOBuffer, X::Tuple{AbstractMatrix, AbstractMatrix}, i::Int)
    X1, X2 = X[1], X[2]
    x1 = @view(X1[i, :])
    x2 = @view(X2[i, :])

    n = length(x1)
    @assert n == length(x2)
    @inbounds for j in 1:n
        if x1[j] == x2[j] == 0
            print(pb, "\t0|0")
        elseif x1[j] == 0 && x2[j] == 1
            print(pb, "\t0|1")
        elseif x1[j] == 1 && x2[j] == 0
            print(pb, "\t1|0")
        elseif x1[j] == 1 && x2[j] == 1
            print(pb, "\t1|1")
        else
            error("phased genotypes can only be 0|0, 0|1, 1|0 or 1|1 but
                got $(x1[j])|$(x2[j])")
        end
    end
    print(pb, "\n")
    nothing
end

"""
    impute!(X1, X2, H, phase, [impute_untyped])

Imputes `X = X1 + X2` completely using haplotype segments of `H`, where segments
information are stored in `phase`. `X1` is strand1 and `X2` is strand 2. If
`impute_untyped=false`, untyped SNPs will not be imputed.

Non-missing entries in `X` can be different after imputation.
"""
function impute!(
    X1::AbstractMatrix,
    X2::AbstractMatrix,
    compressed_Hunique::CompressedHaplotypes,
    phase::Vector{HaplotypeMosaicPair},
    impute_untyped::Bool=true
    )

    fill!(X1, 0)
    fill!(X2, 0)

    # loop over individuals
    for i in 1:size(X1, 2)
        # strand 1
        for s in 1:(length(phase[i].strand1.start) - 1) # first w-1 windows
            X_idx = phase[i].strand1.start[s]:(phase[i].strand1.start[s + 1] - 1)
            w = phase[i].strand1.window[s]
            H = impute_untyped ? compressed_Hunique.CW[w].uniqueH :
                compressed_Hunique.CW_typed[w].uniqueH
            H_start = impute_untyped ? (abs(phase[i].strand1.start[s] -
                compressed_Hunique.Hstart[w]) + 1) : 1
            H_idx = H_start:(H_start + length(X_idx) - 1)
            X1[X_idx, i] = H[H_idx, phase[i].strand1.haplotypelabel[s]]
        end
        w = phase[i].strand1.window[end] # last window
        X_end = impute_untyped ? phase[i].strand1.length : size(X1, 1)
        X_idx = phase[i].strand1.start[end]:X_end
        H_start = impute_untyped ? (abs(phase[i].strand1.start[end] -
            compressed_Hunique.Hstart[w]) + 1) : 1
        H_idx = H_start:(H_start + length(X_idx) - 1)
        H = impute_untyped ? compressed_Hunique.CW[w].uniqueH :
            compressed_Hunique.CW_typed[w].uniqueH
        X1[X_idx, i] = H[H_idx, phase[i].strand1.haplotypelabel[end]]

        # strand 2
        for s in 1:(length(phase[i].strand2.start) - 1) # first w-1 windows
            X_idx = phase[i].strand2.start[s]:(phase[i].strand2.start[s + 1] - 1)
            w = phase[i].strand2.window[s]
            H = impute_untyped ? compressed_Hunique.CW[w].uniqueH :
                compressed_Hunique.CW_typed[w].uniqueH
            H_start = impute_untyped ? (abs(phase[i].strand2.start[s] -
                compressed_Hunique.Hstart[w]) + 1) : 1
            H_idx = H_start:(H_start + length(X_idx) - 1)
            X2[X_idx, i] = H[H_idx, phase[i].strand2.haplotypelabel[s]]
        end
        w = phase[i].strand2.window[end] # last window
        X_end = impute_untyped ? phase[i].strand1.length : size(X2, 1)
        X_idx = phase[i].strand2.start[end]:X_end
        H = impute_untyped ? compressed_Hunique.CW[w].uniqueH :
            compressed_Hunique.CW_typed[w].uniqueH
        H_start = impute_untyped ? (abs(phase[i].strand2.start[end] -
            compressed_Hunique.Hstart[w]) + 1) : 1
        H_idx = H_start:(H_start + length(X_idx) - 1)
        X2[X_idx, i] = H[H_idx, phase[i].strand2.haplotypelabel[end]]
    end
end

function impute!(
    X1::AbstractMatrix,
    X2::AbstractMatrix,
    H::AbstractMatrix,
    phase::Vector{HaplotypeMosaicPair},
    )

    fill!(X1, 0)
    fill!(X2, 0)

    # loop over individuals
    for i in 1:size(X1, 2)
        for s in 1:(length(phase[i].strand1.start) - 1)
            idx = phase[i].strand1.start[s]:(phase[i].strand1.start[s + 1] - 1)
            X1[idx, i] = H[idx, phase[i].strand1.haplotypelabel[s]]
        end
        idx = phase[i].strand1.start[end]:phase[i].strand1.length
        X1[idx, i] = H[idx, phase[i].strand1.haplotypelabel[end]]
        for s in 1:(length(phase[i].strand2.start) - 1)
            idx = phase[i].strand2.start[s]:(phase[i].strand2.start[s + 1] - 1)
            X2[idx, i] += H[idx, phase[i].strand2.haplotypelabel[s]]
        end
        idx = phase[i].strand2.start[end]:phase[i].strand2.length
        X2[idx, i] += H[idx, phase[i].strand2.haplotypelabel[end]]
    end
end

"""
    impute_discard_phase!(X, H, phase, [impute_untyped])

Imputes missing entries of `X` using corresponding haplotypes `H` via `phase`
information. Non-missing entries in `X` will not change, but X and H has to be
aligned. This does NOT preserve phase information.
"""
function impute_discard_phase!(
    X::AbstractMatrix,
    compressed_Hunique::CompressedHaplotypes,
    phase::Vector{HaplotypeMosaicPair},
    impute_untyped::Bool=true
    )

    p, n = size(X)

    # for person in 1:n
    Threads.@threads for person in 1:n
        @inbounds for snp in 1:p
            if ismissing(X[snp, person])
                #find which segment the snp is located
                hap1_segment = searchsortedlast(phase[person].strand1.start,snp)
                hap2_segment = searchsortedlast(phase[person].strand2.start,snp)

                #find haplotype pair in corresponding window for this segment
                h1 = phase[person].strand1.haplotypelabel[hap1_segment]
                h2 = phase[person].strand2.haplotypelabel[hap2_segment]
                w1 = phase[person].strand1.window[hap1_segment]
                w2 = phase[person].strand2.window[hap2_segment]
                i1 = snp - compressed_Hunique.Hstart[w1] + 1
                i2 = snp - compressed_Hunique.Hstart[w2] + 1

                # imputation step
                H1 = impute_untyped ? compressed_Hunique.CW[w1].uniqueH :
                    compressed_Hunique.CW_typed[w1].uniqueH
                H2 = impute_untyped ? compressed_Hunique.CW[w2].uniqueH :
                    compressed_Hunique.CW_typed[w2].uniqueH
                X[snp, person] = H1[i1, h1] + H2[i2, h2]
            end
        end
    end

    return nothing
end

"""
    update_marker_position!(phaseinfo, tgtfile, reffile)

Converts `phaseinfo`'s strand1 and strand2's starting position in
terms of matrix rows of `X` to starting position in terms matrix
rows in `H`.
"""
function update_marker_position!(
    phaseinfo::Vector{HaplotypeMosaicPair},
    XtoH_idx::AbstractVector,
    )
    people = length(phaseinfo)

    @inbounds for j in 1:people
        # update strand1's starting position
        for (i, idx) in enumerate(phaseinfo[j].strand1.start)
            phaseinfo[j].strand1.start[i] = XtoH_idx[idx]
        end
        # update strand2's starting position
        for (i, idx) in enumerate(phaseinfo[j].strand2.start)
            phaseinfo[j].strand2.start[i] = XtoH_idx[idx]
        end
    end

    # update first starting position
    @inbounds for j in 1:people
        phaseinfo[j].strand1.start[1] = 1
        phaseinfo[j].strand2.start[1] = 1
    end

    return nothing
end


text:
#= Given a number N check if it is a power of 2 in O(1). We will be
using bitwise operator to get the calculation.=#

## Function

function isPowerOfTwo(n)
    return (n & (n - 1) == 0)
end

## Input

println("Enter the number : ")
n = readline()
n = parse(Int64, n)

## Calling the function

if (isPowerOfTwo(n))
    println("Its a power of 2")
else
    println("Its not a power of 2")
end

#=
Sample Test Case
Input:
    Enter a number : 4294967296
Output:
    Its a power of 2
Time complexity: O( 1 )
=#


text:
module ThresholdStability

using MathematicalSystems, SwitchOnSafety
using JuMP, SumOfSquares, DynamicPolynomials, MultivariatePolynomials
using CSDP
using Combinatorics, LinearAlgebra

import Reexport
Reexport.@reexport using HybridSystems

const AbstractSwitchedSystem = Union{DiscreteSwitchedLinearSystem, ConstrainedDiscreteSwitchedLinearSystem, StateDepDiscreteSwitchedLinearSystem}

"""
    spectral_radius(A::AbstractMatrix)

Calculates the spectral radius of matrix `A`.
"""
spectral_radius(A::AbstractMatrix) = ρ(A)  # from SwitchOnSafety.jl

include("conversion.jl")

"""
    jsr(Σ::AbstractVector{<:AbstractMatrix}; d = 2)

Calculates an upper bound on the joint spectral radius of a set of matrices `Σ`.

    jsr(s::AbstractSwitchedSystem; d = 2)

Calculates an upper bound on the joint spectral radius for the switched linear system `s`.
"""
jsr(Σ::AbstractVector{<:AbstractMatrix}; d = 2, optimizer = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)) = soslyapb(discreteswitchedsystem(Σ), d, optimizer_constructor = optimizer)[2]
jsr(s::AbstractSwitchedSystem; d = 2, optimizer = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)) = soslyapb(unconstrained(unstatedep(s)), d, optimizer_constructor = optimizer)[2]

"""
    cjsr(Σ::AbstractVector{<:AbstractMatrix}, G::AbstractAutomaton; d = 2)

Calculates an upper bound on the constrained joint spectral radius of the constrained switched system `(Σ, G)`.

    cjsr(s::AbstractSwitchedSystem; d = 2)

Calculates an upper bound on the constrained joint spectral radius of the switched system `s`.
"""
cjsr(Σ::AbstractVector{<:AbstractMatrix}, G::AbstractAutomaton; d = 2, optimizer = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)) = soslyapb(discreteswitchedsystem(Σ, G), d, optimizer_constructor = optimizer)[2]
cjsr(s::AbstractSwitchedSystem; d = 2, optimizer = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)) = soslyapb(unstatedep(s), 2, optimizer_constructor = optimizer)[2]

export spectral_radius, cjsr, jsr

include("indicator.jl")
include("automaton_constructor.jl")

include("cksvar_functions.jl")

include("veronese_lift.jl")
include("midpoints.jl")
include("sdp_lyapunov_feasibility.jl")
include("sdp_gamma_search.jl")
include("sos_lyapunov_feasibility.jl")
include("sos_gamma_search.jl")
end


text:
using IntUtils

# Why can't I write the following on one line?
@enum DIR Right Down Left Up

function next_grid_elt(dir::DIR, r, c)
    if dir == Right
        return (r, c + 1)
    elseif dir == Down
        return (r + 1, c)
    elseif dir == Left
        return (r, c - 1)
    elseif dir == Up
        return (r - 1, c)
    else
        return (0, 0)
    end
end

function next_dir(dir)
    return DIR((Int(dir) + 1) % 4)
end

function lookup_index{T<:Integer}(dim::T, row::T, col::T)::T
    dir = Right
    runlen = 1
    r = div(dim, 2) + 1
    c = div(dim, 2) + 1
    value = 1
    if r == row && c == col
        return value
    end

    runlen = 1
    max_value = dim * dim

    while value < max_value
        for step in 1:2*runlen
            (r, c) = next_grid_elt(dir, r, c)
            value += 1
            if r == row && c == col
                return value
            end
            if step == runlen
                dir = next_dir(dir)
            end
        end
        runlen += 1
        dir = next_dir(dir)
    end
    error(@sprintf("dim=%d, row=%d, col=%d, r=%d, c=%d, value=%d\n",
                   dim, row, col, r, c, value))
end

function main()
    dim = 1001
    ul_diag = 0
    ll_diag = 0
    center = div(dim, 2) + 1
    for d in 1:dim
        ul_diag += lookup_index(dim, d, d)
        ll_diag += lookup_index(dim, dim - d + 1, d)
    end
    v = ul_diag + ll_diag - lookup_index(dim, center , center)
    @printf("array sum over diagonals for dim = %d is %d\n", dim, v)
end

main()


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#profiling speed differences between exponentiation and bitshifting.

const count = 1000000

#enclose in a function body because in julia there are no constant type globals.
function testy()

println("UInt64 1<<3 $(1<<3)")
x::UInt64 = 1
tic()
for(i=1:count)
  x << 3
end
toc()

println("undeclared UInt16 1<<3 $(1<<3)")
z = UInt16(1)
tic()
for(i=1:count)
  z << 3
end
toc()

println("UInt16 1<<3 $(1<<3)")
y::UInt16 = 1
tic()
for(i=1:count)
  y << 3
end
toc()

println("UInt16 2^3 $(2^3)")
q::UInt16 = 2
tic()
for(i=1:count)
  q^3
end
toc()
println("----")

end

#do it three times in case there's something wierd going on with the JIT.
testy()
testy()
testy()

#18 august 2015
#results run on Laptop (Core3 2 Ghz x 2)

#UInt64 1<<3 8
#elapsed time: 1.094e-6 seconds
#undeclared UInt16 1<<3 8
#elapsed time: 1.232e-6 seconds
#UInt16 1<<3 8
#elapsed time: 2.53e-7 seconds
#UInt16 2^3 8
#elapsed time: 0.019976645 seconds
#----
#UInt64 1<<3 8
#elapsed time: 5.39e-7 seconds
#undeclared UInt16 1<<3 8
#elapsed time: 2.75e-7 seconds
#UInt16 1<<3 8
#elapsed time: 2.34e-7 seconds
#UInt16 2^3 8
#elapsed time: 0.020264097 seconds
#----
#UInt64 1<<3 8
#elapsed time: 6.19e-7 seconds
#undeclared UInt16 1<<3 8
#elapsed time: 2.33e-7 seconds
#UInt16 1<<3 8
#elapsed time: 2.35e-7 seconds
#UInt16 2^3 8
#elapsed time: 0.019941948 seconds
#----

#results:  It is a good idea to give the JIT compiler some hints as to what
#you're about to do.  Using UInt16 is much better than using UInt64, bitshifting
#is WAY better than exponentiation.  When profiling, it's a good idea to run
#multiple times, because it seems julia caches its compiled results.

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
"""
    Ensembles

This module provides the main function [`run_ensemble`](@ref Ensembles.run_ensemble).
This serves to run multiple trajectories for a given simulation type, sampling
from an initial distribution.
"""
module Ensembles

using SciMLBase: SciMLBase, EnsembleProblem, solve
using RecursiveArrayTools: ArrayPartition
using TypedTables: TypedTables

using NQCBase: NQCBase
using NQCDynamics:
    AbstractSimulation,
    Simulation,
    RingPolymerSimulation,
    DynamicsUtils,
    DynamicsMethods,
    DynamicsOutputs
using RingPolymerArrays: RingPolymerArray

export run_ensemble

include("selections.jl")
include("outputs.jl")
include("reductions.jl")

"""
    run_ensemble(sim::AbstractSimulation, tspan, distribution;
        selection::Union{Nothing,AbstractVector}=nothing,
        output=(sol,i)->(sol,false),
        reduction=:append
        reduction::Symbol=:append,
        ensemble_algorithm=SciMLBase.EnsembleThreads(),
        algorithm=DynamicsMethods.select_algorithm(sim),
        saveat=[],
        trajectories=1,
        kwargs...
        )

Run multiple trajectories for timespan `tspan` sampling from `distribution`.
The DifferentialEquations ensemble interface is used which allows us to specify
functions to modify the output and how it is reduced across trajectories.

# Keywords

* `selection` should be an `AbstractVector` containing the indices to sample from the `distribution`. By default, `nothing` leads to random sampling.
* `output` can be a function that transforms the DiffEq solution to an output, or a tuple of output quantities as for `run_trajectory`.
* `reduction` defines how the data is reduced across trajectories. Options are `:append`, `:mean` or `:sum`.
* `ensemble_algorithm` is the algorithm from DifferentialEquations which determines which form of parallelism is used.
* `algorithm` is the algorithm used to integrate the equations of motion.
* `saveat` mirrors the DiffEq keyword and labels the time points to save the output.
* `trajectories` is the number of trajectories to perform.

This function wraps the EnsembleProblem from DifferentialEquations and passes the `kwargs`
to the `solve` function.
"""
function run_ensemble(
    sim::AbstractSimulation, tspan, distribution;
    selection::Union{Nothing,AbstractVector}=nothing,
    output=SciMLBase.DEFAULT_OUTPUT_FUNC,
    reduction::Symbol=:append,
    ensemble_algorithm=SciMLBase.EnsembleThreads(),
    algorithm=DynamicsMethods.select_algorithm(sim),
    saveat=[],
    trajectories=1,
    u_init=nothing,
    kwargs...
)

    kwargs = NQCBase.austrip_kwargs(;kwargs...)
    trajectories = convert(Int, trajectories)
    saveat = austrip.(saveat)
    tspan = austrip.(tspan)
    (output isa Symbol) && (output = (output,))
    reduction = Reduction(reduction)

    ensemble_problem = EnsembleProblem(sim, tspan, distribution, selection, output, reduction, u_init, saveat, trajectories, kwargs)

    @info "Performing $trajectories trajectories."
    stats = @timed solve(ensemble_problem, algorithm, ensemble_algorithm; saveat, trajectories, kwargs...)
    @info "Finished after $(stats.time) seconds."

    return extract_output(stats, output, reduction, trajectories)
end

function SciMLBase.EnsembleProblem(sim::AbstractSimulation, tspan, distribution, selection, output, reduction, u_init, saveat, trajectories, kwargs)

    prob_func = Selection(distribution, selection)

    problem = create_problem(sim, distribution, tspan)
    output_func = Output(output, sim)
    if output_func isa DynamicsOutputs.EnsembleSaver
        u_init = get_u_init(reduction, saveat, kwargs, tspan, problem.u0, output_func)
    end

    return EnsembleProblem(problem; prob_func, output_func, reduction, u_init)
end

Output(output::Tuple, sim) = DynamicsOutputs.EnsembleSaver(output, sim)
Output(output, _) = output

function create_problem(sim::AbstractSimulation, distribution, tspan)
    u0 = sample_distribution(sim, distribution)
    return DynamicsMethods.create_problem(u0, tspan, sim)
end

function extract_output(stats, output, reduction, trajectories)
    sol = stats.value

    if output isa Tuple
        process_output(reduction, sol, output, trajectories)
    else
        return sol.u
    end
end

function process_output(::AppendReduction, sol, output, trajectories)
    return [generate_table(output, s) for s in sol.u]
end

function process_output(::SumReduction, sol, output, trajectories)
    for i in eachindex(sol.u)
        sol.u[i][1] /= trajectories
    end
    return generate_table(output, sol.u)
end

function process_output(::MeanReduction, sol, output, trajectories)
    return generate_table(output, sol.u)
end

function generate_table(output, u)
    output = (:t, output...)
    out = [(;zip(output, quantity)...) for quantity in u]
    return TypedTables.Table(out)
end

end # module


text:
"""
    abstract type AbstractOrbital

Abstract supertype of all orbital types.

!!! note "Broadcasting"

    When broadcasting, orbital objects behave like scalars.
"""
abstract type AbstractOrbital end
Base.Broadcast.broadcastable(x::AbstractOrbital) = Ref(x)

"""
    const MQ = Union{Int,Symbol}

Defines the possible types that may represent the main quantum number. It can either be an
non-negative integer or a `Symbol` value (generally used to label continuum electrons).
"""
const MQ = Union{Int,Symbol}

nisless(an::T, bn::T) where T = an < bn
# Our convention is that symbolic main quantum numbers are always
# greater than numeric ones, such that ks appears after 2p, etc.
nisless(an::I, bn::Symbol) where {I<:Integer} = true
nisless(an::Symbol, bn::I) where {I<:Integer} = false

function Base.ascii(o::AbstractOrbital)
    io = IOBuffer()
    ctx = IOContext(io, :ascii=>true)
    show(ctx, o)
    String(take!(io))
end

# * Non-relativistic orbital

"""
    struct Orbital{N <: AtomicLevels.MQ} <: AbstractOrbital

Label for an atomic orbital with a principal quantum number `n::N` and orbital angular
momentum `ℓ`.

The type parameter `N` has to be such that it can represent a proper principal quantum number
(i.e. a subtype of [`AtomicLevels.MQ`](@ref)).

# Properties

The following properties are part of the public API:

* `.n :: N` -- principal quantum number ``n``
* `.ℓ :: Int` -- the orbital angular momentum ``\\ell``

# Constructors

    Orbital(n::Int, ℓ::Int)
    Orbital(n::Symbol, ℓ::Int)

Construct an orbital label with principal quantum number `n` and orbital angular momentum `ℓ`.
If the principal quantum number `n` is an integer, it has to positive and the angular momentum
must satisfy `0 <= ℓ < n`.

```jldoctest
julia> Orbital(1, 0)
1s

julia> Orbital(:K, 2)
Kd
```
"""
struct Orbital{N<:MQ} <: AbstractOrbital
    n::N
    ℓ::Int
    function Orbital(n::Int, ℓ::Int)
        n ≥ 1 || throw(ArgumentError("Invalid principal quantum number $(n)"))
        0 ≤ ℓ && ℓ < n || throw(ArgumentError("Angular quantum number has to be ∈ [0,$(n-1)] when n = $(n)"))
        new{Int}(n, ℓ)
    end
    function Orbital(n::Symbol, ℓ::Int)
        new{Symbol}(n, ℓ)
    end
end

Orbital{N}(n::N, ℓ::Int) where {N<:MQ} = Orbital(n, ℓ)

Base.:(==)(a::Orbital, b::Orbital) =
    a.n == b.n && a.ℓ == b.ℓ

Base.hash(o::Orbital, h::UInt) = hash(o.n, hash(o.ℓ, h))

"""
    mqtype(::Orbital{MQ}) = MQ

Returns the main quantum number type of an [`Orbital`](@ref).
"""
mqtype(::Orbital{MQ}) where MQ = MQ

Base.show(io::IO, orb::Orbital{N}) where N =
    write(io, "$(orb.n)$(spectroscopic_label(orb.ℓ))")

"""
    degeneracy(orbital::Orbital)

Returns the degeneracy of `orbital` which is `2(2ℓ+1)`

# Examples

```jldoctest
julia> degeneracy(o"1s")
2

julia> degeneracy(o"2p")
6
```
"""
degeneracy(orb::Orbital) = 2*(2orb.ℓ + 1)

"""
    isless(a::Orbital, b::Orbital)

Compares the orbitals `a` and `b` to decide which one comes before the
other in a configuration.

# Examples

```jldoctest
julia> o"1s" < o"2s"
true

julia> o"1s" < o"2p"
true

julia> o"ks" < o"2p"
false
```
"""
function Base.isless(a::Orbital, b::Orbital)
    nisless(a.n, b.n) && return true
    a.n == b.n && a.ℓ < b.ℓ && return true
    false
end

"""
    parity(orbital::Orbital)

Returns the parity of `orbital`, defined as `(-1)^ℓ`.

# Examples

```jldoctest
julia> parity(o"1s")
even

julia> parity(o"2p")
odd
```

"""
parity(orb::Orbital) = p"odd"^orb.ℓ

"""
    symmetry(orbital::Orbital)

Returns the symmetry for `orbital` which is simply `ℓ`.
"""
symmetry(orb::Orbital) = orb.ℓ

"""
    isbound(::Orbital)

Returns `true` is the main quantum number is an integer, `false`
otherwise.

```jldoctest
julia> isbound(o"1s")
true

julia> isbound(o"ks")
false
```
"""
function isbound end
isbound(::Orbital{Int}) = true
isbound(::Orbital{Symbol}) = false

"""
    angular_momenta(orbital)

Returns the angular momentum quantum numbers of `orbital`.

# Examples
```jldoctest
julia> angular_momenta(o"2s")
(0, 1/2)

julia> angular_momenta(o"3d")
(2, 1/2)
```
"""
angular_momenta(orbital::Orbital) = (orbital.ℓ,half(1))
angular_momentum_labels(::Orbital) = ("ℓ","s")

"""
    angular_momentum_ranges(orbital)

Return the valid ranges within which projections of each of the
angular momentum quantum numbers of `orbital` must fall.

# Examples
```jldoctest
julia> angular_momentum_ranges(o"2s")
(0:0, -1/2:1/2)

julia> angular_momentum_ranges(o"4f")
(-3:3, -1/2:1/2)
```
"""
angular_momentum_ranges(orbital::AbstractOrbital) =
    map(j -> -j:j, angular_momenta(orbital))

# ** Saving/loading

Base.write(io::IO, o::Orbital{Int}) = write(io, 'i', o.n, o.ℓ)
Base.write(io::IO, o::Orbital{Symbol}) = write(io, 's', sizeof(o.n), o.n, o.ℓ)

function Base.read(io::IO, ::Type{Orbital})
    kind = read(io, Char)
    n = if kind == 'i'
        read(io, Int)
    elseif kind == 's'
        b = Vector{UInt8}(undef, read(io, Int))
        readbytes!(io, b)
        Symbol(b)
    else
        error("Unknown Orbital type $(kind)")
    end
    ℓ = read(io, Int)
    Orbital(n, ℓ)
end

# * Orbital construction from strings

parse_orbital_n(m::RegexMatch,i=1) =
    isnumeric(m[i][1]) ? parse(Int, m[i]) : Symbol(m[i])

function parse_orbital_ℓ(m::RegexMatch,i=2)
    ℓs = strip(m[i], ['[',']'])
    if isnumeric(ℓs[1])
        parse(Int, ℓs)
    else
        ℓi = findfirst(ℓs, spectroscopic)
        isnothing(ℓi) && throw(ArgumentError("Invalid spectroscopic label: $(m[i])"))
        first(ℓi) - 1
    end
end

function Base.parse(::Type{<:Orbital}, orb_str)
    m = match(r"^([0-9]+|.)([a-z]|\[[0-9]+\])$", orb_str)
    isnothing(m) && throw(ArgumentError("Invalid orbital string: $(orb_str)"))
    n = parse_orbital_n(m)
    ℓ = parse_orbital_ℓ(m)
    Orbital(n, ℓ)
end

"""
    @o_str -> Orbital

A string macro to construct an [`Orbital`](@ref) from the canonical string representation.

```jldoctest
julia> o"1s"
1s

julia> o"Fd"
Fd
```
"""
macro o_str(orb_str)
    parse(Orbital, orb_str)
end

function orbitals_from_string(::Type{O}, orbs_str::AbstractString) where {O<:AbstractOrbital}
    map(split(orbs_str)) do orb_str
        m = match(r"^([0-9]+|.)\[([a-z]|[0-9]+)(-([a-z]|[0-9]+)){0,1}\]$", strip(orb_str))
        m === nothing && throw(ArgumentError("Invalid orbitals string: $(orb_str)"))
        n = parse_orbital_n(m)
        ℓs = map(filter(i -> !isnothing(m[i]), [2,4])) do i
            parse_orbital_ℓ(m, i)
        end
        orbs = if O == RelativisticOrbital
            orbs = map(ℓ -> O(n, ℓ, ℓ-1//2), max(first(ℓs),1):last(ℓs))
            append!(orbs, map(ℓ -> O(n, ℓ, ℓ+1//2), first(ℓs):last(ℓs)))
        else
            map(ℓ -> O(n, ℓ), first(ℓs):last(ℓs))
        end
        sort(orbs)
    end |> o -> vcat(o...) |> sort
end

"""
    @os_str -> Vector{Orbital}

Can be used to easily construct a list of [`Orbital`](@ref)s.

# Examples

```jldoctest
julia> os"5[d] 6[s-p] k[7-10]"
7-element Vector{Orbital}:
 5d
 6s
 6p
 kk
 kl
 km
 kn
```
"""
macro os_str(orbs_str)
    orbitals_from_string(Orbital, orbs_str)
end

export AbstractOrbital, Orbital,
    @o_str, @os_str,
    degeneracy, symmetry, isbound, angular_momenta, angular_momentum_ranges