text:
export buildnet

function buildnet(;name="temp",model=NetworkModel(),method="Euler1",duration=0.1)
  iteration=round(Int64,duration/model.dt)
  syncIter=model.synctime>0?round(Int64,model.synctime÷model.dt):1
  iteration=iteration÷syncIter

  # open a file
  fid=open(name*".jl","w")

  str="""
  println("\\nExecuting $(name).jl")

  using JLD

  if $(model.nCores)>0 && (addCores=$(model.nCores)-nworkers())>0
    nworkers()==1?addprocs(addCores):addprocs(addCores-1)
  end
  println("Using \$(nworkers()) workers")
  """
  println(fid,str)

  # units
  unit="""
  metre=1e2;meter=1e2;cm=metre/1e2;mm=metre/1e3;um=metre/1e6;nm=metre/1e9;
  second=1e3;ms=second/1e3;
  Hz=1/second;
  voltage=1e3;mV=voltage/1e3;
  ampere=1e6;mA=ampere/1e3;uA=ampere/1e6;nA=ampere/1e9;
  farad=1e6;uF=ufarad=farad/1e6;
  siemens=1e3;mS=msiemens=siemens/1e3;nS=nsiemens=siemens/1e9"""
  println(fid,"@everywhere begin\n",unit,"\nend\n")

  # parameters
  println(fid,"@everywhere begin\n",model.parameters,"\nend\n")

  # equations
  eqs=map(strip,split(replace(model.equations,";","\n"),"\n"));
  vars=[]
  diffEqs=[]

  for k=1:length(eqs)
    if eqs[k][1]=='d'
      vars=[vars;match(r"(?<=d).*(?=/dt)",eqs[k]).match]
      diffEq=vectorize(match(r"(?<==).*$",eqs[k]).match)
      diffEqs=[diffEqs;diffEq]
    else
      funcName=match(r"^.*(?==)",eqs[k]).match
      funcExpr=vectorize(match(r"(?<==).*$",eqs[k]).match)
      diffEqs=map(diffEq->replace(diffEq,funcName,funcExpr),diffEqs)
    end
  end

  varsType=map(var->var*"::Array{Float64}",vars)
  argumentType=join([varsType;"t::Float64"],",");
  println(fid,"# function")
  for k=1:length(vars)
    println(fid,"@everywhere F",vars[k],"(",argumentType,")=",diffEqs[k])
  end
  println(fid)

  # function
  println(fid,"function $name(iteration::Int64=10)\n")

  # Memory use before allocation
  str="""
  if Sys.OS_NAME == :Linux
    println("Memory use before allocation:")
    run(pipeline(`free -og`, `grep -v Swap`))
  elseif Sys.OS_NAME == :Windows
    # println("Memory use before allocation:")
    # run(`wmic OS get FreePhysicalMemory /Value`)
  end
  """
  println(fid,str)

  println(fid,"t=0.0")
  println(fid,"dt=",model.dt)
  println(fid,"N=",model.groups[:N])
  println(fid,"syncIter=",syncIter,"\n")

  # initialize variables
  println(fid,"# initialize variables")
  for var in vars
    println(fid,var,"=SharedArray(Float64,N)")
  end
  println(fid,model.initialization,"\n")

  # connection
  println(fid,"# initialize connections\n",model.connection,"\n")

  # report
  report="""# report
  if k%(1000÷syncIter)==0
     @printf("%s%4.2f finished\\n","%",100*k/iteration)
  end
  """

  # integrate

  argument=join([map(var->var*"[localindex]",vars);"t"],",")
  splitArg=split(argument,",")
  argumentTemp=join([splitArg[1]*"+Δ$(vars[1])";splitArg[2:end]],",")

  func=map(var->"F$var($argument)",vars)
  tempFunc="F$(vars[1])($argumentTemp)"

  if model.solver=="Euler1"
    integrate="Δ$(vars[1])=dt*$(func[1])\n"
  elseif model.solver=="Euler2"
    integrate="Δ$(vars[1])=dt/2.0*$(func[1])\n"
    integrate*="Δ$(vars[1])+=dt/2.0*$tempFunc\n"
  elseif model.solver=="Euler4"
    integrate="Δ$(vars[1])=dt/4.0*$(func[1])\n"
    integrate*="Δ$(vars[1])+=dt/4.0*$tempFunc\n"
    integrate*="Δ$(vars[1])+=dt/4.0*$tempFunc\n"
    integrate*="Δ$(vars[1])+=dt/4.0*$tempFunc\n"
  end

  for k=2:length(vars)
    integrate*="Δ$(vars[k])=dt*$(func[k])\n"
  end

  for k=1:length(vars)
	integrate*="$(vars[k])[localindex]+=Δ$(vars[k])\n"
  end
  integrate="""# integrate
  @sync for p in workers()
  @spawnat p begin
  for idx=1:syncIter
  localindex=localindexes($(vars[1]))
  $integrate
  end #end of sync
  end #end of begin
  end #end of workers
  """

  for var in vars
  reg= Regex("(?<=\\b)($var)(?=\\b)")
  diffEqs=map(diffEq->replace(diffEq,reg,s"\1[localindex]"),diffEqs)
  end

  # spike
  if !isempty(model.spike)
    spike=split(replace(model.spike,";","\n"),"\n")
    threshold=spike[1]
    length(spike)==2? refractory=spike[2]:refractory=""

    println(fid,"# initialize tSpike")
    println(fid,"tSpike=zeros(round(Int64,N*iteration*dt))")
    println(fid,"ts=zeros(N)")
    println(fid,"spikeNeuron=zeros(Int64,round(Int64,N*iteration*dt))")
    println(fid,"state=1","\n")

    spike="# spike\n"
    spike*="isExceedThreshold=$threshold\n"
    if !isempty(refractory)
      spike*="isNotRefractory=$refractory\n"
      spike*="isSpike=isExceedThreshold & isNotRefractory\n"
    else
      spike*="isSpike=isExceedThreshold\n"
    end
    spike*="spikeIndex=find(isSpike)\n"
    spike*="nSpike=length(spikeIndex)\n"
    spike*="tSpike[state:state+nSpike-1]=t\n"
    spike*="ts[spikeIndex]=t\n"
    spike*="spikeNeuron[state:state+nSpike-1]=spikeIndex\n"
    spike*="state+=nSpike\n"
  else spike=""
  end

  # reset
  if !isempty(model.reset)
    reset="# reset\n"
    for line in split(replace(model.reset,";","\n"),"\n")
      var=match(r"^.*(?==)",line).match
      expr=match(r"(?<==).*",line).match
      reset*="$(var)[spikeIndex]=$(expr)\n"
    end
  else reset=""
  end

  # synapse
  if !isempty(model.synapse)
    synapse="# synapse\n"
    for line in split(replace(model.synapse,";","\n"),"\n")
      var1=match(r"^.*(?=\+)",line).match
      var2=match(r"(?<==).*",line).match
      if model.nCores>0
      synapse*="$(var1)[:]+=sum($(var2)[:,spikeIndex],2)\n"
      else synapse*="$(var1)+=sum($(var2)[:,spikeIndex],2)\n"
      end
    end
  else synapse=""
  end

  # plasticity
  if !isempty(model.plasticity)
    plasticity="# plasticity\n"
    plasticity*=replace(replace(model.plasticity,"j",":"),"i","spikeIndex")*"\n"
  else plasticity=""
  end

  # record
  if !isempty(model.record)
    println(fid,"# initialize Record")
    for k=1:length(vars)
      println(fid,"$(vars[k])Record=zeros($(length(model.record)),iteration)")
    end
    println(fid)

    record="# record\n"
    for k=1:length(vars)
      record*="$(vars[k])Record[:,k]=$(vars[k])[$(model.record)]\n"
    end
  else
    record=""
  end

  # Memory use after allocation
  str="""
  if Sys.OS_NAME == :Linux
    println("Memory use after allocation:")
    run(pipeline(`free -og`, `grep -v Swap`))
  elseif Sys.OS_NAME == :Windows
    # println("Memory use before allocation:")
    # run(`wmic OS get FreePhysicalMemory /Value`)
  end
  """
  println(fid,str)
  # print status
  println(fid,"# print status")
  str="""
  if iteration>=50
      println("simulation starts...")
  else
    println("compiling...")
  end"""
  println(fid,str,"\n")

  # for-loop
  println(fid,"tic()")
  println(fid,"@fastmath @inbounds for k=1:iteration\n")
  println(fid,report)
  println(fid,integrate)
  println(fid,spike)
  println(fid,reset)
  println(fid,synapse)
  println(fid,plasticity)
  println(fid,record)
  println(fid,"t+=syncIter*dt")
  println(fid,"end # end of for loop ")
  println(fid,"toc()","\n")

  # save variables
  str="""# save variables
  try jldopen("$name.jld","w") do fid
    try  write(fid,"spikeNeuron",spikeNeuron[1:state-1])  end
    try  write(fid,"tSpike",tSpike[1:state-1])            end
    try  write(fid,"$(vars[1])Record",$(vars[1])Record)   end
    """
  for var in model.saveVars
  str*="try write(fid,\"$(var)\",$var) end\n"
  end
  str*="end #end of do\n"
  str*="end #end of try jldopen\n"
  println(fid,str)
    # end of function
  println(fid,"end # end of function\n")

  # precompile and run
  println(fid,"# precompile and run")
  println(fid,"precompile($name,($iteration,))")
  println(fid,"$name(20)")
  println(fid,"$name(20)")
  println(fid,"$name($iteration)")

  # post-processing
  println(fid,"\n",model.postProcs)

  # close file, open file
  close(fid)
  try run(`atom $(abspath("$name.jl"))`) end
  # include("$name.jl")
end


text:
using Optim, LineSearches
using LinearAlgebra: I, norm

const σx = Float64[0 1; 1 0]
const σy = ComplexF64[0 -1im; 1im 0]
const σz = Float64[1 0; 0 -1]
const id2 = Float64[1 0; 0 1]

"""
    diaglocalhamiltonian(diag::Vector)

return the 2-site Hamiltonian with single-body terms given
by the diagonal `diag`.
"""
function diaglocalhamiltonian(diag::Vector)
    n = length(diag)
    h = ein"i -> ii"(diag)
    id = Matrix(I,n,n)
    reshape(h,n,n,1,1) .* reshape(id,1,1,n,n) .+ reshape(h,1,1,n,n) .* reshape(id,n,n,1,1)
end

"""
    energy(h, ipeps; χ, tol, maxit)

return the energy of the `ipeps` 2-site hamiltonian `h` and calculated via a
ctmrg with parameters `χ`, `tol` and `maxit`.
"""
function energy(h::AbstractArray{T,4}, ipeps::IPEPS; χ::Int, tol::Real, maxit::Int) where T
    ipeps = indexperm_symmetrize(ipeps)  # NOTE: this is not good
    D = getd(ipeps)^2
    s = gets(ipeps)
    ap = ein"abcdx,ijkly -> aibjckdlxy"(ipeps.bulk, conj(ipeps.bulk))
    ap = reshape(ap, D, D, D, D, s, s)
    a = ein"ijklaa -> ijkl"(ap)

    rt = SquareCTMRGRuntime(a, Val(:raw), χ)
    rt  = ctmrg(rt; tol=tol, maxit=maxit)
    e = expectationvalue(h, ap, rt)
    return e
end

"""
    expectationvalue(h, ap, rt)

return the expectationvalue of a two-site operator `h` with the sites
described by rank-6 tensor `ap` each and an environment described by
a `SquareCTMRGRuntime` `rt`.
"""
function expectationvalue(h, ap, rt::SquareCTMRGRuntime) where T
    corner, edge = rt.corner, rt.edge
    ap /= norm(ap)
    l = ein"ab,ica,bde,cjfdlm,eg,gfk -> ijklm"(corner,edge,edge,ap,corner,edge)
    e = ein"abcij,abckl,ijkl -> "(l,l,h)[]
    n = ein"ijkaa,ijkbb -> "(l,l)[]
    return e/n
end

"""
    optimiseipeps(ipeps, h; χ, tol, maxit, optimargs = (), optimmethod = LBFGS(m = 20))

return the tensor `bulk'` that describes an ipeps that minimises the energy of the
two-site hamiltonian `h`. The minimization is done using `Optim` with default-method
`LBFGS`. Alternative methods can be specified by loading `LineSearches` and
providing `optimmethod`. Other options to optim can be passed with `optimargs`.
The energy is calculated using ctmrg with parameters `χ`, `tol` and `maxit`.
"""
function optimiseipeps(ipeps::IPEPS{LT}, h; χ::Int, tol::Real, maxit::Int,
        optimargs = (),
        optimmethod = LBFGS(m = 20)) where LT
    bulk = ipeps.bulk
    let energy = x -> real(energy(h, IPEPS{LT}(x); χ=χ, tol=tol, maxit=maxit))
        res = optimize(energy,
            Δ -> Zygote.gradient(energy,Δ)[1], bulk, optimmethod, inplace = false, optimargs...)
    end
end

ext: jl


text:
import Base.==

# at some point these should be a pull request for Cairo.jl?

# this is a #define - no idea how to extract it dynamically
# https://developer.gnome.org/pango/stable/pango-Glyph-Storage.html#PANGO-SCALE:CAPS
const PANGO_SCALE = 1024

# these are enums - again, no idea how to read
# https://people.redhat.com/otaylor/pango-mirror/api/pango-fonts.html#PANGOSTYLE
const PANGO_STYLE_NORMAL = 0
const PANGO_STYLE_ITALIC = 1
const PANGO_STYLE_OBLIQUE = 2
# https://people.redhat.com/otaylor/pango-mirror/api/pango-fonts.html#PANGOVARIANT
const PANGO_VARIANT_NORMAL = 0
const PANGO_VARIANT_SMALL_CAPS = 1
# https://developer.gnome.org/pygtk/stable/class-pangofontdescription.html#method-pangofontdescription--set-weight
# any value from 100-900
const PANGO_WEIGHT_ULTRALIGHT = 200
const PANGO_WEIGHT_LIGHT = 300
const PANGO_WEIGHT_NORMAL = 400
const PANGO_WEIGHT_BOLD = 700
const PANGO_WEIGHT_ULTRABOLD = 800
const PANGO_WEIGHT_HEAVY = 900
# https://people.redhat.com/otaylor/pango-mirror/api/pango-fonts.html#PANGOSTRETCH
const PANGO_STRETCH_ULTRA_CONDENSED = 0
const PANGO_STRETCH_EXTRA_CONDENSED = 1
const PANGO_STRETCH_CONDENSED = 2
const PANGO_STRETCH_SEMI_CONDENSED = 3
const PANGO_STRETCH_NORMAL = 4
const PANGO_STRETCH_SEMI_EXPANDED = 5
const PANGO_STRETCH_EXPANDED = 6
const PANGO_STRETCH_EXTRA_EXPANDED = 7
const PANGO_STRETCH_ULTRA_EXPANDED = 8
# http://sourcecodebrowser.com/pango1.0/1.22.1/pango-gravity_8h.html#ab9f9c35c3778c1ca07195d6687a83e9c
const PANGO_GRAVITY_SOUTH = 0
const PANGO_GRAVITY_EAST = 1
const PANGO_GRAVITY_NORTH = 2
const PANGO_GRAVITY_WEST = 3
const PANGO_GRAVITY_AUTO = 4

function has_current_point(ctx::X.CairoContext)
    ccall((:cairo_has_current_point, X._jl_libcairo), Int, (Ptr{Void},), ctx.ptr) != 0
end

# can be called if has_current_point is false (will return 0.0, 0.0)
function get_current_point(ctx::X.CairoContext)
    x = Cdouble[0]
    y = Cdouble[0]
    ccall((:cairo_get_current_point, X._jl_libcairo), Void, (Ptr{Void}, Ptr{Cdouble}, Ptr{Cdouble}), ctx.ptr, x, y)
    (x[1], y[1])
end

immutable FontMap
    ptr::Ptr{Void}
end

function get_font_map_default()
    FontMap(ccall((:pango_cairo_font_map_get_default, X._jl_libpangocairo), Ptr{Void}, ()))
end

immutable FontFamily
    ptr::Ptr{Void}
end

# this assumes the underlying structs that are pointed to are
# persistent which seems to be implied by the api.  so we copy the
# array of pointers (maybe that could be avoied by calling
# pointer_to_array with true)
function list_families(fm::FontMap)
    pff = Ptr{Ptr{Void}}[0]
    pn = Cint[0]
    ccall((:pango_font_map_list_families, X._jl_libpango), Void, (Ptr{Void}, Ptr{Ptr{Ptr{Void}}}, Ptr{Cint}), fm.ptr, pff, pn)
    a = pointer_to_array(pff[1], (pn[1],), false)
    ff = [FontFamily(p) for p in a]
    Libc.free(pff[1])
    ff
end

function get_name(ff::FontFamily)
    bytestring(ccall((:pango_font_family_get_name, X._jl_libpango), Cstring, (Ptr{Void},), ff.ptr))
end

function is_monospace(ff::FontFamily)
    convert(Bool, ccall((:pango_font_family_is_monospace, X._jl_libpango), Cint, (Ptr{Void},), ff.ptr))
end

immutable FontFace
    ptr::Ptr{Void}
end

# this assumes the underlying structs that are pointed to are
# persistent which seems to be implied by the api.  so we copy the
# array of pointers (maybe that could be avoied by calling
# pointer_to_array with true)
function list_faces(ff::FontFamily)
    pfc = Ptr{Ptr{Void}}[0]
    pn = Cint[0]
    ccall((:pango_font_family_list_faces, X._jl_libpango), Void, (Ptr{Void}, Ptr{Ptr{Ptr{Void}}}, Ptr{Cint}), ff.ptr, pfc, pn)
    a = pointer_to_array(pfc[1], (pn[1],), false)
    fc = [FontFace(p) for p in a]
    Libc.free(pfc[1])
    fc
end

function get_name(fc::FontFace)
    bytestring(ccall((:pango_font_face_get_face_name, X._jl_libpango), Cstring, (Ptr{Void},), fc.ptr))
end

# cannot use immutable here as finalizer() complains
type FontDescription
    ptr::Ptr{Void}
    function FontDescription(ptr)
        this = new(ptr)
        finalizer(this, fd -> ccall((:pango_font_description_free, X._jl_libpango), Void, (Ptr{Void},), fd.ptr))
        this
    end
end

function describe(fc::FontFace)
    FontDescription(ccall((:pango_font_face_describe, X._jl_libpango), Ptr{Void}, (Ptr{Void},), fc.ptr))
end

# non-standard function names here to support round-tripping to strings

function Base.show(io::IO, fd::FontDescription)
    s = bytestring(ccall((:pango_font_description_to_string, X._jl_libpango), Cstring, (Ptr{Void},), fd.ptr))
    print(io, s)
end

function Base.convert(::Type{FontDescription}, s::AbstractString)
    FontDescription(ccall((:pango_font_description_from_string, X._jl_libpango), Ptr{Void}, (Ptr{Cchar},), s))
end

function ==(fd1::FontDescription, fd2::FontDescription)
    convert(Bool, ccall((:pango_font_description_equal, X._jl_libpango), Cint, (Ptr{Void}, Ptr{Void}), fd1.ptr, fd2.ptr))
end

function Base.hash(fd::FontDescription)
    ccall((:pango_font_description_hash, X._jl_libpango), Culong, (Ptr{Void},), fd.ptr)
end

function Base.copy(fd::FontDescription)
    FontDescription(ccall((:pango_font_description_copy, X._jl_libpango), Ptr{Void}, (Ptr{Void},), fd.ptr))
end

for name in [:style, :variant, :weight, :stretch, :gravity, :size]
    g1 = symbol(:get_, name)
    g2 = symbol(:pango_font_description_, g1)
    s1 = symbol(:set_, name)
    s2 = symbol(:pango_font_description_, g1)
    @eval begin
        $g1(fd::FontDescription) = ccall(($(Expr(:quote, g2)), X._jl_libpango), Culong, (Ptr{Void},), fd.ptr)
        $s1(fd::FontDescription, x) = ccall(($(Expr(:quote, s2)), X._jl_libpango), Void, (Ptr{Void}, Culong), fd.ptr, x)
    end
end

function set_absolute_size(fd::FontDescription, size)
    ccall((:pango_font_description_set_absolute_size, X._jl_libpango), Void, (Ptr{Void}, Cdouble), fd.ptr, size * PANGO_SCALE)
end

function get_absolute_size(fd::FontDescription, size)
    ccall((:pango_font_description_get_absolute_size, X._jl_libpango), Cdouble, (Ptr{Void},), fd.ptr) / PANGO_SCALE
end

immutable Layout
    ptr::Ptr{Void}
end

function set_text(l::Layout, text)
    ccall((:pango_layout_set_text, X._jl_libpango), Void, (Ptr{Void}, Ptr{Cchar}, Cint), l.ptr, text, length(text))
end

function set_description(l::Layout, fd::FontDescription)
    ccall((:pango_layout_set_font_description, X._jl_libpango), Void, (Ptr{Void}, Ptr{Void}), l.ptr, fd.ptr)
end

function update_layout(c::X.CairoContext, l::Layout)
    ccall((:pango_cairo_update_layout, X._jl_libpangocairo), Void, (Ptr{Void}, Ptr{Void}), c.ptr, l.ptr)
end

function show_path(c::X.CairoContext, l::Layout)
    ccall((:pango_cairo_layout_path, X._jl_libpangocairo), Void, (Ptr{Void}, Ptr{Void}), c.ptr, l.ptr)
end

immutable Rectangle
    x::Cint
    y::Cint
    width::Cint
    height::Cint
end

function get_pixel_extents(l::Layout)
    ink = Rectangle[Rectangle(0, 0, 0, 0)]
    logical = Rectangle[Rectangle(0, 0, 0, 0)]
    ccall((:pango_layout_get_pixel_extents, X._jl_libpango), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void}), l.ptr, ink, logical)
    ink[1], logical[1]
end


text:
abstract type AbstractMapper end

mapper(mesh) = Dict((c,i) for (i,c) in enumerate(cells(mesh)))


text:
using Documenter
using PackedVectorsOfVectors

DocMeta.setdocmeta!(
    PackedVectorsOfVectors,
    :DocTestSetup,
    quote
        using PackedVectorsOfVectors
        using PackedVectorsOfVectors: cumsum1
    end
)

makedocs(
    modules = PackedVectorsOfVectors,
    sitename="PackedVectorsOfVectors.jl"
)
deploydocs(repo = "github.com/synchronoustechnologies/PackedVectorsOfVectors.jl.git")


text:
# The basis definition
@variables u[1:2]
basis = Basis([polynomial_basis(u, 5); sin.(u); cos.(u)], u)

function pendulum(u, p, t)
    x = u[2]
    y = -9.81sin(u[1]) - 0.1u[2]^3 - 0.2 * cos(u[1])
    return [x;y]
end

u0 = [0.99π; -1.0]
tspan = (0.0, 20.0)
dt = 0.1
prob = ODEProblem(pendulum, u0, tspan)
sol = solve(prob, Tsit5(), saveat=dt)

X = sol[:,:]
t = sol.t

DX = similar(sol[:,:])

for (i, xi) in enumerate(eachcol(sol[:,:]))
    DX[:,i] = pendulum(xi, [], 0.0)
end

##
@testset "Ideal data" begin
    dd_prob = ContinuousDataDrivenProblem(
        sol
        )

    opts = [
    STLSQ(1e-2), ADMM(1e-2, 1e-2), SR3(1e-2, SoftThreshold()),
    SR3(1e-4)
    ]

res = solve(dd_prob, basis, opts[1], maxiter = 10000)

    for opt in opts
        m = DataDrivenDiffEq.metrics(res)
        @test all(m[:L₂] .< 1e-1*size(X, 2))
        @test all(m[:AIC] .>= 200) # Perfect Match or close
        @test all(m[:R²] .>= 0.9)
    end
end

Random.seed!(1234)
X = X .+ 1e-1*randn(size(X))

@testset "Noisy data" begin

    dd_prob_noisy = ContinuousDataDrivenProblem(
        X, t, GaussianKernel()
        )

    opts = [
            STLSQ(1e-1:5e-1:1e3), ADMM(1e-1:5e-1:1e5, 1.0), SR3(1e-1:5e-1:1e5, 2.0, SoftThreshold()),
            SR3((1e-1:5e-1:1e5), 2.0, HardThreshold()), SR3(1e-3:5e-1:1e5, 2.0, ClippedAbsoluteDeviation())
        ]

    for opt in opts
        res = solve(dd_prob_noisy, basis, opt, maxiter = 50000, denoise = true, normalize = true)
        m = DataDrivenDiffEq.metrics(res)
        @test all(m[:L₂]./size(X,2) .< [30.; 800])
        @test all(-500.0 .< m[:AIC])
        @test all(m[:R²] .>= 0.9)
    end

end


text:
include("MMKP.jl")

mp(a, b) = joinpath("MMKP", "benchmark_problems", a, b * ".txt")
flat(a) = vcat(a...)

test_cases = [
	['A', "Class A", ["I07", "I13"]],
	['B', "Class B", ["INST01", "INST11", "INST20"]],
	['C', "Class C", ["INST21", "INST30"]],
	['D', "Class D", ["Instance1", "Instance256"]],
	['D', "Class D", ["Instance100"]]
]
test_problems = flat(map(tc->[[tc[1], mp(tc[2], file)] for file in tc[3]], test_cases))

loaders = Dict('A'=>load_catA, 'B'=>load_catB, 'C'=>load_catC, 'D'=>load_catD)

function load_model(tc)
	println("testing $(tc[2])")
	problem = loaders[tc[1]](tc[2])
	model = jump_model(problem, 4)
	model
end

models = load_model.(test_problems)

load_model(test_problems[end])

function test_optimize(m)
	set_optimizer(m, CPLEX.Optimizer)
	JOH.Matheur.set_time!(m, 5)
	optimize!(m)
	m
end

test_optimize.(models)

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
# ----------------------------------------
# FUNCTIONS
# ----------------------------------------

function distance_matrix(left, right; deletion = 1, insertion = 1, substitution = 1)
    function complete(left, right)
        distance = zeros(Int64, length(left) + 1, length(right) + 1)

        for l ∈ 1:length(left)
            @inbounds distance[l + 1, 1] = l * deletion
        end

        for r ∈ 1:length(right)
            @inbounds distance[1, r + 1] = r * insertion
        end

        for l ∈ 1:length(left)
            for r ∈ 1:length(right)
                @inbounds distance[l + 1, r + 1] = min(distance[l, r + 1] + deletion,
                                                       distance[l + 1, r] + insertion,
                                                       distance[l, r] + (left[l] == right[r] ? 0 : substitution))
            end
        end

        return distance
    end

    if length(left) == 0
        return length(right)

    elseif length(right) == 0
        return length(left)

    else
        return complete(left, right)[end]
    end
end

max_stars_repo_name: jlegare/LevenshteinToolkit

max_issues_repo_name: jlegare/LevenshteinToolkit

max_forks_repo_name: jlegare/LevenshteinToolkit


text:
function Ephemeris(arg0::List, arg1::jint)
    return Ephemeris((List, jint), arg0, arg1)
end

function Ephemeris(arg0::List, arg1::jint, arg2::jdouble)
    return Ephemeris((List, jint, jdouble), arg0, arg1, arg2)
end

function Ephemeris(arg0::List, arg1::jint, arg2::jdouble, arg3::AttitudeProvider)
    return Ephemeris((List, jint, jdouble, AttitudeProvider), arg0, arg1, arg2, arg3)
end

function basic_propagate(obj::Ephemeris, arg0::AbsoluteDate)
    return jcall(obj, "basicPropagate", SpacecraftState, (AbsoluteDate,), arg0)
end

function get_extrapolation_threshold(obj::Ephemeris)
    return jcall(obj, "getExtrapolationThreshold", jdouble, ())
end

function get_frame(obj::Ephemeris)
    return jcall(obj, "getFrame", Frame, ())
end

function get_initial_state(obj::Ephemeris)
    return jcall(obj, "getInitialState", SpacecraftState, ())
end

function get_managed_additional_states(obj::Ephemeris)
    return jcall(obj, "getManagedAdditionalStates", Vector{JString}, ())
end

function get_max_date(obj::Ephemeris)
    return jcall(obj, "getMaxDate", AbsoluteDate, ())
end

function get_min_date(obj::Ephemeris)
    return jcall(obj, "getMinDate", AbsoluteDate, ())
end

function get_pv_coordinates(obj::Ephemeris, arg0::AbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, (AbsoluteDate, Frame), arg0, arg1)
end

function is_additional_state_managed(obj::Ephemeris, arg0::JString)
    return jcall(obj, "isAdditionalStateManaged", jboolean, (JString,), arg0)
end

function reset_initial_state(obj::Ephemeris, arg0::SpacecraftState)
    return jcall(obj, "resetInitialState", void, (SpacecraftState,), arg0)
end


text:
module TestLexAndParse
    using Base.Test
    using Brainfuck

    s1 = "<-->"

    @test Brainfuck.lexandparse(s1) == [
        Brainfuck.OP2,
        Brainfuck.OP4,
        Brainfuck.OP4,
        Brainfuck.OP1,
    ]

    s2 = """
        >--< ++ -- ><
        >--< ++ -- ><
    """

    @test Brainfuck.lexandparse(s2) == [
        Brainfuck.OP1,
        Brainfuck.OP4,
        Brainfuck.OP4,
        Brainfuck.OP2,
        Brainfuck.OP3,
        Brainfuck.OP3,
        Brainfuck.OP4,
        Brainfuck.OP4,
        Brainfuck.OP1,
        Brainfuck.OP2,

        Brainfuck.OP1,
        Brainfuck.OP4,
        Brainfuck.OP4,
        Brainfuck.OP2,
        Brainfuck.OP3,
        Brainfuck.OP3,
        Brainfuck.OP4,
        Brainfuck.OP4,
        Brainfuck.OP1,
        Brainfuck.OP2,
    ]
end


text:
# Autogenerated wrapper script for FSPS_jll for x86_64-linux-gnu-libgfortran5
export fsps

using sed_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("FSPS")
JLLWrappers.@declare_library_product(fsps, "fsps.so")
function __init__()
    JLLWrappers.@generate_init_header(sed_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        fsps,
        "lib/fsps.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
include("../src/Includes.jl")

function thread_printer(this::Thread)
  println("thread_printer")
  i=0
  while true
    global Messages

    println(this, i ;title="PRINT")
    showMessages()

    sleep(1)
    i+=1
  end
  nothing
end

function thread_compute(this::Thread)
  println("thread_compute")
  timeRef = Ref(0.0)
  i=0
  while true
    if OnTime(0.25, timeRef)
      info(this, i ;title="CALC")
      i+=1
    end
    sleep(0.001)
  end
end

function thread_renderer(this::Thread)
  println("thread_renderer")
  timeRef = Ref(0.0)
  i=0
  while true
    if OnTime(0.5, timeRef)
      debug(this, i ;title="RENDER")
      i+=1
    end
    sleep(0.01)
  end
  nothing
end

function thread_sound(this::Thread)
  println("thread_sound")
  timeRef = Ref(0.0)
  i=0
  while true
    if OnTime(0.75, timeRef)
      warn(this, i ;title="SOUND")
      i+=1
    end
    sleep(0.1)
  end
  nothing
end

function thread_else(this::Thread)
  println("thread_else")
  timeRef = Ref(0.0)
  i=0
  while true
    if OnTime(2.0, timeRef)
      error(this, i ;title="ELSE")
      i+=1
    end
    sleep(0.5)
  end
  nothing
end

function start_pool()
  pool = thread_pool()

  set(pool, thread_printer, "Printer")
  set(pool, thread_compute, "Compute")
  set(pool, thread_renderer, "Renderer")
  set(pool, thread_sound, "Sound")
  #set(pool, thread_else, "Else")

  #close(pool)

  start(pool) # anything below this line will be paused until threads are closed
end

function start_list()
  list=Function[]
  t = ThreadManager.Thread()

  Base.println(Threads.nthreads())

  Base.push!(list, init(t, thread_printer, "Printer"))
  Base.push!(list, init(t, thread_compute, "Compute"))
  Base.push!(list, init(t, thread_renderer, "Renderer"))
  Base.push!(list, init(t, thread_sound, "Sound"))
  #Base.push!(list, init(t, thread_else, "Else"))

  start(t, list)
end

function main()
  #LoggerManager.log(()->begin
  println("Test Threads")

  start_pool()
  #start_list()

  #end)
end

main()

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
module Examples

export default_delays, mk_cf, ExampleSetting
export Car1DExample, Car2DExamples
export CFName, naive_cf, local_cf, const_u_cf, onevision_cf

using OneVision
using OneVision: ℝ, 𝕋, @kwdef

default_delays = DelayModel(obs = 3, act = 4, com = 6, ΔT = 5)

@enum CFName naive_cf local_cf const_u_cf onevision_cf

@inline function mk_cf(
    name, world_model::WorldDynamics{N}, central, delays,
    loss_model; X, Z, H
) where N
    if name == naive_cf
        NaiveCF(X, Z, N, central, msg_queue_length(delays), delays.ΔT)
    elseif name == local_cf
        LocalCF(central, world_model, delays; X, Z, conpensate_comm = false)
    elseif name == const_u_cf
        LocalCF(central, world_model, delays; X, Z, conpensate_comm = true)
    elseif name == onevision_cf
        OvCF(loss_model, delays; Z, H)
    else
        error("Unexpected CF name: $name")
    end
end

@kwdef mutable struct ExampleSetting
    "The length of the simulation (in seconds)"
    time_end::ℝ
    "physics freqency (in Hz)"
    freq::ℝ
    "Sensor noise strength"
    sensor_noise::ℝ
    "external disturbance strength"
    noise::ℝ
    "Delay model"
    delays::DelayModel
    "Prediction horizon"
    H::𝕋
    "The error ratio between the modeled and true dynamics"
    model_error::ℝ
end

include("Car1DExample.jl")
include("Car2DExamples.jl")

end # module Examples


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

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

# ╔═╡ 86e1ee96-f314-11ea-03f6-0f549b79e7c9
begin
	using Pkg
	Pkg.activate(mktempdir())
end

# ╔═╡ a4937996-f314-11ea-2ff9-615c888afaa8
begin
	Pkg.add([
			"Compose",
			"Colors",
			"PlutoUI",
			])

	using Colors
	using PlutoUI
	using Compose
	using LinearAlgebra
end

# ╔═╡ e6b6760a-f37f-11ea-3ae1-65443ef5a81a
md"_homework 3, version 3_"

# ╔═╡ 85cfbd10-f384-11ea-31dc-b5693630a4c5
md"""

# **Homework 3**: _Structure and Language_
`18.S191`, fall 2020

This notebook contains _built-in, live answer checks_! In some exercises you will see a coloured box, which runs a test case on your code, and provides feedback based on the result. Simply edit the code, run it, and the check runs again.

_For MIT students:_ there will also be some additional (secret) test cases that will be run as part of the grading process, and we will look at your notebook and write comments.

Feel free to ask questions!
"""

# ╔═╡ 33e43c7c-f381-11ea-3abc-c942327456b1
# edit the code below to set your name and kerberos ID (i.e. email without @mit.edu)

student = (name = "Panagiotis Georgakopoulos", kerberos_id = "jazz")

# you might need to wait until all other cells in this notebook have completed running.
# scroll around the page to see what's up

# ╔═╡ ec66314e-f37f-11ea-0af4-31da0584e881
md"""

Submission by: **_$(student.name)_** ($(student.kerberos_id)@mit.edu)
"""

# ╔═╡ 938185ec-f384-11ea-21dc-b56b7469f798
md"_Let's create a package environment:_"

# ╔═╡ b49a21a6-f381-11ea-1a98-7f144c55c9b7
html"""

"""

# ╔═╡ 6f9df800-f92d-11ea-2d49-c1aaabd2d012
md"""
## **Exercise 1:** _Language detection_

In this exercise, we are going to create some super simple _Artificial Intelligence_. Natural language can be quite messy, but hidden in this mess is _structure_, which we are going to look for today.

Let's start with some obvious structure in English text: the set of characters that we write the language in. If we generate random text by sampling _random Unicode characters_, it does not look like English:
"""

# ╔═╡ b61722cc-f98f-11ea-22ae-d755f61f78c3
String(rand(Char, 40))

# ╔═╡ f457ad44-f990-11ea-0e2d-2bb7627716a8
md"""
Instead, let's define an _alphabet_, and only use those letters to sample from. To keep things simple, we ignore punctuation, capitalization, etc, and only use these 27 characters:
"""

# ╔═╡ 4efc051e-f92e-11ea-080e-bde6b8f9295a
alphabet = ['a':'z' ; ' ']   # includes the space character

# ╔═╡ 38d1ace8-f991-11ea-0b5f-ed7bd08edde5
md"""
Let's sample random characters from our alphabet:
"""

# ╔═╡ ddf272c8-f990-11ea-2135-7bf1a6dca0b7
String(rand(alphabet, 40)) |> Text

# ╔═╡ 3cc688d2-f996-11ea-2a6f-0b4c7a5b74c2
md"""
That already looks a lot better than our first attempt! But still, this does not look like English text -- we can do better.

$(html"
")

English words are not well modelled by this random-Latin-characters model. Our first observation is that **some letters are more common than others**. To put this observation into practice, we would like to have the **frequency table** of the Latin alphabet. We could search for it online, but it is actually very simple to calculate ourselves! The only thing we need is a _representative sample_ of English text.

The following samples are from Wikipedia, but feel free to type in your own sample! You can also enter a sample of a different language, if that language can be expressed in the Latin alphabet.

Remember that the $(html"") button on the left of a cell will show or hide the code.

We also include a sample of Spanish, which we'll use later!
"""

# ╔═╡ a094e2ac-f92d-11ea-141a-3566552dd839
md"""
#### Exercise 1.1 - _Data cleaning_
Looking at the sample, we see that it is quite _messy_: it contains punctuation, accented letters and numbers. For our analysis, we are only interested in our 27-character alphabet (i.e. `'a'` through `'z'` plus `' '`). We are going to clean the data using the Julia function `filter`.
"""

# ╔═╡ 27c9a7f4-f996-11ea-1e46-19e3fc840ad9
filter(isodd, [6, 7, 8, 9, -5])

# ╔═╡ f2a4edfa-f996-11ea-1a24-1ba78fd92233
md"""
`filter` takes two arguments: a **function** and a **collection**. The function is applied to each element of the collection, and it returns either `true` or `false`. If the result is `true`, then that element is included in the final collection.

Did you notice something cool? Functions are also just _objects_ in Julia, and you can use them as arguments to other functions! _(Fons thinks this is super cool.)_

$(html"
")

We have written a function `isinalphabet`, which takes a character and returns a boolean:
"""

# ╔═╡ 5c74a052-f92e-11ea-2c5b-0f1a3a14e313
function isinalphabet(character)
	character ∈ alphabet
end

# ╔═╡ dcc4156c-f997-11ea-3e6f-057cd080d9db
isinalphabet('a'), isinalphabet('+')

# ╔═╡ 129fbcfe-f998-11ea-1c96-0fd3ccd2dcf8
md"👉 Use `filter` to extract just the characters from our alphabet out of `messy_sentence_1`:"

# ╔═╡ 3a5ee698-f998-11ea-0452-19b70ed11a1d
messy_sentence_1 = "#wow 2020 ¥500 (blingbling!)"

# ╔═╡ 75694166-f998-11ea-0428-c96e1113e2a0
cleaned_sentence_1 = filter(isinalphabet, messy_sentence_1)

# ╔═╡ 05f0182c-f999-11ea-0a52-3d46c65a049e
md"""
$(html"
")

We are not interested in the case of letters (i.e. `'A'` vs `'a'`), so we want to map these to lower case _before_ we apply our filter. If we don't, all upper case letters would get deleted.
"""

# ╔═╡ 98266882-f998-11ea-3270-4339fb502bc7
md"👉 Use the function `lowercase` to convert `messy_sentence_2` into a lower case string, and then use `filter` to extract only the characters from our alphabet."

# ╔═╡ d3c98450-f998-11ea-3caf-895183af926b
messy_sentence_2 = "Awesome! 😍"

# ╔═╡ d3a4820e-f998-11ea-2a5c-1f37e2a6dd0a
cleaned_sentence_2 = filter(isinalphabet, lowercase.(messy_sentence_2))

# ╔═╡ aad659b8-f998-11ea-153e-3dae9514bfeb
md"""
$(html"
")

Finally, we need to deal with **accents**: simply deleting accented characters from the source text might deform it too much. We can add accented letters to our alphabet, but a simpler solution is to *replace* accented letters with the corresponding unaccented base character. We have written a function `unaccent` that does just that.
"""

# ╔═╡ d236b51e-f997-11ea-0c55-abb11eb35f4d
french_word = "Égalité!"

# ╔═╡ 24860970-fc48-11ea-0009-cddee695772c
import Unicode

# ╔═╡ 734851c6-f92d-11ea-130d-bf2a69e89255
"""
Turn `"áéíóúüñ asdf"` into `"aeiouun asdf"`.
"""
unaccent(str) = Unicode.normalize(str, stripmark=true)

# ╔═╡ d67034d0-f92d-11ea-31c2-f7a38ebb412f
samples = (
	English = """
Although the word forest is commonly used, there is no universally recognised precise definition, with more than 800 definitions of forest used around the world.[4] Although a forest is usually defined by the presence of trees, under many definitions an area completely lacking trees may still be considered a forest if it grew trees in the past, will grow trees in the future,[9] or was legally designated as a forest regardless of vegetation type.[10][11]

The word forest derives from the Old French forest (also forès), denoting "forest, vast expanse covered by trees"; forest was first introduced into English as the word denoting wild land set aside for hunting[14] without the necessity in definition of having trees on the land.[15] Possibly a borrowing, probably via Frankish or Old High German, of the Medieval Latin foresta, denoting "open wood", Carolingian scribes first used foresta in the Capitularies of Charlemagne specifically to denote the royal hunting grounds of the King. The word was not endemic to Romance languages, e. g. native words for forest in the Romance languages derived from the Latin silva, which denoted "forest" and "wood(land)" (confer the English sylva and sylvan); confer the Italian, Spanish, and Portuguese selva; the Romanian silvă; and the Old French selve, and cognates in Romance languages, e. g. the Italian foresta, Spanish and Portuguese floresta, etc., are all ultimately derivations of the French word.
""",
	Spanish =  """
Un bosque es un ecosistema donde la vegetación predominante la constituyen los árboles y matas.1 Estas comunidades de plantas cubren grandes áreas del globo terráqueo y funcionan como hábitats para los animales, moduladores de flujos hidrológicos y conservadores del suelo, constituyendo uno de los aspectos más importantes de la biosfera de la Tierra. Aunque a menudo se han considerado como consumidores de dióxido de carbono atmosférico, los bosques maduros son prácticamente neutros en cuanto al carbono, y son solamente los alterados y los jóvenes los que actúan como dichos consumidores.23 De cualquier manera, los bosques maduros juegan un importante papel en el ciclo global del carbono, como reservorios estables de carbono y su eliminación conlleva un incremento de los niveles de dióxido de carbono atmosférico.

Los bosques pueden hallarse en todas las regiones capaces de mantener el crecimiento de árboles, hasta la línea de árboles, excepto donde la frecuencia de fuego natural es demasiado alta, o donde el ambiente ha sido perjudicado por procesos naturales o por actividades humanas. Los bosques a veces contienen muchas especies de árboles dentro de una pequeña área (como la selva lluviosa tropical y el bosque templado caducifolio), o relativamente pocas especies en áreas grandes (por ejemplo, la taiga y bosques áridos montañosos de coníferas). Los bosques son a menudo hogar de muchos animales y especies de plantas, y la biomasa por área de unidad es alta comparada a otras comunidades de vegetación. La mayor parte de esta biomasa se halla en el subsuelo en los sistemas de raíces y como detritos de plantas parcialmente descompuestos. El componente leñoso de un bosque contiene lignina, cuya descomposición es relativamente lenta comparado con otros materiales orgánicos como la celulosa y otros carbohidratos. Los bosques son áreas naturales y silvestre
""" |> unaccent,
)

# ╔═╡ a56724b6-f9a0-11ea-18f2-991e0382eccf
unaccent(french_word)

# ╔═╡ 8d3bc9ea-f9a1-11ea-1508-8da4b7674629
md"""
$(html"
")

👉 Let's put everything together. Write a function `clean` that takes a string, and returns a _cleaned_ version, where:
- accented letters are replaced by their base characters;
- upper-case letters are converted to lower case;
- it is filtered to only contain characters from `alphabet`
"""

# ╔═╡ 4affa858-f92e-11ea-3ece-258897c37e51
function clean(text)
	# we turn everything to lowercase to keep the number of letters small
	filter(isinalphabet, unaccent(lowercase.(text)))
end

# ╔═╡ e00d521a-f992-11ea-11e0-e9da8255b23b
clean("Crème brûlée est mon plat préféré.")

# ╔═╡ 2680b506-f9a3-11ea-0849-3989de27dd9f
first_sample = clean(first(samples))

# ╔═╡ 571d28d6-f960-11ea-1b2e-d5977ecbbb11
function letter_frequencies(txt)
	ismissing(txt) && return missing
	f = count.(string.(alphabet), txt)
	f ./ sum(f)
end

# ╔═╡ 6a64ab12-f960-11ea-0d92-5b88943cdb1a
sample_freqs = letter_frequencies(first_sample)

# ╔═╡ 603741c2-f9a4-11ea-37ce-1b36ecc83f45
md"""
The result is a 27-element array, with values between `0.0` and `1.0`. These values correspond to the _frequency_ of each letter.

`sample_freqs[i] == 0.0` means that the $i$th letter did not occur in your sample, and
`sample_freqs[i] == 0.1` means that 10% of the letters in the sample are the $i$th letter.

To make it easier to convert between a character from the alphabet and its index, we have the following function:
"""

# ╔═╡ b3de6260-f9a4-11ea-1bae-9153a92c3fe5
index_of_letter(letter) = findfirst(isequal(letter), alphabet)

# ╔═╡ a6c36bd6-f9a4-11ea-1aba-f75cecc90320
index_of_letter('a'), index_of_letter('b'), index_of_letter(' ')

# ╔═╡ 6d3f9dae-f9a5-11ea-3228-d147435e266d
md"""
$(html"
")

👉 Which letters from the alphabet did not occur in the sample?
"""

# ╔═╡ 92bf9fd2-f9a5-11ea-25c7-5966e44db6c6
unused_letters = let
	[l for l ∈ alphabet if sample_freqs[index_of_letter(l)] === 0.0]
end

# ╔═╡ 01215e9a-f9a9-11ea-363b-67392741c8d4
md"""
**Random letters at the correct frequencies:**
"""

# ╔═╡ 8ae13cf0-f9a8-11ea-3919-a735c4ed9e7f
md"""
By considering the _frequencies_ of letters in English, we see that our model is already a lot better!

Our next observation is that **some letter _combinations_ are more common than others**. Our current model thinks that `potato` is just as 'English' as `ooaptt`. In the next section, we will quantify these _transition frequencies_, and use it to improve our model.
"""

# ╔═╡ 343d63c2-fb58-11ea-0cce-efe1afe070c2

# ╔═╡ b5b8dd18-f938-11ea-157b-53b145357fd1
function rand_sample(frequencies)
	x = rand()
	findfirst(z -> z >= x, cumsum(frequencies ./ sum(frequencies)))
end

# ╔═╡ 0e872a6c-f937-11ea-125e-37958713a495
function rand_sample_letter(frequencies)
	alphabet[rand_sample(frequencies)]
end

# ╔═╡ fbb7c04e-f92d-11ea-0b81-0be20da242c8
function transition_counts(cleaned_sample)
	[count(string(a, b), cleaned_sample)
		for a in alphabet,
			b in alphabet]
end

# ╔═╡ 80118bf8-f931-11ea-34f3-b7828113ffd8
normalize_array(x) = x ./ sum(x)

# ╔═╡ 7f4f6ce4-f931-11ea-15a4-b3bec6a7e8b6
transition_frequencies = normalize_array ∘ transition_counts;

# ╔═╡ d40034f6-f9ab-11ea-3f65-7ffd1256ae9d
transition_frequencies(first_sample)

# ╔═╡ 689ed82a-f9ae-11ea-159c-331ff6660a75
md"What we get is a **27 by 27 matrix**. Each entry corresponds to a character pair. The _row_ corresponds to the first character, the _column_ is the second character. Let's visualize this:"

# ╔═╡ 0b67789c-f931-11ea-113c-35e5edafcbbf
md"""
Answer the following questions with respect to the **cleaned English sample text**, which we called `first_sample`. Let's also give the transition matrix a name:
"""

# ╔═╡ 6896fef8-f9af-11ea-0065-816a70ba9670
sample_freq_matrix = transition_frequencies(first_sample);

# ╔═╡ 39152104-fc49-11ea-04dd-bb34e3600f2f
if first_sample === missing
	md"""
	!!! danger "Don't worry!"
	    👆 These errors will disappear automatically once you have completed the earlier exercises!
	"""
end

# ╔═╡ e91c6fd8-f930-11ea-01ac-476bbde79079
md"""👉 What is the frequency of the combination `"th"`?"""

# ╔═╡ 1b4c0c28-f9ab-11ea-03a6-69f69f7f90ed
th_frequency = sample_freq_matrix[index_of_letter('t'), index_of_letter('h')]

# ╔═╡ 1f94e0a2-f9ab-11ea-1347-7dd906ebb09d
md"""👉 What about `"ht"`?"""

# ╔═╡ 41b2df7c-f931-11ea-112e-ede3b16f357a
ht_frequency = sample_freq_matrix[index_of_letter('h'), index_of_letter('t')]

# ╔═╡ 1dd1e2f4-f930-11ea-312c-5ff9e109c7f6
md"""
👉 Which le**tt**ers appeared double in our sample?
"""

# ╔═╡ 65c92cac-f930-11ea-20b1-6b8f45b3f262
double_letters = [l for l ∈ 'a':'z' if sample_freq_matrix[index_of_letter(l), index_of_letter(l)] > 0]

# ╔═╡ 4582ebf4-f930-11ea-03b2-bf4da1a8f8df
md"""
👉 Which letter is most likely to follow a **W**?
"""

# ╔═╡ 6d8f47d2-2108-11eb-105e-63c5f7d74c6e
getfreq = (p, n) -> sample_freq_matrix[index_of_letter(p), index_of_letter(n)]

# ╔═╡ b5fcd572-2108-11eb-3581-b18f259b42f6
index_of_max = a -> begin
	m = max(a...)
	return findfirst(i -> i == m, a)
end

# ╔═╡ 7898b76a-f930-11ea-2b7e-8126ec2b8ffd
most_likely_to_follow_w = alphabet[index_of_max([getfreq('w', l) for l ∈ alphabet])]

# ╔═╡ 458cd100-f930-11ea-24b8-41a49f6596a0
md"""
👉 Which letter is most likely to precede a **W**?
"""

# ╔═╡ bc401bee-f931-11ea-09cc-c5efe2f11194
most_likely_to_precede_w = alphabet[index_of_max([getfreq(l, 'w') for l ∈ alphabet])]

# ╔═╡ 45c20988-f930-11ea-1d12-b782d2c01c11
md"""
👉 What is the sum of each row? What is the sum of each column? How can we interpret these values?"
"""

# ╔═╡ 432276c6-2109-11eb-04e8-0bb0716eb660
sum(sample_freq_matrix, dims=2)[:, 1]

# ╔═╡ cc62929e-f9af-11ea-06b9-439ac08dcb52
row_col_answer = md"""

"""

# ╔═╡ 2f8dedfc-fb98-11ea-23d7-2159bdb6a299
md"""
We can use the measured transition frequencies to generate text in a way that it has **the same transition frequencies** as our original sample. Our generated text is starting to look like real language!
"""

# ╔═╡ b7446f34-f9b1-11ea-0f39-a3c17ba740e5
@bind ex23_sample Select([v => String(k) for (k,v) in zip(fieldnames(typeof(samples)), samples)])

# ╔═╡ 4f97b572-f9b0-11ea-0a99-87af0797bf28
md"""
**Random letters from the alphabet:**
"""

# ╔═╡ 4e8d327e-f9b0-11ea-3f16-c178d96d07d9
md"""
**Random letters at the correct frequencies:**
"""

# ╔═╡ d83f8bbc-f9af-11ea-2392-c90e28e96c65
md"""
**Random letters at the correct transition frequencies:**
"""

# ╔═╡ 0e465160-f937-11ea-0ebb-b7e02d71e8a8
function sample_text(A, n)

	first_index = rand_sample(vec(sum(A, dims=1)))

	indices = reduce(1:n; init=[first_index]) do word, _
		prev = last(word)
		freq = normalize_array(A[prev, :])
		next = rand_sample(freq)
		[word..., next]
	end

	String(alphabet[indices])
end

# ╔═╡ 141af892-f933-11ea-1e5f-154167642809
md"""
It looks like we have a decent language model, in the sense that it understands _transition frequencies_ in the language. In the demo above, try switching the language between $(join(string.(fieldnames(typeof(samples))), " and ")) -- the generated text clearly looks more like one or the other, demonstrating that the model can capture differences between the two languages. What's remarkable is that our "training data" was just a single paragraph per language.

In this exercise, we will use our model to write a **classifier**: a program that automatically classifies a text as either $(join(string.(fieldnames(typeof(samples))), " or ")).

This is not a difficult task - you can get dictionaries for both languages, and count matches - but we are doing something much more cool: we only use a single paragraph of each language, and we use a _language model_ as classifier.
"""

# ╔═╡ 7eed9dde-f931-11ea-38b0-db6bfcc1b558
html"
Mystery sample

Enter some text here -- we will detect in which language it is written!
" # dont delete me

# ╔═╡ 7e3282e2-f931-11ea-272f-d90779264456
@bind mystery_sample TextField((70,8), default="""
Small boats are typically found on inland waterways such as rivers and lakes, or in protected coastal areas. However, some boats, such as the whaleboat, were intended for use in an offshore environment. In modern naval terms, a boat is a vessel small enough to be carried aboard a ship. Anomalous definitions exist, as lake freighters 1,000 feet (300 m) long on the Great Lakes are called "boats".
""")

# ╔═╡ 7df55e6c-f931-11ea-33b8-fdc3be0b6cfa
mystery_sample

# ╔═╡ 292e0384-fb57-11ea-0238-0fbe416fc976
md"""
Let's compute the transition frequencies of our mystery sample! Type some text in the box below, and check whether the frequency matrix updates.
"""

# ╔═╡ 7dabee08-f931-11ea-0cb2-c7d5afd21551
transition_frequencies(mystery_sample)

# ╔═╡ 3736a094-fb57-11ea-1d39-e551aae62b1d
md"""
Our model will **compare the transition frequencies of our mystery sample** to those of our two language sample. The closest match will be our detected language.

The only question left is: How do we compare two matrices? When two matrices are almost equal, but not exactly, we want to quantify their _distance_.

👉 Write a function called `matrix_distance` which takes 2 matrices of the same size and finds the distance between them by:

1. Subtracting corresponding elements
2. Finding the absolute value of the difference
3. Summing the differences
"""

# ╔═╡ 13c89272-f934-11ea-07fe-91b5d56dedf8
function matrix_distance(A, B)
	sum(abs.(A .- B))
end

# ╔═╡ 7d60f056-f931-11ea-39ae-5fa18a955a77
distances = map(samples) do sample
	matrix_distance(transition_frequencies(mystery_sample), transition_frequencies(sample))
end

# ╔═╡ 7d1439e6-f931-11ea-2dab-41c66a779262
try
	@assert !ismissing(distances.English)
	"""
It looks like this text is **$(argmin(distances))**!
""" |> HTML
catch
end

# ╔═╡ 8c7606f0-fb93-11ea-0c9c-45364892cbb8
md"""
We have written a cell that selects the language with the _smallest distance_ to the mystery language. Make sure sure that `matrix_distance` is working correctly, and [scroll up](#mystery-detect) to the mystery text to see it in action!

#### Further reading
It turns out that the SVD of the transition matrix can mysteriously group the alphabet into vowels and consonants, without any extra information. See [this paper](http://languagelog.ldc.upenn.edu/myl/Moler1983.pdf) if you want to try it yourself! We found that removing the space from `alphabet` (to match the paper) gave better results.
"""

# ╔═╡ 82e0df62-fb54-11ea-3fff-b16c87a7d45b
md"""
## **Exercise 2** - _Language generation_

Our model from Exercise 1 has the property that it can easily be 'reversed' to _generate_ text. While this is useful to demonstrate its structure, the produced text is mostly meaningless: it fails to model words, let alone sentence structure.

To take our model one step further, we are going to _generalize_ what we have done so far. Instead of looking at _letter combinations_, we will model _word combinations_.  And instead of analyzing the frequencies of bigrams (combinations of two letters), we are going to analyze _$n$-grams_.

#### Dataset
This also means that we are going to need a larger dataset to train our model on: the number of English words (and their combinations) is much higher than the number of letters.

We will train our model on the novel [_Emma_ (1815), by Jane Austen](https://en.wikipedia.org/wiki/Emma_(novel)). This work is in the public domain, which means that we can download the whole book as a text file from `archive.org`. We've done the process of downloading and cleaning already, and we have split the text into word and punctuation tokens.
"""

# ╔═╡ b7601048-fb57-11ea-0754-97dc4e0623a1
emma = let
	raw_text = read(download("https://ia800303.us.archive.org/24/items/EmmaJaneAusten_753/emma_pdf_djvu.txt"), String)

	first_words = "Emma Woodhouse"
	last_words = "THE END"
	start_index = findfirst(first_words, raw_text)[1]
	stop_index = findlast(last_words, raw_text)[end]

	raw_text[start_index:stop_index]
end;

# ╔═╡ cc42de82-fb5a-11ea-3614-25ef961729ab
function splitwords(text)
	# clean up whitespace
	cleantext = replace(text, r"\s+" => " ")

	# split on whitespace or other word boundaries
	tokens = split(cleantext, r"(\s|\b)")
end

# ╔═╡ d66fe2b2-fb5a-11ea-280f-cfb12b8296ac
emma_words = splitwords(emma)

# ╔═╡ 4ca8e04a-fb75-11ea-08cc-2fdef5b31944
forest_words = splitwords(samples.English)

# ╔═╡ 6f613cd2-fb5b-11ea-1669-cbd355677649
md"""
#### Exercise 2.1 - _bigrams revisited_

The goal of the upcoming exercises is to **generalize** what we have done in Exercise 1. To keep things simple, we _split up our problem_ into smaller problems. (The solution to any computational problem.)

First, here is a function that takes an array, and returns the array of all **neighbour pairs** from the original. For example,

```julia
bigrams([1, 2, 3, 42])
```
gives

```julia
[[1,2], [2,3], [3,42]]
```

(We used integers as "words" in this example, but our function works with any type of word.)
"""

# ╔═╡ 91e87974-fb78-11ea-3ce4-5f64e506b9d2
function bigrams(words)
	map(1:length(words)-1) do i
		words[i:i+1]
	end
end

# ╔═╡ 9f98e00e-fb78-11ea-0f6c-01206e7221d6
bigrams([1, 2, 3, 42])

# ╔═╡ d7d8cd0c-fb6a-11ea-12bf-2d1448b38162
md"""
👉 Next, it's your turn to write a more general function `ngrams` that takes an array and a number $n$, and returns all **subsequences of length $n$**. For example:

```julia
ngrams([1, 2, 3, 42], 3)
```
should give

```julia
[[1,2,3], [2,3,42]]
```

and

```julia
ngrams([1, 2, 3, 42], 2) == bigrams([1, 2, 3, 42])
```
"""

# ╔═╡ 7be98e04-fb6b-11ea-111d-51c48f39a4e9
function ngrams(words, n)
	map(1:length(words) - n + 1) do i
		words[i: (i + n - 1)]
	end
end

# ╔═╡ c560ec86-211c-11eb-06e9-b3cc1c80e46d

# ╔═╡ 052f822c-fb7b-11ea-382f-af4d6c2b4fdb
ngrams([1, 2, 3, 42], 3)

# ╔═╡ 067f33fc-fb7b-11ea-352e-956c8727c79f
ngrams(forest_words, 4)

# ╔═╡ 7b10f074-fb7c-11ea-20f0-034ddff41bc3
md"""
If you are stuck, you can write `ngrams(words, n) = bigrams(words)` (ignoring the true value of $n$), and continue with the other exercises.

#### Exercise 2.2 - _frequency matrix revisisted_
In Exercise 1, we use a 2D array to store the bigram frequencies, where each column or row corresponds to a character from the alphabet. If we use trigrams, we could store the frequencies in a 3D array, and so on.

However, when counting words instead of letters, we run into a problem. A 3D array with one row, column and layer per word has too many elements to store on our computer.
"""

# ╔═╡ 24ae5da0-fb7e-11ea-3480-8bb7b649abd5
md"""
_Emma_ consists of $(
	length(Set(emma_words))
) unique words. This means that there are $(
	Int(floor(length(Set(emma_words))^3 / 10^9))
) billion possible trigrams - that's too much!
"""

# ╔═╡ 47836744-fb7e-11ea-2305-3fa5819dc154
md"""
$(html"
")

Although the frequency array would be very large, most entries are zero. For example, _"Emma"_ is a common word, but _"Emma Emma Emma"_ does not occur in the novel. This _sparsity_ of non-zero entries can be used to **store the same information in a more efficient structure**.

Julia's [`SparseArrays.jl` package](https://docs.julialang.org/en/v1/stdlib/SparseArrays/index.html) might sound like a logical choice, but these arrays only support 1D and 2D types, and we also want to directly index using strings, not just integers. So instead, we will use `Dict`, the dictionary type.
"""

# ╔═╡ df4fc31c-fb81-11ea-37b3-db282b36f5ef
healthy = Dict("fruits" => ["🍎", "🍊"], "vegetables" => ["🌽", "🎃", "🍕"])

# ╔═╡ c83b1770-fb82-11ea-20a6-3d3a09606c62
healthy["fruits"]

# ╔═╡ 52970ac4-fb82-11ea-3040-8bd0590348d2
md"""
(Did you notice something funny? The dictionary is _unordered_, this is why the entries were printed in reverse from the definition.)

You can dynamically add or change values of a `Dict` by assigning to `my_dict[key]`. You can check whether a key already exists using `haskey(my_dict, key)`.

👉 Use these two techniques to write a function `word_counts` that takes an array of words, and returns a `Dict` with entries `word => number_of_occurences`.

For example:
```julia
word_counts(["to", "be", "or", "not", "to", "be"])
```
should return
```julia
Dict(
	"to" => 2,
	"be" => 2,
	"or" => 1,
	"not" => 1
)
```
"""

# ╔═╡ 8ce3b312-fb82-11ea-200c-8d5b12f03eea
function word_counts(words::Vector)
	counts = Dict()
	for word ∈ words
		if haskey(counts, word)
			counts[word] = counts[word] + 1
		else
			counts[word] = 1
		end
	end
	return counts
end

# ╔═╡ a2214e50-fb83-11ea-3580-210f12d44182
word_counts(["to", "be", "or", "not", "to", "be"])

# ╔═╡ 808abf6e-fb84-11ea-0785-2fc3f1c4a09f
md"""
How many times does `"Emma"` occur in the book?
"""

# ╔═╡ 953363dc-fb84-11ea-1128-ebdfaf5160ee
emma_count = word_counts(emma_words)["Emma"]

# ╔═╡ 294b6f50-fb84-11ea-1382-03e9ab029a2d
md"""
Great! Let's get back to our ngrams. For the purpose of generating text, we are going to store a _completions cache_. This is a dictionary where the keys are $(n-1)$-grams, and the values are all found words that complete it to an $n$-gram. Let's look at an example:

```julia
let
	trigrams = ngrams(split("to be or not to be that is the question", " "), 3)
	cache = completions_cache(trigrams)
	cache == Dict(
		["to", "be"] => ["or", "that"],
		["be", "or"] => ["not"],
		["or", "not"] => ["to"],
		...
	)
end
```

So for trigrams, our keys are the first $2$ words of each trigram, and the values are arrays containing every third word of those trigrams.

If the same ngram occurs multiple times (e.g. "said Emma laughing"), then the last word ("laughing") should also be stored multiple times. This will allow us to generate trigrams with the same frequencies as the original text.

👉 Write the function `completions_cache`, which takes an array of ngrams (i.e. an array of arrays of words, like the result of your `ngram` function), and returns a dictionary like described above.
"""

# ╔═╡ b726f824-fb5e-11ea-328e-03a30544037f
function completions_cache(grams)
	cache = Dict()
	for gram ∈ grams
    	L = length(gram)
		key = gram[1:L-1]
		val = gram[L]
		if haskey(cache, key)
			push!(cache[key], val)
		else
			cache[key] = [val]
		end
	end
	cache
end

# ╔═╡ 18355314-fb86-11ea-0738-3544e2e3e816
let
	trigrams = ngrams(split("to be or not to be or is the question", " "), 3)
	completions_cache(trigrams)
end

# ╔═╡ 3d105742-fb8d-11ea-09b0-cd2e77efd15c
md"""
#### Exercise 2.4 - _write a novel_

We have everything we need to generate our own novel! The final step is to sample random ngrams, in a way that each next ngram overlaps with the previous one. We've done this in the function `generate_from_ngrams` below - feel free to look through the code, or to implement your own version.
"""

# ╔═╡ a72fcf5a-fb62-11ea-1dcc-11451d23c085
"""
	generate_from_ngrams(grams, num_words)

Given an array of ngrams (i.e. an array of arrays of words), generate a sequence of `num_words` words by sampling random ngrams.
"""
function generate_from_ngrams(grams, num_words)
	n = length(first(grams))
	cache = completions_cache(grams)

	# we need to start the sequence with at least n-1 words.
	# a simple way to do so is to pick a random ngram!
	sequence = [rand(grams)...]

	# we iteratively add one more word at a time
	for i ∈ n+1:num_words
		# the previous n-1 words
		tail = sequence[end-(n-2):end]

		# possible next words
		completions = cache[tail]

		choice = rand(completions)
		push!(sequence, choice)
	end
	sequence
end

# ╔═╡ f83991c0-fb7c-11ea-0e6f-1f80709d00c1
"Compute the ngrams of an array of words, but add the first n-1 at the end, to ensure that every ngram ends in the the beginning of another ngram."
function ngrams_circular(words, n)
	ngrams([words..., words[1:n-1]...], n)
end

# ╔═╡ abe2b862-fb69-11ea-08d9-ebd4ba3437d5
completions_cache(ngrams_circular(forest_words, 3))

# ╔═╡ 4b27a89a-fb8d-11ea-010b-671eba69364e
"""
	generate(source_text::AbstractString, num_token; n=3, use_words=true)

Given a source text, generate a `String` that "looks like" the original text by satisfying the same ngram frequency distribution as the original.
"""
function generate(source_text::AbstractString, s; n=3, use_words=true)
	preprocess = if use_words
		splitwords
	else
		collect
	end

	words = preprocess(source_text)
	if length(words) < n
		""
	else
		grams = ngrams_circular(words, n)
		result = generate_from_ngrams(grams, s)
		if use_words
			join(result, " ")
		else
			String(result)
		end
	end
end

# ╔═╡ d7b7a14a-fb90-11ea-3e2b-2fd8f379b4d8
md"
#### Interactive demo

Enter your own text in the box below, and use that as training data to generate anything!
"

# ╔═╡ 1939dbea-fb63-11ea-0bc2-2d06b2d4b26c
@bind generate_demo_sample TextField((50,5), default=samples.English)

# ╔═╡ 70169682-fb8c-11ea-27c0-2dad2ff3080f
md"""Using $(@bind generate_sample_n_letters NumberField(1:5))grams for characters"""

# ╔═╡ 402562b0-fb63-11ea-0769-375572cc47a8
md"""Using $(@bind generate_sample_n_words NumberField(1:5))grams for words"""

# ╔═╡ 2521bac8-fb8f-11ea-04a4-0b077d77529e
md"""
### Automatic Jane Austen

Uncomment the cell below to generate some Jane Austen text:
"""

# ╔═╡ cc07f576-fbf3-11ea-2c6f-0be63b9356fc
if student.name == "Jazzy Doe"
	md"""
	!!! danger "Before you submit"
	    Remember to fill in your **name** and **Kerberos ID** at the top of this notebook.
	"""
end

# ╔═╡ 6b4d6584-f3be-11ea-131d-e5bdefcc791b
md"## Function library

Just some helper functions used in the notebook."

# ╔═╡ 54b1e236-fb53-11ea-3769-b382ef8b25d6
function Quote(text::AbstractString)
	text |> Markdown.Paragraph |> Markdown.BlockQuote |> Markdown.MD
end

# ╔═╡ b3dad856-f9a7-11ea-1552-f7435f1cb605
String(rand(alphabet, 400)) |> Quote

# ╔═╡ be55507c-f9a7-11ea-189c-4ffe8377212e
if sample_freqs !== missing
	String([rand_sample_letter(sample_freqs) for _ in 1:400]) |> Quote
end

# ╔═╡ 46c905d8-f9b0-11ea-36ed-0515e8ed2621
String(rand(alphabet, 400)) |> Quote

# ╔═╡ 489b03d4-f9b0-11ea-1de0-11d4fe4e7c69
String([rand_sample_letter(letter_frequencies(ex23_sample)) for _ in 1:400]) |> Quote

# ╔═╡ fd202410-f936-11ea-1ad6-b3629556b3e0
sample_text(transition_frequencies(clean(ex23_sample)), 400) |> Quote

# ╔═╡ b5dff8b8-fb6c-11ea-10fc-37d2a9adae8c
generate(
	generate_demo_sample, 400;
	n=generate_sample_n_letters,
	use_words=false
) |> Quote

# ╔═╡ ee8c5808-fb5f-11ea-19a1-3d58217f34dc
generate(
	generate_demo_sample, 100;
	n=generate_sample_n_words,
	use_words=true
) |> Quote

# ╔═╡ 49b69dc2-fb8f-11ea-39af-030b5c5053c3
generate(emma, 100; n=4) |> Quote

# ╔═╡ ddef9c94-fb96-11ea-1f17-f173a4ff4d89
function compimg(img, labels=[c*d for c in replace(alphabet, ' ' => "_"), d in replace(alphabet, ' ' => "_")])
	xmax, ymax = size(img)
	xmin, ymin = 0, 0
	arr = [(j-1, i-1) for i=1:ymax, j=1:xmax]

	compose(context(units=UnitBox(xmin, ymin, xmax, ymax)),
		fill(vec(img)),
		compose(context(),
			fill("white"), font("monospace"),
			text(first.(arr) .+ .1, last.(arr) .+ 0.6, labels)),
		rectangle(
			first.(arr),
			last.(arr),
			fill(1.0, length(arr)),
			fill(1.0, length(arr))))
end

# ╔═╡ b7803a28-fb96-11ea-3e30-d98eb322d19a
function show_pair_frequencies(A)
	imshow = let
		to_rgb(x) = RGB(0.36x, 0.82x, 0.8x)
		to_rgb.(A ./ maximum(abs.(A)))
	end
	compimg(imshow)
end

# ╔═╡ ace3dc76-f9ae-11ea-2bee-3d0bfa57cfbc
show_pair_frequencies(transition_frequencies(first_sample))

# ╔═╡ ffc17f40-f380-11ea-30ee-0fe8563c0eb1
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]))

# ╔═╡ 7df7ab82-f9ad-11ea-2243-21685d660d71
hint(md"You can answer this question without writing any code: have a look at the values of `sample_freqs`.")

# ╔═╡ e467c1c6-fbf2-11ea-0d20-f5798237c0a6
hint(md"Start out with the same code as `bigrams`, and use the Julia documentation to learn how it works. How can we generalize the `bigram` function into the `ngram` function? It might help to do this on paper first.")

# ╔═╡ ffc40ab2-f380-11ea-2136-63542ff0f386
almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]))

# ╔═╡ ffceaed6-f380-11ea-3c63-8132d270b83f
still_missing(text=md"Replace `missing` with your answer.") = Markdown.MD(Markdown.Admonition("warning", "Here we go!", [text]))

# ╔═╡ ffde44ae-f380-11ea-29fb-2dfcc9cda8b4
keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]))

# ╔═╡ ffe326e0-f380-11ea-3619-61dd0592d409
yays = [md"Fantastic!", md"Splendid!", md"Great!", md"Yay ❤", md"Great! 🎉", md"Well done!", md"Keep it up!", md"Good job!", md"Awesome!", md"You got the right answer!", md"Let's move on to the next section."]

# ╔═╡ fff5aedc-f380-11ea-2a08-99c230f8fa32
correct(text=rand(yays)) = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]))

# ╔═╡ 954fc466-fb7b-11ea-2724-1f938c6b93c6
let
	output = ngrams([1, 2, 3, 42], 2)

	if output isa Missing
		still_missing()
	elseif !(output isa Vector{<:Vector})
		keep_working(md"Make sure that `ngrams` returns an array of arrays.")
	elseif output == [[1,2], [2,3], [3,42]]
		if ngrams([1,2,3], 1) == [[1],[2],[3]]
			if ngrams([1,2,3], 3) == [[1,2,3]]
				if ngrams(["a"],1) == [["a"]]
					correct()
				else
					keep_working(md"`ngrams` should work with any type, not just integers!")
				end
			else
				keep_working(md"`ngrams(x, 3)` did not give a correct result.")
			end
		else
			keep_working(md"`ngrams(x, 1)` did not give a correct result.")
		end
	else
		keep_working(md"`ngrams(x, 2)` did not give the correct bigrams. Start out with the same code as `bigrams`.")
	end
end

# ╔═╡ a9ffff9a-fb83-11ea-1efd-2fc15538e52f
let
	output = word_counts(["to", "be", "or", "not", "to", "be"])

	if output === nothing
		keep_working(md"Did you forget to write `return`?")
	elseif output == Dict()
		still_missing(md"Write your function `word_counts` above.")
	elseif !(output isa Dict)
		keep_working(md"Make sure that `word_counts` returns a `Dict`.")
	elseif output == Dict("to" => 2, "be" => 2, "or" => 1, "not" => 1)
		correct()
	else
		keep_working()
	end
end

# ╔═╡ 00026442-f381-11ea-2b41-bde1fff66011
not_defined(variable_name) = Markdown.MD(Markdown.Admonition("danger", "Oopsie!", [md"Make sure that you define a variable called **$(Markdown.Code(string(variable_name)))**"]))

# ╔═╡ 6fe693c8-f9a1-11ea-1983-f159131880e9
if !@isdefined(messy_sentence_1)
	not_defined(:messy_sentence_1)
elseif !@isdefined(cleaned_sentence_1)
	not_defined(:cleaned_sentence_1)
else
	if cleaned_sentence_1 isa Missing
		still_missing()
	elseif cleaned_sentence_1 isa Vector{Char}
		keep_working(md"Use `String(x)` to turn an array of characters `x` into a `String`.")
	elseif cleaned_sentence_1 == filter(isinalphabet, messy_sentence_1)
		correct()
	else
		keep_working()
	end
end

# ╔═╡ cee0f984-f9a0-11ea-2c3c-53fe26156ea4
if !@isdefined(messy_sentence_2)
	not_defined(:messy_sentence_2)
elseif !@isdefined(cleaned_sentence_2)
	not_defined(:cleaned_sentence_2)
else
	if cleaned_sentence_2 isa Missing
		still_missing()
	elseif cleaned_sentence_2 isa Vector{Char}
		keep_working(md"Use `String(x)` to turn an array of characters `x` into a `String`.")
	elseif cleaned_sentence_2 == filter(isinalphabet, lowercase(messy_sentence_2))
		correct()
	else
		keep_working()
	end
end

# ╔═╡ ddfb1e1c-f9a1-11ea-3625-f1170272e96a
if !@isdefined(clean)
	not_defined(:clean)
else
	let
		input = "Aè !!!  x1"
		output = clean(input)

		if output isa Missing
			still_missing()
		elseif output isa Vector{Char}
			keep_working(md"Use `String(x)` to turn an array of characters `x` into a `String`.")
		elseif output == "ae   x"
			correct()
		else
			keep_working()
		end
	end
end

# ╔═╡ 95b81778-f9a5-11ea-3f51-019430bc8fa8
if !@isdefined(unused_letters)
	not_defined(:unused_letters)
else
	if sample_freqs === missing
		md"""
		!!! warning "Oopsie!"
		    You need to complete the previous exercises first.
		"""
	elseif unused_letters isa Missing
		still_missing()
	elseif unused_letters isa String
		keep_working(md"Use `collect` to turn a string into an array of characters.")
	elseif Set(index_of_letter.(unused_letters)) == Set(findall(isequal(0.0), sample_freqs))
		correct()
	else
		keep_working()
	end
end

# ╔═╡ 489fe282-f931-11ea-3dcb-35d4f2ac8b40
if !@isdefined(th_frequency)
	not_defined(:th_frequency)
elseif !@isdefined(ht_frequency)
	not_defined(:ht_frequency)
else
	if th_frequency isa Missing  || ht_frequency isa Missing
		still_missing()
	elseif th_frequency < ht_frequency
		keep_working(md"Looks like your answers should be flipped. Which combination is more frequent in English?")
	elseif th_frequency == sample_freq_matrix[index_of_letter('t'), index_of_letter('h')] && ht_frequency == sample_freq_matrix[index_of_letter('h'), index_of_letter('t')]
		correct()
	else
		keep_working()
	end
end

# ╔═╡ 671525cc-f930-11ea-0e71-df9d4aae1c05
if !@isdefined(double_letters)
	not_defined(:double_letters)
end

# ╔═╡ a5fbba46-f931-11ea-33e1-054be53d986c
if !@isdefined(most_likely_to_follow_w)
	not_defined(:most_likely_to_follow_w)
end

# ╔═╡ ba695f6a-f931-11ea-0fbb-c3ef1374270e
if !@isdefined(most_likely_to_precede_w)
	not_defined(:most_likely_to_precede_w)
end

# ╔═╡ b09f5512-fb58-11ea-2527-31bea4cee823
if !@isdefined(matrix_distance)
	not_defined(:matrix_distance)
else
	try
	let
		A = rand(Float64, (5,4))
		B = rand(Float64, (5,4))

		output = matrix_distance(A,B)

		if output isa Missing
			still_missing()
		elseif !(output isa Number)
			keep_working(md"Make sure that `matrix_distance` returns a nunmber.")
		elseif output == 0.0
			keep_working(md"Two different matrices should have non-zero distance.")
		else
			if matrix_distance(A,B) < 0 || matrix_distance(B,A) < 0
				keep_working(md"The distance between two matrices should always be positive.")
			elseif matrix_distance(A,A) != 0
				almost(md"The distance between two identical matrices should be zero.")
			elseif matrix_distance([1 -1], [0 0]) == 0.0
				almost(md"`matrix_distance([1 -1], [0 0])` should not be zero.")
			else
				correct()
			end
		end
	end
	catch
		keep_working(md"The function errored.")
	end
end

# ╔═╡ 00115b6e-f381-11ea-0bc6-61ca119cb628
bigbreak = html"

";

# ╔═╡ c086bd1e-f384-11ea-3b26-2da9e24360ca
bigbreak

# ╔═╡ eaa8c79e-f9a2-11ea-323f-8bb2bd36e11c
md"""
$(bigbreak)
#### Exercise 1.2 - _Letter frequencies_

We are going to count the _frequency_ of each letter in this sample, after applying your `clean` function. Can you guess which character is most frequent?
"""

# ╔═╡ dcffd7d2-f9a6-11ea-2230-b1afaecfdd54
md"""
$(bigbreak)
Now that we know the frequencies of letters in English, we can generate random text that already looks closer to English!

**Random letters from the alphabet:**
"""

# ╔═╡ 77623f3e-f9a9-11ea-2f46-ff07bd27cd5f
md"""
$(bigbreak)
#### Exercise 1.3 - _Transition frequencies_
In the previous exercise we computed the frequency of each letter in the sample by _counting_ their occurences, and then dividing by the total number of counts.

In this exercise, we are going to count _letter transitions_, such as `aa`, `as`, `rt`, `yy`. Two letters might both be common, like `a` and `e`, but their combination, `ae`, is uncommon in English.

To quantify this observation, we will do the same as in our last exercise: we count occurences in a _sample text_, to create the **transition frequency matrix**.
"""

# ╔═╡ d3d7bd9c-f9af-11ea-1570-75856615eb5d
bigbreak

# ╔═╡ 6718d26c-f9b0-11ea-1f5a-0f22f7ddffe9
md"""
$(bigbreak)

#### Exercise 1.4 - _Language detection_
"""

# ╔═╡ 568f0d3a-fb54-11ea-0f77-171718ef12a5
bigbreak

# ╔═╡ 7f341c4e-fb54-11ea-1919-d5421d7a2c75
bigbreak

# ╔═╡ Cell order:
# ╟─e6b6760a-f37f-11ea-3ae1-65443ef5a81a
# ╟─ec66314e-f37f-11ea-0af4-31da0584e881
# ╟─85cfbd10-f384-11ea-31dc-b5693630a4c5
# ╠═33e43c7c-f381-11ea-3abc-c942327456b1
# ╟─938185ec-f384-11ea-21dc-b56b7469f798
# ╠═86e1ee96-f314-11ea-03f6-0f549b79e7c9
# ╠═a4937996-f314-11ea-2ff9-615c888afaa8
# ╟─b49a21a6-f381-11ea-1a98-7f144c55c9b7
# ╟─c086bd1e-f384-11ea-3b26-2da9e24360ca
# ╟─6f9df800-f92d-11ea-2d49-c1aaabd2d012
# ╠═b61722cc-f98f-11ea-22ae-d755f61f78c3
# ╟─f457ad44-f990-11ea-0e2d-2bb7627716a8
# ╠═4efc051e-f92e-11ea-080e-bde6b8f9295a
# ╟─38d1ace8-f991-11ea-0b5f-ed7bd08edde5
# ╠═ddf272c8-f990-11ea-2135-7bf1a6dca0b7
# ╟─3cc688d2-f996-11ea-2a6f-0b4c7a5b74c2
# ╟─d67034d0-f92d-11ea-31c2-f7a38ebb412f
# ╟─a094e2ac-f92d-11ea-141a-3566552dd839
# ╠═27c9a7f4-f996-11ea-1e46-19e3fc840ad9
# ╟─f2a4edfa-f996-11ea-1a24-1ba78fd92233
# ╟─5c74a052-f92e-11ea-2c5b-0f1a3a14e313
# ╠═dcc4156c-f997-11ea-3e6f-057cd080d9db
# ╟─129fbcfe-f998-11ea-1c96-0fd3ccd2dcf8
# ╠═3a5ee698-f998-11ea-0452-19b70ed11a1d
# ╠═75694166-f998-11ea-0428-c96e1113e2a0
# ╟─6fe693c8-f9a1-11ea-1983-f159131880e9
# ╟─05f0182c-f999-11ea-0a52-3d46c65a049e
# ╟─98266882-f998-11ea-3270-4339fb502bc7
# ╠═d3c98450-f998-11ea-3caf-895183af926b
# ╠═d3a4820e-f998-11ea-2a5c-1f37e2a6dd0a
# ╟─cee0f984-f9a0-11ea-2c3c-53fe26156ea4
# ╟─aad659b8-f998-11ea-153e-3dae9514bfeb
# ╠═d236b51e-f997-11ea-0c55-abb11eb35f4d
# ╠═a56724b6-f9a0-11ea-18f2-991e0382eccf
# ╟─24860970-fc48-11ea-0009-cddee695772c
# ╟─734851c6-f92d-11ea-130d-bf2a69e89255
# ╟─8d3bc9ea-f9a1-11ea-1508-8da4b7674629
# ╠═4affa858-f92e-11ea-3ece-258897c37e51
# ╠═e00d521a-f992-11ea-11e0-e9da8255b23b
# ╟─ddfb1e1c-f9a1-11ea-3625-f1170272e96a
# ╟─eaa8c79e-f9a2-11ea-323f-8bb2bd36e11c
# ╠═2680b506-f9a3-11ea-0849-3989de27dd9f
# ╟─571d28d6-f960-11ea-1b2e-d5977ecbbb11
# ╠═6a64ab12-f960-11ea-0d92-5b88943cdb1a
# ╟─603741c2-f9a4-11ea-37ce-1b36ecc83f45
# ╟─b3de6260-f9a4-11ea-1bae-9153a92c3fe5
# ╠═a6c36bd6-f9a4-11ea-1aba-f75cecc90320
# ╟─6d3f9dae-f9a5-11ea-3228-d147435e266d
# ╠═92bf9fd2-f9a5-11ea-25c7-5966e44db6c6
# ╟─95b81778-f9a5-11ea-3f51-019430bc8fa8
# ╟─7df7ab82-f9ad-11ea-2243-21685d660d71
# ╟─dcffd7d2-f9a6-11ea-2230-b1afaecfdd54
# ╟─b3dad856-f9a7-11ea-1552-f7435f1cb605
# ╟─01215e9a-f9a9-11ea-363b-67392741c8d4
# ╟─be55507c-f9a7-11ea-189c-4ffe8377212e
# ╟─8ae13cf0-f9a8-11ea-3919-a735c4ed9e7f
# ╟─343d63c2-fb58-11ea-0cce-efe1afe070c2
# ╟─b5b8dd18-f938-11ea-157b-53b145357fd1
# ╟─0e872a6c-f937-11ea-125e-37958713a495
# ╟─77623f3e-f9a9-11ea-2f46-ff07bd27cd5f
# ╠═fbb7c04e-f92d-11ea-0b81-0be20da242c8
# ╠═80118bf8-f931-11ea-34f3-b7828113ffd8
# ╠═7f4f6ce4-f931-11ea-15a4-b3bec6a7e8b6
# ╠═d40034f6-f9ab-11ea-3f65-7ffd1256ae9d
# ╟─689ed82a-f9ae-11ea-159c-331ff6660a75
# ╠═ace3dc76-f9ae-11ea-2bee-3d0bfa57cfbc
# ╟─0b67789c-f931-11ea-113c-35e5edafcbbf
# ╠═6896fef8-f9af-11ea-0065-816a70ba9670
# ╟─39152104-fc49-11ea-04dd-bb34e3600f2f
# ╟─e91c6fd8-f930-11ea-01ac-476bbde79079
# ╠═1b4c0c28-f9ab-11ea-03a6-69f69f7f90ed
# ╟─1f94e0a2-f9ab-11ea-1347-7dd906ebb09d
# ╠═41b2df7c-f931-11ea-112e-ede3b16f357a
# ╟─489fe282-f931-11ea-3dcb-35d4f2ac8b40
# ╟─1dd1e2f4-f930-11ea-312c-5ff9e109c7f6
# ╠═65c92cac-f930-11ea-20b1-6b8f45b3f262
# ╟─671525cc-f930-11ea-0e71-df9d4aae1c05
# ╟─4582ebf4-f930-11ea-03b2-bf4da1a8f8df
# ╠═6d8f47d2-2108-11eb-105e-63c5f7d74c6e
# ╠═b5fcd572-2108-11eb-3581-b18f259b42f6
# ╠═7898b76a-f930-11ea-2b7e-8126ec2b8ffd
# ╟─a5fbba46-f931-11ea-33e1-054be53d986c
# ╟─458cd100-f930-11ea-24b8-41a49f6596a0
# ╠═bc401bee-f931-11ea-09cc-c5efe2f11194
# ╟─ba695f6a-f931-11ea-0fbb-c3ef1374270e
# ╠═45c20988-f930-11ea-1d12-b782d2c01c11
# ╠═432276c6-2109-11eb-04e8-0bb0716eb660
# ╠═cc62929e-f9af-11ea-06b9-439ac08dcb52
# ╟─d3d7bd9c-f9af-11ea-1570-75856615eb5d
# ╟─2f8dedfc-fb98-11ea-23d7-2159bdb6a299
# ╟─b7446f34-f9b1-11ea-0f39-a3c17ba740e5
# ╟─4f97b572-f9b0-11ea-0a99-87af0797bf28
# ╟─46c905d8-f9b0-11ea-36ed-0515e8ed2621
# ╟─4e8d327e-f9b0-11ea-3f16-c178d96d07d9
# ╟─489b03d4-f9b0-11ea-1de0-11d4fe4e7c69
# ╟─d83f8bbc-f9af-11ea-2392-c90e28e96c65
# ╟─fd202410-f936-11ea-1ad6-b3629556b3e0
# ╟─0e465160-f937-11ea-0ebb-b7e02d71e8a8
# ╟─6718d26c-f9b0-11ea-1f5a-0f22f7ddffe9
# ╟─141af892-f933-11ea-1e5f-154167642809
# ╟─7eed9dde-f931-11ea-38b0-db6bfcc1b558
# ╟─7e3282e2-f931-11ea-272f-d90779264456
# ╟─7d1439e6-f931-11ea-2dab-41c66a779262
# ╠═7df55e6c-f931-11ea-33b8-fdc3be0b6cfa
# ╟─292e0384-fb57-11ea-0238-0fbe416fc976
# ╠═7dabee08-f931-11ea-0cb2-c7d5afd21551
# ╟─3736a094-fb57-11ea-1d39-e551aae62b1d
# ╠═13c89272-f934-11ea-07fe-91b5d56dedf8
# ╟─7d60f056-f931-11ea-39ae-5fa18a955a77
# ╟─b09f5512-fb58-11ea-2527-31bea4cee823
# ╟─8c7606f0-fb93-11ea-0c9c-45364892cbb8
# ╟─568f0d3a-fb54-11ea-0f77-171718ef12a5
# ╟─82e0df62-fb54-11ea-3fff-b16c87a7d45b
# ╠═b7601048-fb57-11ea-0754-97dc4e0623a1
# ╟─cc42de82-fb5a-11ea-3614-25ef961729ab
# ╠═d66fe2b2-fb5a-11ea-280f-cfb12b8296ac
# ╠═4ca8e04a-fb75-11ea-08cc-2fdef5b31944
# ╟─6f613cd2-fb5b-11ea-1669-cbd355677649
# ╠═91e87974-fb78-11ea-3ce4-5f64e506b9d2
# ╠═9f98e00e-fb78-11ea-0f6c-01206e7221d6
# ╟─d7d8cd0c-fb6a-11ea-12bf-2d1448b38162
# ╠═7be98e04-fb6b-11ea-111d-51c48f39a4e9
# ╠═c560ec86-211c-11eb-06e9-b3cc1c80e46d
# ╠═052f822c-fb7b-11ea-382f-af4d6c2b4fdb
# ╠═067f33fc-fb7b-11ea-352e-956c8727c79f
# ╟─954fc466-fb7b-11ea-2724-1f938c6b93c6
# ╟─e467c1c6-fbf2-11ea-0d20-f5798237c0a6
# ╟─7b10f074-fb7c-11ea-20f0-034ddff41bc3
# ╟─24ae5da0-fb7e-11ea-3480-8bb7b649abd5
# ╟─47836744-fb7e-11ea-2305-3fa5819dc154
# ╠═df4fc31c-fb81-11ea-37b3-db282b36f5ef
# ╠═c83b1770-fb82-11ea-20a6-3d3a09606c62
# ╟─52970ac4-fb82-11ea-3040-8bd0590348d2
# ╠═8ce3b312-fb82-11ea-200c-8d5b12f03eea
# ╠═a2214e50-fb83-11ea-3580-210f12d44182
# ╟─a9ffff9a-fb83-11ea-1efd-2fc15538e52f
# ╟─808abf6e-fb84-11ea-0785-2fc3f1c4a09f
# ╠═953363dc-fb84-11ea-1128-ebdfaf5160ee
# ╟─294b6f50-fb84-11ea-1382-03e9ab029a2d
# ╠═b726f824-fb5e-11ea-328e-03a30544037f
# ╠═18355314-fb86-11ea-0738-3544e2e3e816
# ╠═abe2b862-fb69-11ea-08d9-ebd4ba3437d5
# ╟─3d105742-fb8d-11ea-09b0-cd2e77efd15c
# ╟─a72fcf5a-fb62-11ea-1dcc-11451d23c085
# ╟─f83991c0-fb7c-11ea-0e6f-1f80709d00c1
# ╟─4b27a89a-fb8d-11ea-010b-671eba69364e
# ╟─d7b7a14a-fb90-11ea-3e2b-2fd8f379b4d8
# ╟─1939dbea-fb63-11ea-0bc2-2d06b2d4b26c
# ╟─70169682-fb8c-11ea-27c0-2dad2ff3080f
# ╟─b5dff8b8-fb6c-11ea-10fc-37d2a9adae8c
# ╟─402562b0-fb63-11ea-0769-375572cc47a8
# ╟─ee8c5808-fb5f-11ea-19a1-3d58217f34dc
# ╟─2521bac8-fb8f-11ea-04a4-0b077d77529e
# ╠═49b69dc2-fb8f-11ea-39af-030b5c5053c3
# ╟─7f341c4e-fb54-11ea-1919-d5421d7a2c75
# ╟─cc07f576-fbf3-11ea-2c6f-0be63b9356fc
# ╟─6b4d6584-f3be-11ea-131d-e5bdefcc791b
# ╟─54b1e236-fb53-11ea-3769-b382ef8b25d6
# ╟─b7803a28-fb96-11ea-3e30-d98eb322d19a
# ╟─ddef9c94-fb96-11ea-1f17-f173a4ff4d89
# ╟─ffc17f40-f380-11ea-30ee-0fe8563c0eb1
# ╟─ffc40ab2-f380-11ea-2136-63542ff0f386
# ╟─ffceaed6-f380-11ea-3c63-8132d270b83f
# ╟─ffde44ae-f380-11ea-29fb-2dfcc9cda8b4
# ╟─ffe326e0-f380-11ea-3619-61dd0592d409
# ╟─fff5aedc-f380-11ea-2a08-99c230f8fa32
# ╟─00026442-f381-11ea-2b41-bde1fff66011
# ╟─00115b6e-f381-11ea-0bc6-61ca119cb628

max_stars_repo_name: pankgeorg/learning

max_issues_repo_name: pankgeorg/learning

max_forks_repo_name: pankgeorg/learning


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
function g!(x, y)
    for i in 1:length(x)
        x[i] *= 3
        y[i] += x[i]
    end
    x, y
end

function h!()
    for i in 1:length(x_global)
        x_global[i] *= 3
        y_global[i] += x_global[i]
    end
end

function f0()
    x = randn(100, 100)
    y = randn(100, 100)
    x, y
end

function f1()
    x = randn(100, 100)
    y = randn(100, 100)
    g!(x, y)
end

function f2!(x = randn(100, 100), y = randn(100, 100))
    g!(x, y)
end

function f3!()
   g!(x_global, y_global)
end

function f4!()
   h!()
end

using BenchmarkTools

x_global = randn(100, 100)
y_global = randn(100, 100)

@show VERSION
@btime g!($x_global, $y_global)
@btime h!()
@btime f0()
@btime f1()
@btime f2!()
@btime f3!()
@btime f4!()

# %%
@code_warntype h!()

# %%
@code_warntype g!(x_global, y_global)

# %%
@code_typed f3!()

# %%
@code_typed g!(x_global, y_global)

# %%
@code_warntype f3!()

# %%

max_stars_repo_path: 0005/benchmark example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0005/benchmark example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0005/benchmark example.jl

max_forks_repo_name: genkuroki/public


text:
# 3D Objects
struct Cube <: Terminal
    size::AbstractArray{Real, 1}
    center::Bool
end

struct Sphere <: Terminal
    r::Real
end

struct Cylinder <: Terminal
    height::Real
    r1::Real
    r2::Real
    center::Bool
end

struct Polyhedron <: Terminal
    fields
end

toString(c::Cube) = "cube(size = [$(c.size[1]), $(c.size[2]), $(c.size[3])]);"
toString(s::Sphere) = "sphere($(s.r));"


text:
module Chrono

using TimeZones

export parsedate, parsedatetz
export Years, Months, Weeks, Days, Hours, Minutes, Seconds, Milliseconds
export Year, Month, Day, Weekday, NthWeekday, WeekdayLast, Monthday, MonthdayLast, MonthWeekday, MonthWeekdayLast
export YearMonth, YearMonthday, YearMonthdayLast, YearMonthWeekday, YearMonthWeekdayLast
export isok
export jan, feb, mar, apr, may, jun, jul, aug, sep, oct, nov, dec
export sun, mon, tue, wed, thu, fri, sat

const depfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depfile)
    include(depfile)
else
    error("Chrono not properly installed. Please run Pkg.build(\"Chrono\")")
end

include("parse.jl")
include("utils.jl")

end # module

max_stars_repo_name: tanmaykm/Chrono

max_issues_repo_name: tanmaykm/Chrono

max_forks_repo_name: tanmaykm/Chrono


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

using Markdown
using InteractiveUtils

# ╔═╡ 31b32634-6730-11eb-1eb6-d512e0afd93c
using Pkg, DrWatson

# ╔═╡ 417ce078-6730-11eb-0150-d515235783f0
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ 63ddbd40-6730-11eb-3a30-61121b8061b0
df = DataFrame(
	speed = [4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13,
		13, 14, 14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18,
		19, 19, 19, 20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25],
	dist = [2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34,
		34, 46, 26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84,
		36, 46, 68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85]
);

# ╔═╡ 85013b04-6731-11eb-3df7-4f9ba267f65b
PRECIS(df)

# ╔═╡ 97755284-6731-11eb-3a87-51384900e50b
cars_stan = "
data {
    int N;
    vector[N] speed;
    vector[N] dist;
}
parameters {
    real a;
    real b;
    real sigma;
}
transformed parameters{
    vector[N] mu;
    mu = a + b * speed;
}
model {
    a ~ normal(0, 100);
    b ~ normal(0, 10);
    sigma ~ exponential(1);
    dist ~ normal(mu, sigma)    ;
}
generated quantities {
    vector[N] log_lik;
    for (i in 1:N)
        log_lik[i] = normal_lpdf(dist[i] | mu[i], sigma);
}
";

# ╔═╡ a3a5ba84-6732-11eb-3051-afe0b875f51b
begin
cars_stan_model = SampleModel("cars.model", cars_stan)
	data = (N = size(df, 1), speed = df.speed, dist = df.dist)
	rc = stan_sample(cars_stan_model; data)

	if success(rc)
		stan_summary(cars_stan_model, true)
		nt_cars = read_samples(cars_stan_model);
	end
end;

# ╔═╡ 1d99951e-673c-11eb-3920-49db0451c1cd
begin
	post_df = read_samples(cars_stan_model, :dataframe)
	PRECIS(post_df[:, [:a, :b, :sigma]])
end

# ╔═╡ a3a5feae-6732-11eb-2893-7397c30b2daa
begin
	log_lik = nt_cars.log_lik'
	n_sam, n_obs = size(log_lik)
	lppds = reshape(lppd(log_lik), n_obs)
	sum(lppds)
end

# ╔═╡ 41909984-6734-11eb-1cd7-79ae89926d10
size(log_lik)

# ╔═╡ a3a6bcfe-6732-11eb-0b05-db55dac22d9c
begin
	pwaic = [var2(log_lik[:, i]) for i in 1:n_obs]
	-2(sum(lppds) - sum(pwaic))
end

# ╔═╡ 17b3bcb0-6741-11eb-1d7a-69218af68695
md"
!!! note
	Below WAIC value is identical to the value obtained by switching to rethinking's ulam() for the cars data."

# ╔═╡ a3b8541e-6732-11eb-14a2-515449707918
waic(log_lik)

# ╔═╡ a1ef0db2-6732-11eb-05e9-67cda4f8de09
begin
	loo, loos, pk = psisloo(log_lik)
	loo
end

# ╔═╡ 995fbd44-6733-11eb-216c-47405debe19d
sum(loos)

# ╔═╡ bb0ef2de-6733-11eb-0b0e-5339d80ffdfd
-2(loo - sum(pwaic))

# ╔═╡ 7b69ec38-6733-11eb-02a2-25ea26f08d22
pk_qualify(pk)

# ╔═╡ 83adace8-6733-11eb-33af-c3d49d9385b1
pk_plot(pk)

# ╔═╡ 0552e952-8385-11eb-1d2f-154960f187ba
md" ### End of clip-07-19-24s.jl"

# ╔═╡ Cell order:
# ╠═31b32634-6730-11eb-1eb6-d512e0afd93c
# ╠═417ce078-6730-11eb-0150-d515235783f0
# ╠═63ddbd40-6730-11eb-3a30-61121b8061b0
# ╠═85013b04-6731-11eb-3df7-4f9ba267f65b
# ╠═97755284-6731-11eb-3a87-51384900e50b
# ╠═a3a5ba84-6732-11eb-3051-afe0b875f51b
# ╠═1d99951e-673c-11eb-3920-49db0451c1cd
# ╠═a3a5feae-6732-11eb-2893-7397c30b2daa
# ╠═41909984-6734-11eb-1cd7-79ae89926d10
# ╠═a3a6bcfe-6732-11eb-0b05-db55dac22d9c
# ╟─17b3bcb0-6741-11eb-1d7a-69218af68695
# ╠═a3b8541e-6732-11eb-14a2-515449707918
# ╠═a1ef0db2-6732-11eb-05e9-67cda4f8de09
# ╠═995fbd44-6733-11eb-216c-47405debe19d
# ╠═bb0ef2de-6733-11eb-0b0e-5339d80ffdfd
# ╠═7b69ec38-6733-11eb-02a2-25ea26f08d22
# ╠═83adace8-6733-11eb-33af-c3d49d9385b1
# ╟─0552e952-8385-11eb-1d2f-154960f187ba


text:
include("../intcode_v1.jl")

orig_prog = parse_data("1,0,0,3,1,1,2,3,1,3,4,3,1,5,0,3,2,10,1,19,1,19,5,23,1,23,9,27,2,27,6,31,1,31,6,35,2,35,9,39,1,6,39,43,2,10,43,47,1,47,9,51,1,51,6,55,1,55,6,59,2,59,10,63,1,6,63,67,2,6,67,71,1,71,5,75,2,13,75,79,1,10,79,83,1,5,83,87,2,87,10,91,1,5,91,95,2,95,6,99,1,99,6,103,2,103,6,107,2,107,9,111,1,111,5,115,1,115,6,119,2,6,119,123,1,5,123,127,1,127,13,131,1,2,131,135,1,135,10,0,99,2,14,0,0")

prog = copy(orig_prog)
prog[2] = 12
prog[3] = 2
res=run_program!(prog)[1]
println(res)

for i = 1:100
   for j = 1:100
       prog = copy(orig_prog)
       prog[2]=i
       prog[3]=j
       res=run_program!(prog)[1]
       if res == 19690720
           println(i, j)
	   end
   end
end

max_stars_repo_name: st--/advent-of-code

max_issues_repo_name: st--/advent-of-code

max_forks_repo_name: st--/advent-of-code


text:
@testset "FM90" begin

    x_inv_microns = [3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 9.0, 10.0]
    wave = 1e4 ./ x_inv_microns

    model = FM90()
    model1 = FM90([0.10 0.70 3.23 0.41])

    # Test out of bounds
    bad_waves = [100, 4e4]
    @test model.(bad_waves) == zeros(length(bad_waves))

    ref_values = [2.9829317, 4.1215415, 6.4135842, 5.6574243, 4.7573250, 5.4905843, 9.2853567, 12.462238]
    @test model.(wave) ≈ ref_values rtol = 1e-4
    @test model.(wave) == model1.(wave)

    # uncertainties
    noise = randn(length(wave)) .* 0.01
    wave_unc = wave .± noise
    reddening = @inferred broadcast(model, wave_unc)
    reddening1 = @inferred broadcast(model1, wave_unc)
    @test Measurements.value.(reddening) ≈ ref_values rtol = 1e-4
    @test Measurements.value.(reddening1) ≈ ref_values rtol = 1e-4

    # Unitful
    wave_u = wave * u"angstrom"
    reddening = @inferred broadcast(model, wave_u)
    reddening1 = @inferred broadcast(model1, wave_u)
    @test eltype(reddening) <: Gain
    @test eltype(reddening1) <: Gain
    @test ustrip.(reddening) ≈ ref_values rtol = 1e-4
    @test ustrip.(reddening1) ≈ ref_values rtol = 1e-4
end


text:
struct ExtensionalArcFunction{T<:Real}
    components      ::Dict{Tuple{Int, Int}, T}
end

struct ExtensionalArcObjective<:ObjectiveFunction
    f           ::ExtensionalArcFunction
    type        ::ObjectiveType
end

"""
    evaluateDecision(obj::ExtensionalArcObjective, i1::Int, i2::Int)

Takes in a decision and returns the impact on the function

# Arguments
- `obj::ExtensionalArcObjective`: the function to be evaluated
- `i1::Int`: the index of the first point
- `i2::Int`: the index of the second point
"""
function evaluateDecision(obj::ExtensionalArcObjective, i1::Int, i2::Int)
    return obj.f.components[(i1, i2)]
end

"""
    evaluateDecisions(obj::ExtensionalArcObjective,decisions::Vector{Int})

Takes in a list of decisions and returns the impact on the function

# Arguments
- `obj::ExtensionalArcObjective`: the function to be evaluated
- `decisions::Vector{Int}`: the list
"""
function evaluateDecisions(obj::ExtensionalArcObjective,decisions::Vector{Int})
    sum = 0
    for i in 1:length(decisions)-1
        sum += evaluateDecision(obj, decisions[i], decisions[i+1])
    end
    return sum
end
"""
    hasKey(obj::ExtensionalArcObjective,i1::Int, i2::Int)

Takes in an arc and return its validity

# Arguments
- `obj::ExtensionalArcObjective`: the function to be evaluated
- `i1::Int`: the index of the first point
- `i2::Int`: the index of the second point
"""
function hasKey(obj::ExtensionalArcObjective,i1::Int, i2::Int)
    return haskey(obj.f.components,(i1, i2))
end


text:
function predict(x_curr, Ah, Qh, PI=I)
    mean = Ah * x_curr.μ
    cov = Ah * x_curr.Σ * Ah' .+ Qh

    return Gaussian(mean, cov)
end

function update(x_pred::Gaussian, h::AbstractVector, H::AbstractMatrix, R::AbstractMatrix, PI=I)
    m_p, P_p = x_pred.μ, x_pred.Σ

    # If the predicted covariance is zero, the prediction will not be adjusted!
    v = 0 .- h
    S = Symmetric(H * P_p * H' .+ R)
    S_inv = inv(S)
    K = P_p * H' * S_inv

    filt_mean = m_p .+ P_p * H' * (S\v)
    # KSK = P_p * H' * (S \ (H * P_p'))
    # filt_cov = P_p .- KSK
    filt_cov = X_A_Xt(PDMat(Symmetric(P_p)), (I-K*H))
    if !iszero(R)
        filt_cov .+= Symmetric(X_A_Xt(PDMat(R), K))
    end

    assert_nonnegative_diagonal(filt_cov)

    return Gaussian(filt_mean, filt_cov)
end

function smooth(x_curr::Gaussian, x_next_smoothed::Gaussian, Ah::AbstractMatrix, Qh::AbstractMatrix,
                PI=I)
    # @info "smooth" x_curr x_next_pred x_next_smoothed Ah
    P_p = Symmetric(Ah * x_curr.Σ * Ah' .+ Qh)
    P_p_inv = inv(P_p)

    G = x_curr.Σ * Ah' * P_p_inv

    smoothed_mean = x_curr.μ + G * (x_next_smoothed.μ - Ah*x_curr.μ)

    # x_next_diff_cov = x_next_smoothed.Σ - x_next_pred.Σ
    # GDG = Gain * (x_next_diff_cov) * Gain'
    # smoothed_cov = x_curr.Σ + GDG

    P = copy(x_curr.Σ)
    C_tilde = Ah
    K_tilde = P * Ah' * P_p_inv
    smoothed_cov = ((I - K_tilde*C_tilde) * P * (I - K_tilde*C_tilde)'
                    + K_tilde * Qh * K_tilde' + G * x_next_smoothed.Σ * G')

    assert_nonnegative_diagonal(smoothed_cov)

    x_curr_smoothed = Gaussian(smoothed_mean, smoothed_cov)
    return x_curr_smoothed, G
end

function assert_nonnegative_diagonal(cov)
    if !all(diag(cov) .>= 0)
        @error "Non-positive variances" cov diag(cov)
        error("The provided covariance has non-positive entries on the diagonal!")
    end
end

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
include("../src/FDTD2D.jl")
import .FDTD2D
using Plots
using Test

# Test FDTD2D.calculate_params
function laser_pulse_gauss(x_min, width, duration)
  function f(t, x, y)
    return exp(-(y/width)^2)*exp(-((t - (x - x_min) - duration)/duration)^2)*sin(2.0*pi*(t - x))
  end
  return f
end

params_correct = Dict(
  "x_bounds" => (0.0, 12.0),
  "y_bounds" => (-4.0, 4.0),
  "matrix_size" => Dict("x" => 96, "y" => 64),
  "time_steps" => 81,

  "output" => Dict(
    "iteration_pass" => 10,
    "directory_name" => "test/output"
  )
)
params_correct["laser_pulse_y_shape"] = laser_pulse_gauss(params_correct["x_bounds"][1], 1., 2.)
params_correct["laser_pulse_z_shape"] = laser_pulse_gauss(params_correct["x_bounds"][1], 1., 2.)

params_calculated = FDTD2D.calculate_params(params_correct)

@testset "Calculation of parameters" begin
    @test params_calculated["box_size"]["x"] == 12.0
    @test params_calculated["box_size"]["y"] == 8.0
    @test params_calculated["space_step"]["x"] == 0.125
    @test params_calculated["space_step"]["y"] == 0.125
    @test params_calculated["time_step"] == sqrt(2.0)/16.0
    @test params_calculated["half_cfl"]["x"] == 0.5*sqrt(0.5)
    @test params_calculated["half_cfl"]["y"] == 0.5*sqrt(0.5)
end

# Test FDTD2D.init_data
data_empty = Dict()
for name in ("ex","ey","ezx","ezy","ez")
  data_empty[name] = zeros(Float64, params_calculated["matrix_size"]["x"], params_calculated["matrix_size"]["y"])
end
for name in ("hx","hy","hzx","hzy","hz")
  data_empty[name] = zeros(Float64, params_calculated["matrix_size"]["x"]-1, params_calculated["matrix_size"]["y"]-1)
end

data_test = FDTD2D.init_data(params_calculated)
@test isequal(data_test,data_empty)

@time d = FDTD2D.run(params_correct)

heatmap(d["ez"])


text:
# Inverse and soln of matrix
# Gaussian
# generally JL handles eqns and invs

# we use Finverse to denote function to find inverse
# 🙈

function Finverse_Gaussian(x)
    n = Int(sqrt(length(x)))
    x = reshape(x,n,n)
    I = zeros(n,n)
    for i in 1:n
        I[i,i] =1
    end
    for i in 1:(n-1)
        for j in (i+1):n
            m = x[j,i]/x[i,i]
            x[j,i:end] = x[j,i:end].-x[i,i:end].*m
            I[j,:] = I[j,:].-I[i,:].*m
        end
    end
    for i in n:-1:2
        for j in (i-1):-1:1
            m = x[j,i]/x[i,i]
            x[j,i]=0
            I[j,:] =I[j,:] .- I[i,:].*m
        end
    end
    for i in 1:n
        I[i,:] = I[i,:]./x[i,i]
    end
    return I
end

function Fsolve_Gaussian(x,b)
    n=length(b)
    x = reshape(x,n,n)
    for i in 1:(n-1)
        for j in (i+1):n
            m = x[j,i]/x[i,i]
            x[j,i:end] = x[j,i:end].-x[i,i:end].*m
            b[j] = b[j]-b[i]*m
        end
    end
    return b
end

Fsolve_Gaussian([1 2 3; 4 5 6;7 8 9],[1;2;3])
[1 2 3; 4 5 6;7 8 9]*[1 ,-2, 0]

#Finverse_Gaussian([1 2;6 7])

#[1 2;6 7]*Finverse_Gaussian([1 2;6 7])

# safer version of Gaussian solve -judge and shift terms with 0

#function Fsolve_Gaussian_safely(x,b)
#    n=length(b)
#    for j in 1:(n-1)
#    end
#end

x = [1 2 3;4 5 6;7 8 9]

function LU(x)
    n=Int(sqrt(length(x)))
    for j in 1:(n-1)
        for i in (j+1):n
            m = x[i,j] / x[j,j]
            x[i,(j+1):n] = x[i,(j+1):n].-m.*x[j,(j+1):n]
            x[i,j] = m
        end
    end
    return x
end
LU(x)


text:
# Define a generic function
function nanoread end

# SequencingSummary Reader
function nanoread(input::DataFrame)
	return DataFrame(quality=round.(input.mean_qscore_template, digits=1), length=input.sequence_length_template)
end

# Fastq File Reader
function nanoread(input::FASTX.FASTQ.Reader)
	records = collect(input);
	remove_bad_read!(records);
	recordnum = length(records);
	readsinfo = Array{Union{FastqInfo, Missing},1}(missing, recordnum);
	Threads.@threads for i in 1:recordnum
		@inbounds readsinfo[i] = get_info(records[i])
	end
	return extract_info(readsinfo)
end

# BAM File Reader
function nanoread(input::XAM.BAM.Reader)
    readsinfo = BAMInfo[]
    record = BAM.Record()
    while !eof(input)
        empty!(record)
        read!(input, record)
        push!(readsinfo, get_info(record))
    end
	return extract_info(readsinfo)
end

# Miltiple-reads FAST5 File Reader
function nanoread(input::HDF5File, h5ver::String, basecallGroup::String)
	dirpath = "Analyses/$basecallGroup/Summary/basecall_1d_template"
    readIDs = names(input)
	readsnum = length(readIDs)
	readsinfo = Array{Union{FastqInfo, Missing},1}(missing, readsnum)
	@inbounds for i in 1:readsnum
		thisread = readIDs[i]
		readrecord = input[thisread]
		summaryInfo = attrs(readrecord[dirpath])
		readsinfo[i] = get_info(summaryInfo, readrecord, basecallGroup)
	end
	filter!(e -> e !== missing, readsinfo)
	return extract_info(readsinfo)
end

# Single-read FAST5 File Reader
function nanoread(input::HDF5File, h5ver::Float64, basecallGroup::String)
	dirpath = "Analyses/$basecallGroup/Summary/basecall_1d_template"
	summaryInfo = attrs(input[dirpath])
	readinfo = get_info(summaryInfo, input, h5ver, basecallGroup)
	return DataFrame(quality = readinfo.quality, length = readinfo.length, gc_content = readinfo.gc)
end


text:
function correct_guess1( nlp, pars, timer, x, a, J, g, s, y )
      start_advanced_timer(timer, "INIT/centre")

      ais = cons_indicies(nlp)

      y = y / pars.mu_scale
      mu = dot(s,y)[1] / length(y)

      for i = 1:10
        @show size(s), size(y)
        @show mu, norm(s,Inf)

        y_c = mu ./ s
        buffer = 2.0
        y = min.( y_c / (pars.ls.comp_feas * buffer), max.(y, pars.ls.comp_feas * y_c * buffer)) # project onto complementarity constraints

        #Delta_s = norm(g - J' * y,1) / (length(y) + norm(y,1))
        Delta_s = norm(g - J' * y,Inf) / (1.0 + norm(y,Inf))

        primal_feas_avg = norm(a - s,1) / length(s)
        @show primal_feas_avg, Delta_s, norm(y,Inf)
        if Delta_s > primal_feas_avg
          if pars.start_satisfying_bounds
            s[ais] = s[ais] + Delta_s
          else
            s = s + Delta_s
          end

          mu *= 10.0
        else
          break
        end
      end

      pause_advanced_timer(timer, "INIT/centre")

      #Delta_s = norm(g - J' * y,Inf) / (1.0 + 0.01 * norm(y,Inf))
      #@show Delta_s

      start_advanced_timer(timer, "INIT/construct_class")
      init_point = Class_point(x, y, s, mu)
      check_for_nan(init_point)
      init_it = Class_iterate(init_point, nlp, Class_local_info(), timer, pars);
      pause_advanced_timer(timer, "INIT/construct_class")

      return init_it
end

function correct_guess2( nlp, pars, timer, x, a, J, g, y, mu, conWeight )
      start_advanced_timer(timer, "INIT/centre")

      ais = cons_indicies(nlp)
      s = zeros(length(a))

      i = 1
      min_s = deepcopy(s)
      for i = 1:100
        @show mu
        s = a + mu * conWeight
        #@assert()
        y_c = mu ./ s
        buffer = 2.0
        y = min.( y_c / (pars.ls.comp_feas * buffer), max.(y, pars.ls.comp_feas * y_c * buffer)) # project onto complementarity constraints

        min_s = (mu / 20.0) * conWeight #+ conWeight * max(0.0,-minimum(a))

        if all(s .>= min_s) && norm(g - J' * y,1) / (length(s) + norm(y,1)) < mu * pars.init.dual_threshold
          break
        else
          mu *= 2.0
        end
      end

      if i == 100
        @show mu, norm(conWeight), all(s .>= min_s), minimum(a)
        @show norm(g - J' * y,1) / (length(s) + norm(y,1))
        error("init error max it")
      end

      pause_advanced_timer(timer, "INIT/centre")

      start_advanced_timer(timer, "INIT/construct_class")
      init_point = Class_point(x, y, s, mu)
      check_for_nan(init_point)
      init_it = Class_iterate(init_point, nlp, Class_local_info(), timer, pars);
      pause_advanced_timer(timer, "INIT/construct_class")

      return init_it
end

function correct_guess3( nlp, pars, timer, x::Vector, a::Vector, J, g::Vector, y::Vector, mu::Float64, conWeight )
      ais = cons_indicies(nlp)
      bis = bound_indicies(nlp)

      #if norm(g - J' * y,1) >= norm(g,1)
      #  y = 1e-10 * ones(length(y))
      #end

      #n = length(x)
      #tau = norm(g - J' * y,1) / ( n * mu)
      #tau = norm(g - J' * y,1) / ( ncon(nlp) * mu)
      #if tau > 1.0
      #  mu = mu * tau
      #end

      s = a + mu * conWeight
      #s[bis] = a[bis]
      #@show bis, ais
      @assert(all(s[bis] .>= 0.0))
      @assert(all(s[ais] .>= 0.0))
      mu *= pars.init.mu_scale

      #@assert(mu >= 0.0)

      start_advanced_timer(timer, "INIT/construct_class")
      init_point = Class_point(x, y, s, mu)
      check_for_nan(init_point)
      init_it = Class_iterate(init_point, nlp, Class_local_info(), timer, pars);
      pause_advanced_timer(timer, "INIT/construct_class")

      center_dual!(init_it, pars)

      return init_it
end