text:
# Cartesian projection (concrete examples include Lambert or
# EquiRect), defined by the fact that the pixels are stored in a
# matrix. The projection type must have a constructor accepting at
# least (Ny,Nx,T,storage) keyword arguments.
abstract type CartesianProj <: Proj end
make_field_aliases("Cartesian", CartesianProj)

### constructors

_reshape_batch(arr::AbstractArray{T,3}) where {T} = reshape(arr, size(arr,1), size(arr,2), 1, size(arr,3))
_reshape_batch(arr) = arr

## constructing from arrays
# spin-0
function (::Type{F})(Ix::A; Proj=default_proj(F), kwargs...) where {F<:BaseField{Map,<:CartesianProj}, T, A<:AbstractArray{T}}
    BaseField{Map}(_reshape_batch(Ix), Proj(;Ny=size(Ix,1), Nx=size(Ix,2), T, storage=basetype(A), kwargs...))
end
function (::Type{F})(Il::A; Ny, Proj=default_proj(F), kwargs...) where {F<:BaseField{Fourier,<:CartesianProj}, T, A<:AbstractArray{T}}
    BaseField{Fourier}(_reshape_batch(Il), Proj(;Ny, Nx=size(Il,2), T, storage=basetype(A), kwargs...))
end
# spin-2
function (::Type{F})(X::A, Y::A, metadata::CartesianProj) where {B<:Basis2Prod{<:Union{𝐐𝐔,𝐄𝐁}}, F<:BaseField{B,<:CartesianProj}, T, A<:AbstractArray{T}}
    BaseField{B}(cat(_reshape_batch(X), _reshape_batch(Y), dims=Val(3)), metadata)
end
function (::Type{F})(X::A, Y::A; Proj=default_proj(F), kwargs...) where {B<:Basis2Prod{<:Union{𝐐𝐔,𝐄𝐁}}, F<:BaseField{B,<:CartesianProj}, T, A<:AbstractArray{T}}
    BaseField{B,<:CartesianProj}(X, Y, Proj(;Ny=size(X,1), Nx=size(X,2), T, storage=basetype(A), kwargs...))
end
# spin-(0,2)
function (::Type{F})(X::A, Y::A, Z::A, metadata::CartesianProj) where {B<:Basis3Prod{𝐈,<:Union{𝐐𝐔,𝐄𝐁}}, F<:BaseField{B,<:CartesianProj}, T, A<:AbstractArray{T}}
    BaseField{B}(cat(_reshape_batch(X), _reshape_batch(Y), _reshape_batch(Z), dims=Val(3)), metadata)
end
function (::Type{F})(X::A, Y::A, Z::A; Proj=default_proj(F), kwargs...) where {B<:Basis3Prod{𝐈,<:Union{𝐐𝐔,𝐄𝐁}}, F<:BaseField{B,<:CartesianProj}, T, A<:AbstractArray{T}}
    BaseField{B,<:CartesianProj}(X, Y, Z, Proj(;Ny=size(X,1), Nx=size(X,2), T, storage=basetype(A), kwargs...))
end

## constructing from other fields
function (::Type{F})(X::BaseField{B₀,P}, Y::BaseField{B₀,P}) where {B₀<:Union{Map,Fourier}, B<:Basis2Prod{<:Union{𝐐𝐔,𝐄𝐁},B₀}, P<:CartesianProj, F<:BaseField{B,<:CartesianProj}}
    BaseField{B,P}(cat(X.arr, Y.arr, dims=Val(3)), get_metadata_strict(X, Y)) :: F
end
function (::Type{F})(X::BaseField{B₀,P}, Y::BaseField{Basis2Prod{Pol,B₀},P}) where {B₀<:Union{Map,Fourier}, Pol<:Union{𝐐𝐔,𝐄𝐁}, B<:Basis3Prod{𝐈,Pol,B₀}, P<:CartesianProj, F<:BaseField{B,<:CartesianProj}}
    BaseField{B,P}(cat(X.arr, Y.arr, dims=Val(3)), get_metadata_strict(X, Y)) :: F
end
function (::Type{F})(X::BaseField{B₀,P}, Y::BaseField{B₀,P}, Z::BaseField{B₀,P}) where {B₀<:Union{Map,Fourier}, B<:Basis3Prod{𝐈,<:Union{𝐐𝐔,𝐄𝐁},B₀}, P<:CartesianProj, F<:BaseField{B,<:CartesianProj}}
    BaseField{B,P}(cat(X.arr, Y.arr, Z.arr, dims=Val(3)), get_metadata_strict(X, Y, Z)) :: F
end

# consistency checks called from BaseField constructor
function check_field_consistency(::SpatialBasis{Map}, arr, proj::P) where {P<:CartesianProj}
    size(arr)[1:2] == (proj.Ny, proj.Nx) || error("$(basetype(P)) metadata (Ny=$(proj.Ny), Nx=$(proj.Nx)) is inconsistent with pixel array size $(size(arr)).")
end
function check_field_consistency(::SpatialBasis{Fourier}, arr, proj::P) where {P<:CartesianProj}
    size(arr)[1:2] == (proj.Ny÷2+1, proj.Nx) || error("$(basetype(P)) metadata (Ny=$(proj.Ny), Nx=$(proj.Nx)) is inconsistent with half-plane Fourier array size $(size(arr)).")
end

### array interface
# most is inherited from BaseField. the main thing we have specify
# here has to do with which dimension is the "batch" dimension
# (dimension 4), since that is not assumed in BaseField
similar(f::CartesianField{B}, Nbatch::Int) where {B} = CartesianField{B}(similar(f.arr, size(f.arr,1), size(f.arr,2), size(f.arr,3), Nbatch), f.metadata)
nonbatch_dims(f::CartesianField) = ntuple(identity,min(3,ndims(f.arr)))
require_unbatched(f::CartesianField) = (f.Nbatch==1) || error("This function not implemented for batched fields.")
pol_slice(f::CartesianField, i) = (:, :, i, ..)

### properties
# generic
getproperty(f::CartesianField, ::Val{:Nbatch}) = size(getfield(f,:arr), 4)
getproperty(f::CartesianField, ::Val{:Npol})   = size(getfield(f,:arr), 3)
getproperty(f::CartesianField, ::Val{:T})      = eltype(f)
getproperty(f::CartesianField, ::Val{:proj})   = getfield(f, :metadata)

### indices
function getindex(f::CartesianS0, k::Symbol; full_plane=false)
    maybe_unfold = full_plane ? x->unfold(x,fieldinfo(f).Ny) : identity
    @match k begin
        :I  => f
        :Ix => Map(f).Ix
        :Il => maybe_unfold(Fourier(f).Il)
        _   => throw(ArgumentError("Invalid CartesianS0 index: $k"))
    end
end
function getindex(f::CartesianS2{Basis2Prod{B₁,B₂}}, k::Symbol; full_plane=false) where {B₁,B₂}
    maybe_unfold = (full_plane && k in [:El,:Bl,:Ql,:Ul]) ? x->unfold(x,fieldinfo(f).Ny) : identity
    B = @match k begin
        (:P)         => identity
        (:E  || :B)  => Basis2Prod{𝐄𝐁,B₂}
        (:Q  || :U)  => Basis2Prod{𝐐𝐔,B₂}
        (:Ex || :Bx) => Basis2Prod{𝐄𝐁,Map}
        (:El || :Bl) => Basis2Prod{𝐄𝐁,Fourier}
        (:Qx || :Ux) => Basis2Prod{𝐐𝐔,Map}
        (:Ql || :Ul) => Basis2Prod{𝐐𝐔,Fourier}
        _ => throw(ArgumentError("Invalid CartesianS2 index: $k"))
    end
    maybe_unfold(getproperty(B(f),k))
end
function getindex(f::CartesianS02{Basis3Prod{B₁,B₂,B₃}}, k::Symbol; full_plane=false) where {B₁,B₂,B₃}
    maybe_unfold = (full_plane && k in [:Il,:El,:Bl,:Ql,:Ul]) ? x->unfold(x,fieldinfo(f).Ny) : identity
    B = @match k begin
        (:I  || :P)  => identity
        (:E  || :B)  => Basis3Prod{𝐈,𝐄𝐁,B₃}
        (:Q  || :U)  => Basis3Prod{𝐈,𝐐𝐔,B₃}
        (:Ix)        => Basis3Prod{𝐈,B₂,Map}
        (:Il)        => Basis3Prod{𝐈,B₂,Fourier}
        (:Ex || :Bx) => Basis3Prod{𝐈,𝐄𝐁,Map}
        (:El || :Bl) => Basis3Prod{𝐈,𝐄𝐁,Fourier}
        (:Qx || :Ux) => Basis3Prod{𝐈,𝐐𝐔,Map}
        (:Ql || :Ul) => Basis3Prod{𝐈,𝐐𝐔,Fourier}
        _ => throw(ArgumentError("Invalid CartesianS02 index: $k"))
    end
    maybe_unfold(getproperty(B(f),k))
end


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
#
# Contributed by Kristoffer Carlsson
# Updated by Adam Beckmeyer

import .Threads: @spawn

const LINESIZE = 60
const CHUNKSIZE = 512LINESIZE

const IM = 139968 % Int32
const IA = 3877 % Int32
const IC = 29573 % Int32
const SEED = Ref(42 % Int32)

gen_random() = (SEED[] = ((SEED[] * IA) + IC) % IM)

# Vector{UInt8} faster indexing than Base.CodeUnits{UInt8,String} (julia 1.2)
const alu = Vector{UInt8}(string(
   "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG",
   "GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA",
   "CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT",
   "ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA",
   "GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG",
   "AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC",
   "AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA"))

const iub1 = Vector{UInt8}("acgtBDHKMNRSVWY")
# Instead of the probability, iub2 is the integer upper-bounding the random
# integer to the equivalent index in iub1. cumsum is not used here as
# it adds almost 10% of total time to inference time.
const iub2 = [floor(Int32, IM * x) for x in
              (0.27, 0.39, 0.51, 0.78, 0.8, 0.82, 0.84, 0.86, 0.88, 0.9, 0.92,
               0.94, 0.96, 0.98, 1.0)]
const homosapiens1 = Vector{UInt8}("acgt")
const homosapiens2 = [floor(Int32, IM * x) for x in
                      (0.3029549426680, 0.5009432431601, 0.6984905497992, 1.0)]

# Defining a function with successive comparisons is faster than getindex lookup
for (op, symb, pr) in ((:iub, iub1, iub2),
                       (:homosapiens, homosapiens1, homosapiens2))
    # successively compares x with all values in pr and returns value from symb
    @eval @inline function $op(x)
        $([:(x <= $(pr[i]) && return $(symb[i])) for i=1:length(symb)]...)
        'a' % UInt8
    end
end

function repeat_fasta(io, src, n)
    len = length(src)
    pos = 1
    GC.@preserve src for i=1:LINESIZE:n
        # l is the number of characters we need to write this iteration.
        # m is less than l if not enough characters are left to write from src.
        l = min(LINESIZE, n - i + 1)
        m = min(l, len - pos + 1)

        unsafe_write(io, pointer(src, pos), m)
        if l != m # reached end of src
            # Loop back around to the start of src and write remaining chars
            pos = 1
            m = l - m
            unsafe_write(io, pointer(src, pos), m)
        end
        pos += m

        write(io, '\n')
    end
end

# Writes chunks of random characters from c to io with newlines added
# every LINESIZE. c is filled with Task which must be ordered. The
# chunks yielded from each Task must be a multiple in length of
# LINESIZE unless they are the last chunk in the channel.
function print_chunks(io, c)
    for vt in c
        v = fetch(vt)::Vector{UInt8}
        n = length(v)
        for i=1:LINESIZE:n
            GC.@preserve v unsafe_write(io, pointer(v, i),
                                        min(LINESIZE, n - i + 1))
            write(io, '\n')
        end
    end
end

# Put tasks that yield chunks of random characters into c. The total
# number of characters before c is closed will be n.
function generate_chunks(flookup, c, n)
    for p=1:CHUNKSIZE:n
        # Create a chunk of random Int32 of length l
        l = min(CHUNKSIZE, n - p + 1)
        v = Vector{Int32}(undef, l)
        # Fill chunk with random numbers
        for i=1:l
            @inbounds v[i] = gen_random()
        end

        # Spawn task to convert random Int32 to characters. We must
        # put this task in the channel synchronously to make sure we
        # preserve the order of writing.
        t = @spawn @inbounds map!(flookup, Vector{UInt8}(undef, l), v)
        put!(c, t)
    end
    close(c)
end

function random_fasta(io, flookup, n, firstline, iotask, rtask=nothing)
    c = Channel{Task}(24)

    # We cannot proceed with generating random numbers until rtask is
    # complete otherwise the random numbers will be interleaved.
    rtask === nothing || wait(rtask)
    rtask = @spawn generate_chunks(flookup, c, n)
    # Similarily, we can't start writing to io until the previous
    # iotask has completed.
    iotask2 = @spawn begin
        wait(iotask)
        write(io, firstline)
        print_chunks(io, c)
    end

    iotask2, rtask
end

function fasta(n=25000000, io=stdout)
    write(io, ">ONE Homo sapiens alu\n")
    io1 = @spawn repeat_fasta(io, alu, 2n)

    io2, r1 = random_fasta(io, iub, 3n, ">TWO IUB ambiguity codes\n", io1)
    io3, _ = random_fasta(io, homosapiens, 5n, ">THREE Homo sapiens frequency\n",
                          io2, r1)
    wait(io3)
end

function real_main()
    isinteractive() || fasta(parse(Int,ARGS[1]))
end

if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end

ext: jl


text:
@testset "irtracer: calls" begin
    val, tape = irtrace(inc_mul, 2.0, 3.0)
    @test val == inc_mul(2.0, 3.0)
    @test length(tape) == 5
    @test tape[3] isa Constant
end

@testset "irtracer: bcast" begin
    A = rand(3)
    B = rand(3)
    val, tape = irtrace(inc_mul, A, B)
    @test val == inc_mul(A, B)
    # broadcasting may be lowered to different forms,
    # so making no assumptions regarding the tape

    val, tape = irtrace(inc_mul2, A, B)
    @test val == inc_mul2(A, B)
end

@testset "irtracer: primitives" begin
    x = 3.0
    val1, tape1 = irtrace(non_primitive_caller, x)
    val2, tape2 = irtrace(non_primitive_caller, x; primitives=Set([non_primitive, sin]))

    @test val1 == val2
    @test any(op isa Call && op.fn == (*) for op in tape1)
    @test tape2[2].fn == non_primitive
    @test tape2[3].fn == sin

end


text:
# This file was generated, do not modify it. # hide
r1 = range(xgb, :max_depth, lower=3, upper=10)
r2 = range(xgb, :min_child_weight, lower=0, upper=5)

tm = TunedModel(model=xgb, tuning=Grid(resolution=8),
                resampling=CV(rng=11), ranges=[r1,r2],
                measure=cross_entropy)
mtm = machine(tm, X, y)
fit!(mtm, rows=train)


text:
root = joinpath(@__DIR__, "..")
using Pkg; Pkg.activate(root)

src = joinpath(root, "src")
out = joinpath(root, "notebooks")

using Literate

mkpath(out)

for f in ["Project.toml", "Manifest.toml"]
  cp(joinpath(root, f), joinpath(out, f), force = true)
end

function preprocess(s)
  s = "using Pkg; Pkg.activate(\".\"); Pkg.instantiate()\n#-\n" * s
end

for f in ["utils.jl"]
  cp(joinpath(src, f), joinpath(out, f), force = true)
end

for f in ["intro.jl", "backandforth.jl", "forward.jl", "tracing.jl", "reverse.jl"]
  Literate.notebook(joinpath(src, f), out,
                    preprocess = preprocess,
                    credit = false)
end


text:
stress_scale = 1.0e+5
strain_scale = 1

force_scale = 5.0
tid = 300

testtype = "NeuralNetwork2D"

nntype = "piecewise"

include("nnutil.jl")
printstyled("force_scale=$force_scale, tid=$tid\n", color=:green)

H0 = [1335174.0968380707 326448.3267263398   0.0
      326448.3267263398  1326879.2022994285  0.0
      0.0                0.0                 526955.763626241]/stress_scale

H0inv = inv(H0)

nn_file = "Data/$(nntype)/nn_train_false_$(idx)_$(H_function)_from3_ite50.mat"
#nn_file = "Data/$(nntype)/NNPreLSfit_$(idx)_$(H_function)_2.mat"

@show nn_file
s = ae_to_code(nn_file, nntype)

eval(Meta.parse(s))

# density 4.5*(1 - 0.25) + 3.2*0.25
fiber_fraction = 0.25
#todo
#fiber_fraction = 1.0
prop = Dict("name"=> testtype, "rho"=> 4.5*(1 - fiber_fraction) + 3.2*fiber_fraction, "nn"=>post_nn)

T = 200.0
NT = 200

# nx_f, ny_f = 12, 4
# homogenized computaional domain
# number of elements in each directions
nx, ny = 10, 5

porder = 2

nodes, EBC, g, gt, FBC, fext, ft = BoundaryCondition(tid, nx, ny,porder; force_scale=force_scale)

ndofs=2
elements = []
for j = 1:ny
    for i = 1:nx
        n = (nx*porder+1)*(j-1)*porder + (i-1)porder+1
        #element (i,j)
        if porder == 1
            #   4 ---- 3
            #
            #   1 ---- 2

            elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        elseif porder == 2
            #   4 --7-- 3
            #   8   9   6
            #   1 --5-- 2
            elnodes = [n, n + 2, n + 2 + 2*(2*nx+1),  n + 2*(2*nx+1), n+1, n + 2 + (2*nx+1), n + 1 + 2*(2*nx+1), n + (2*nx+1), n+1+(2*nx+1)]
        else
            error("polynomial order error, porder= ", porder)
        end
        coords = nodes[elnodes,:]
        push!(elements,SmallStrainContinuum(coords,elnodes, prop, 3))
    end
end

domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs)
, zeros(domain.neqs),∂u, domain.neqs, gt, ft)

assembleMassMatrix!(globdat, domain)
updateStates!(domain, globdat)

ts = LinRange(0, T, NT+1)
adaptive_solver_args = Dict("Newmark_rho"=> 0.0,
                          "Newton_maxiter"=>10,
                          "Newton_Abs_Err"=>1e-4,
                          "Newton_Rel_Err"=>1e-6,
                          "damped_Newton_eta" => 1.0)

globdat, domain, ts = AdaptiveSolver("NewmarkSolver", globdat, domain, T, NT, adaptive_solver_args)

# plot
#=
close("all")
visσ(domain)
savefig("Debug/test$tid.png")

close("all")
visσ(domain)
axis("equal")
savefig("Debug/order$porder/test_stress$(tid)_$force_scale.png")

close("all")
ux = [reshape(domain.history["state"][i][1:(nx*porder+1)*(ny*porder+1)], ny*porder+1, nx*porder+1)[1,end] for i = 1:length(domain.history["state"])]
plot(ts, ux)
savefig("Debug/order$porder/test_ux_$(use_reg)_$(idx)_$H_function.png")

close("all")
uy = [reshape(domain.history["state"][i][(nx*porder+1)*(ny*porder+1)+1:end], ny*porder+1, nx*porder+1)[1,end] for i = 1:length(domain.history["state"])]
plot(ts, uy)
savefig("Debug/order$porder/test_uy_$(use_reg)_$(idx)_$H_function.png")
=#


text:
using MeshArrays, MITgcmTools
# this is needed cause writeNetCDFtiles uses `name` etc

"""
    get_testcases_if_needed(pth)

Download and decompress input files if needed.
"""
function get_testcases_if_needed(pth)
    if !isdir(pth)
        run(`git clone https://github.com/gaelforget/nctiles-testcases $pth`)
        run(`gunzip $pth/diags/state_3d_set1.0000000732.data.gz`)
        run(`gunzip $pth/diags/state_3d_set1.0000001428.data.gz`)
        run(`gunzip $pth/diags/state_3d_set1.0000002172.data.gz`)
        run(`gunzip $pth/diags/trsp_3d_set1.0000000732.data.gz`)
        run(`gunzip $pth/diags/trsp_3d_set1.0000001428.data.gz`)
        run(`gunzip $pth/diags/trsp_3d_set1.0000002172.data.gz`)
    end
end

"""
    input_file_paths(inputs)

Put all folder and file paths in a dictionnary.
"""
function input_file_paths(inputs)
    diaglist = joinpath(inputs,"available_diagnostics.log")
    readme = joinpath(inputs,"README")
    grid = joinpath(inputs,"grid_float32/")
    native = joinpath(inputs,"diags/")
    interp = joinpath(inputs,"diags_interp/")
    return Dict("diaglist" => diaglist, "readme" => readme, "grid" => grid,
    "native" => native, "interp" => interp)
end

"""
    dimensions_etc_interp(pth::Dict)

Format grid variables as NCvar instances.
"""
function dimensions_etc_interp(pth::Dict)
    lon_c=-179.75:0.5:179.75
    lat_c=-89.75:0.5:89.75
    n1,n2 = (length(lon_c),length(lat_c))

    n3=50
    prec = Float32
    dep_l=-readbin(joinpath(pth["grid"],"RF.data"),prec,(n3+1,1))[2:end]
    dep_c=-readbin(joinpath(pth["grid"],"RC.data"),prec,(n3,1))[:]

    timeinterval = 3
    nsteps = length(readdir(pth["interp"] * "ETAN"))
    time_steps = timeinterval/2:timeinterval:timeinterval*nsteps
    time_units = "days since 1992-01-01 0:0:0"

    lon_c = NCvar("lon_c","degrees_east",size(lon_c),lon_c,Dict("long_name" => "longitude"),nc)
    lat_c = NCvar("lat_c","degrees_north",size(lat_c),lat_c,Dict("long_name" => "latitude"),nc)
    dep_l = NCvar("dep_l","m",size(dep_l),dep_l,Dict(["long_name" => "depth","positive"=>"down","standard_name"=>"depth"]),nc)
    dep_c = NCvar("dep_c","m",size(dep_c),dep_c,Dict(["long_name" => "depth","positive"=>"down","standard_name"=>"depth"]),nc)
    tim = NCvar("tim",time_units,Inf,time_steps,Dict(("long_name" => "time","standard_name" => "time")),nc)

    readme = readlines(pth["readme"])

    return Dict("lon_c" => lon_c, "lat_c" => lat_c, "dep_c" => dep_c,
    "dep_l" => dep_l, "tim" => tim, "readme" => readme,
    "n1" => n1, "n2" => n2, "n3" => n3)
end

"""
    dimensions_etc_native(pth::Dict)

Format grid variables as NCvar instances via MeshArrays.
"""
function dimensions_etc_native(pth::Dict)
    n3=50
    prec = Float32
    dep_l=-readbin(joinpath(pth["grid"],"RF.data"),prec,(n3+1,1))[2:end]
    dep_c=-readbin(joinpath(pth["grid"],"RC.data"),prec,(n3,1))[:]

    timeinterval = 3
    nsteps = length(filter(x -> (occursin(".data",x) && occursin("state_2d_set1",x)), readdir(pth["native"])))
    time_steps = timeinterval/2:timeinterval:timeinterval*nsteps
    time_units = "days since 1992-01-01 0:0:0"

    dep_l = NCvar("dep_l","m",size(dep_l),dep_l,Dict(["long_name" => "depth","positive"=>"down","standard_name"=>"depth"]),nc)
    dep_c = NCvar("dep_c","m",size(dep_c),dep_c,Dict(["long_name" => "depth","positive"=>"down","standard_name"=>"depth"]),nc)
    tim = NCvar("tim",time_units,Inf,time_steps,Dict(("long_name" => "time","standard_name" => "time")),nc)

    γ = GridSpec("LatLonCap",pth["grid"])
    γ = gcmgrid(pth["grid"],γ.class,γ.nFaces,γ.fSize, γ.ioSize, prec, γ.read, γ.write)
    Γ = GridLoad(γ;option="full")

    tilesize = (90,90)
    (n1,n2) = (90,1170)

    # First two dimensions
    i_c = NCvar("i_c","1",tilesize[1],1:tilesize[1],Dict("long_name" => "Cartesian coordinate 1"),nc)
    j_c = NCvar("j_c","1",tilesize[2],1:tilesize[2],Dict("long_name" => "Cartesian coordinate 2"),nc)
    i_w = NCvar("i_w","1",tilesize[1],1:tilesize[1],Dict("long_name" => "Cartesian coordinate 1"),nc)
    j_w = NCvar("j_w","1",tilesize[2],1:tilesize[2],Dict("long_name" => "Cartesian coordinate 2"),nc)
    i_s = NCvar("i_s","1",tilesize[1],1:tilesize[1],Dict("long_name" => "Cartesian coordinate 1"),nc)
    j_s = NCvar("j_s","1",tilesize[2],1:tilesize[2],Dict("long_name" => "Cartesian coordinate 2"),nc)

    # Land masks indicate which points are land, which are ocean
    land_c = Γ.hFacC
    land_w = Γ.hFacW
    land_s = Γ.hFacS
    for f in land_c.fIndex
        for d in 1:size(land_c,2)
            land_c[f,d][land_c[f,d].==0] .= NaN
            land_c[f,d][land_c[f,d].>0] .= 1

            land_w[f,d][land_w[f,d].==0] .= NaN
            land_w[f,d][land_w[f,d].>0] .= 1

            land_s[f,d][land_s[f,d].==0] .= NaN
            land_s[f,d][land_s[f,d].>0] .= 1
        end
    end
    land2d_c=land_c[:,1]

    # Variable indicating the depth/thickness of each cell
    thic_c = Γ.DRF
    thic_l = Γ.DRC[2:end]

    # TileData struct- calculates the locations of each tile in the
    # data to retrieve when needed for writing
    area_c = TileData(Γ.RAC,tilesize,γ)
    tileinfo = area_c.tileinfo; numtiles = area_c.numtiles
    area_w = TileData(Γ.RAW,tileinfo,tilesize,prec,numtiles)
    area_s = TileData(Γ.RAS,tileinfo,tilesize,prec,numtiles)
    land_c = TileData(land_c,tileinfo,tilesize,prec,numtiles)
    land2d_c = TileData(land2d_c,tileinfo,tilesize,prec,numtiles)
    land_w = TileData(land_w,tileinfo,tilesize,prec,numtiles)
    land_s = TileData(land_s,tileinfo,tilesize,prec,numtiles)
    lat_c = TileData(Γ.YC,tileinfo,tilesize,prec,numtiles)
    lon_c = TileData(Γ.XC,tileinfo,tilesize,prec,numtiles)
    lat_w = TileData(Γ.YW,tileinfo,tilesize,prec,numtiles)
    lon_w = TileData(Γ.XW,tileinfo,tilesize,prec,numtiles)
    lat_s = TileData(Γ.YS,tileinfo,tilesize,prec,numtiles)
    lon_s = TileData(Γ.XS,tileinfo,tilesize,prec,numtiles)

    # NCvar structs outline fields and their metadata to be written to the file
    lon_c = NCvar("lon","degrees_east",[i_c,j_c],lon_c,Dict("long_name" => "longitude"),nc)
    lat_c = NCvar("lat","degrees_north",[i_c,j_c],lat_c,Dict("long_name" => "latitude"),nc)
    lon_w = NCvar("lon","degrees_east",[i_w,j_w],lon_w,Dict("long_name" => "longitude"),nc)
    lat_w = NCvar("lat","degrees_north",[i_w,j_w],lat_w,Dict("long_name" => "latitude"),nc)
    lon_s = NCvar("lon","degrees_east",[i_s,j_s],lon_s,Dict("long_name" => "longitude"),nc)
    lat_s = NCvar("lat","degrees_north",[i_s,j_s],lat_s,Dict("long_name" => "latitude"),nc)
    area_c = NCvar("area","m^2",[i_c,j_c],area_c,Dict(["long_name" => "grid cell area", "standard_name" => "cell_area"]),nc)
    area_w = NCvar("area","m^2",[i_w,j_w],area_w,Dict(["long_name" => "grid cell area", "standard_name" => "cell_area"]),nc)
    area_s = NCvar("area","m^2",[i_s,j_s],area_s,Dict(["long_name" => "grid cell area", "standard_name" => "cell_area"]),nc)
    land_c = NCvar("land","1",[i_c,j_c,dep_c],land_c,Dict(["long_name" => "land mask", "standard_name" => "land_binary_mask"]),nc)
    land2d_c = NCvar("land","1",[i_c,j_c],land2d_c,Dict(["long_name" => "land mask", "standard_name" => "land_binary_mask"]),nc)
    land_w = NCvar("land","1",[i_w,j_w,dep_c],land_w,Dict(["long_name" => "land mask", "standard_name" => "land_binary_mask"]),nc)
    land_s = NCvar("land","1",[i_s,j_s,dep_c],land_s,Dict(["long_name" => "land mask", "standard_name" => "land_binary_mask"]),nc)
    thic_c = NCvar("thic","m",dep_c,thic_c,Dict("standard_name" => "cell_thickness"),nc)
    thic_l = NCvar("thic","m",dep_l,thic_l,Dict("standard_name" => "cell_thickness"),nc)

    readme = readlines(pth["readme"])

    return Dict(
    "i_c" => i_c, "i_w" => i_w, "i_s" => i_s,
    "j_c" => j_c, "j_w" => j_w, "j_s" => j_s,
    "lon_c" => lon_c, "lon_w" => lon_w, "lon_s" => lon_s,
    "lat_c" => lat_c, "lat_w" => lat_w, "lat_s" => lat_s,
    "area_c" => area_c, "area_w" => area_w, "area_s" => area_s,
    "land_c" => land_c, "land2d_c" => land2d_c,
    "land_w" => land_w, "land_s" => land_s,
    "thic_c" => thic_c, "thic_l" => thic_l,
    "dep_c" => dep_c,"dep_l" => dep_l,"tim" => tim,
    "n1" => n1, "n2" => n2, "n3" => n3, "γ" => γ,
    "tilesize" => tilesize, "numtiles" => numtiles, "readme" => readme)
end

"""
    prep_nctiles_interp(inputs,set,nam,prc)

Setup NCvar instance for chosen variable that can then be written to file.
"""
function prep_nctiles_interp(inputs,set,nam,prc)
    pth=input_file_paths(inputs)
    Λ = dimensions_etc_interp(pth) #dimensions, sizes, and meta data

    flddatadir = joinpath(pth["interp"],nam)
    fnames = joinpath.(Ref(flddatadir),filter(x -> occursin(".data",x), readdir(flddatadir)))
    diaginfo = read_available_diagnostics(nam, filename=pth["diaglist"])
    if diaginfo["levs"]==1 && diaginfo["code"][2]=='M'
        flddata = BinData(fnames,prc,(Λ["n1"],Λ["n2"]))
        dims = [Λ["lon_c"],Λ["lat_c"],Λ["tim"]]
    elseif diaginfo["levs"]>1 && diaginfo["code"][2]=='M' && diaginfo["code"][9]=='L'
        flddata = BinData(fnames,prc,(Λ["n1"],Λ["n2"],Λ["n3"]))
        dims = [Λ["lon_c"],Λ["lat_c"],Λ["dep_l"],Λ["tim"]]
    else
        error("other cases remain to be implemented...")
    end
    field = NCvar(nam,diaginfo["units"],dims,flddata,
    Dict("long_name" => diaginfo["title"]),nc)
    savename=joinpath(writedir,nam*".nc")
    readme=Λ["readme"]

    #nct=...
    return field,savename,readme
end

"""
    prep_nctiles_native(inputs,set,nam,prc)

Setup NCvar array for chosen variable that can then be written to file. The
    NCvar array includes ancillary variables like lon and lat.
"""
function prep_nctiles_native(inputs,set,nam,prc)
    pth=input_file_paths(inputs)
    Λ = dimensions_etc_native(pth) #dimensions, sizes, and meta data

    fnames = pth["native"]*'/'.*filter(x -> (occursin(".data",x) && occursin(set,x)), readdir(pth["native"]))
    savepath = joinpath(writedir,nam)
    if ~ispath(savepath); mkpath(savepath); end
    savename = joinpath.(Ref(savepath),nam)
    diaginfo = read_available_diagnostics(nam, filename=pth["diaglist"])
    if diaginfo["levs"]==1 && diaginfo["code"][2]=='M'
        flddata = BinData(fnames,prc,(Λ["n1"],Λ["n2"]))
        dims = [Λ["i_c"],Λ["j_c"],Λ["tim"]]
        coords = join(replace([dim.name for dim in dims],"i_c" => "lon", "j_c" => "lat")," ")
        tmp=Dict("lon" => Λ["lon_c"],"lat" => Λ["lat_c"],
        "area" => Λ["area_c"],"land" => Λ["land2d_c"])
    elseif diaginfo["levs"]>1 && diaginfo["code"][2]=='M' && diaginfo["code"][9]=='M'
        flddata = BinData(fnames,prc,(Λ["n1"],Λ["n2"],Λ["n3"]))
        dims = [Λ["i_c"],Λ["j_c"],Λ["dep_c"],Λ["tim"]]
        coords = join(replace([dim.name for dim in dims],"i_c" => "lon", "j_c" => "lat")," ")
        tmp=Dict(    "lon" => Λ["lon_c"],"lat" => Λ["lat_c"],
        "area" => Λ["area_c"],"land" => Λ["land_c"],"thic" => Λ["thic_c"])
    elseif diaginfo["levs"]>1 && diaginfo["code"][2]=='U' && diaginfo["code"][9]=='M'
        flddata = BinData(fnames,prc,(Λ["n1"],Λ["n2"],Λ["n3"]))
        dims = [Λ["i_w"],Λ["j_w"],Λ["dep_c"],Λ["tim"]]
        coords = join(replace([dim.name for dim in dims],"i_w" => "lon", "j_w" => "lat")," ")
        tmp=Dict(    "lon" => Λ["lon_w"],"lat" => Λ["lat_w"],
        "area" => Λ["area_w"],"land" => Λ["land_w"],"thic" => Λ["thic_c"])
    elseif diaginfo["levs"]>1 && diaginfo["code"][2]=='V' && diaginfo["code"][9]=='M'
        flddata = BinData(fnames,prc,(Λ["n1"],Λ["n2"],Λ["n3"]))
        dims = [Λ["i_s"],Λ["j_s"],Λ["dep_c"],Λ["tim"]]
        coords = join(replace([dim.name for dim in dims],"i_s" => "lon", "j_s" => "lat")," ")
        tmp=Dict(    "lon" => Λ["lon_s"],"lat" => Λ["lat_s"],
        "area" => Λ["area_s"],"land" => Λ["land_s"],"thic" => Λ["thic_c"])
    elseif diaginfo["levs"]>1 && diaginfo["code"][2]=='M' && diaginfo["code"][9]=='L'
        flddata = BinData(fnames,prc,(Λ["n1"],Λ["n2"],Λ["n3"]))
        dims = [Λ["i_c"],Λ["j_c"],Λ["dep_l"],Λ["tim"]]
        coords = join(replace([dim.name for dim in dims],"i_c" => "lon", "j_c" => "lat")," ")
        tmp=Dict(    "lon" => Λ["lon_c"],"lat" => Λ["lat_c"],
        "area" => Λ["area_c"],"land" => Λ["land_c"],"thic" => Λ["thic_l"])
    else
        error("other cases remain to be implemented...")
    end

    tilfld = TileData(flddata,Λ["tilesize"],Λ["γ"])
    flds = Dict([nam => NCvar(nam,diaginfo["units"],dims,tilfld,Dict(["long_name" => diaginfo["title"], "coordinates" => coords]),nc)])
    merge!(flds,tmp)

    #nct=...
    return flds,savename,readme
end

max_stars_repo_name: gaelforget/Notebooks

max_issues_repo_name: gaelforget/Notebooks

max_forks_repo_name: gaelforget/Notebooks


text:
module LogisticRegression

using DataModel
using StatsBase

export LogisticRegressionModel
"""
Datamodel for logistic regression

Fields
    - d: dimension
    - x: covariantes N x d
    - y: labels (±1)
    - priorPrec: prior precision (normal prior)

Constructors
    `LogisticRegressionModel(d,x,y,priorPrec;intercept=false, intercept_priorprec = 1.0)`
    `intercept = true` extend x to include an intercept term and adds
    `intercept_priorprec` to the prior precision matrix
"""
type LogisticRegressionModel <: AbstractDataModel
    d::Int64
    x::Array{Float64,2}
    y::Array{Float64}
    priorPrec::Array{Float64,2}
    true_weights::Nullable{Array{Float64,1}}
    function LogisticRegressionModel(x,y,priorPrec::Array{Float64,1};intercept=false, intercept_priorprec= 1.0)
        priorPrec = diagm(priorPrec)
        d = size(x)[2]
        N = size(x)[1]
        @assert N == length(y)
        @assert size(priorPrec)[1] == size(priorPrec)[2] == d
        @assert isposdef(priorPrec)
        if intercept
            d += 1
            x = hcat(ones(N),x)
            priorExtended = zeros(d,d)
            priorExtended[1,1] = intercept_priorprec
            priorExtended[2:end,2:end] = priorPrec
            priorPrec = priorExtended
        end
        new(d,x,y,priorPrec,Nullable{Array{Float64,1}}())
    end
    function LogisticRegressionModel(x,y,priorPrec::Array{Float64,2};intercept=false, intercept_priorprec= 1.0)
        d = size(x)[2]
        N = size(x)[1]
        @assert N == length(y)
        @assert size(priorPrec)[1] == size(priorPrec)[2] == d
        @assert isposdef(priorPrec)
        @assert intercept_priorprec > 0.
        if intercept
            d += 1
            x = hcat(ones(N),x)
            priorExtended = zeros(d,d)
            priorExtended[1,1] = intercept_priorprec
            priorExtended[2:end,2:end] = priorPrec
            priorPrec = priorExtended
        end
        new(d,x,y,priorPrec,Nullable{Array{Float64,1}}())
    end
    function LogisticRegressionModel(x,y,priorPrec,true_weights::Array{Float64,1};intercept=false,intercept_priorprec= 1.0)
        @assert size(x)[1] == length(y) == size(priorPrec)[1]
        d = size(x)[2]
        N = size(x)[1]
        @assert isposdef(priorPrec)
        @assert intercept_priorprec > 0.
        if intercept
            d += 1
            x = hcat(ones(N),x)
            precExtended = zeros(d,d)
            priorExtended[1,1] = intercept_priorprec
            priorExtended[2:end,2:end] = priorPrec
        end
        new(d,x,y,priorPrec,Nullable(vcat(0,true_weights)))
    end
end

#-------------------------
# Helper functions
#-------------------------
logit(z) = 1.0./(1.0.+exp(-z))
log_logit(z)= -log(1.0 .+ exp(-z))
grad_log_logit(z)=1.0-logit(z)
function ddlogit(x)
    return exp(-x)/(1+exp(-x))^2
end

"""

"""
function DataModel.getN(dm::LogisticRegressionModel)
  return length(dm.y)
end

"""
`DataModel.getllik(dm::LogisticRegressionModel;nobs=length(dm.y))`
returns a full loglikelihood function. The optional `nobs` keyword argument restrict the dataset to the first `nobs` data points.
"""
function DataModel.getllik(dm::LogisticRegressionModel;nobs::Int64=length(dm.y))
    function logdensity(beta)
        d=dm.d
        x=dm.x
        y=dm.y
        priorPrec=dm.priorPrec
        log_prior = -0.5 * dot(beta,priorPrec*beta)
        log_like= sum(log_logit((x[1:nobs,:]* beta).*y[1:nobs]))
        return log_prior+log_like
    end
    return logdensity
end

"""
`DataModel.getgrad(dm::LogisticRegressionModel,batchsize;nobs=length(dm.y))`
returns a stochastic gradient function using batches of size `batchsize`. The optional `nobs` keyword argument restrict the dataset to the first `nobs` data points.
"""
function DataModel.getgrad(dm::LogisticRegressionModel,batchsize::Int64;nobs::Int64=length(dm.y))
    function grad_log_posterior_sub(beta)
        d=dm.d
        x=dm.x
        y=dm.y
        priorPrec=dm.priorPrec
        chosen_indices=sample(1:nobs,batchsize,replace=false)
        log_prior_gradient= -priorPrec* beta

        #first step: compute y*grad_log_logit(y*beta*x)
        weights = y[chosen_indices].*grad_log_logit((x[chosen_indices,:] * beta ).*y[chosen_indices])
        #second ztep: compute y*grad_log_logit(y*beta*x)*x

        return log_prior_gradient+(1.0*nobs)/batchsize*(reshape(weights,1,batchsize)* x[chosen_indices,:])[:]
    end
    return grad_log_posterior_sub
end
"""
`getgrad_i(dm::LogisticRegressionModel,batchsize;nobs=length(dm.y))`
returns a stochastic gradient function where index can be specified in array
"""
function getgrad_i(dm::LogisticRegressionModel;nobs=length(dm.y))
    function grad_log_posterior_sub(beta::Array{Float64,1},chosen_indices)
        d=dm.d
        x=dm.x
        y=dm.y
        priorPrec=dm.priorPrec
        batchsize=length(chosen_indices)
        log_prior_gradient= -priorPrec* beta

        #first step: compute y*grad_log_logit(y*beta*x)
        weights = y[chosen_indices].*grad_log_logit((x[chosen_indices,:] * beta ).*y[chosen_indices])
        #second ztep: compute y*grad_log_logit(y*beta*x)*x

        return log_prior_gradient+(1.0*nobs)/batchsize*(reshape(weights,1,batchsize)* x[chosen_indices,:])[:]
    end
    return grad_log_posterior_sub
end

"""
`DataModel.getfullgrad(dm::LogisticRegressionModel;nobs=length(dm.y))`
returns the full gradient function. The optional `nobs` keyword argument restrict the dataset to the first `nobs` data points.
"""
function DataModel.getfullgrad(dm::LogisticRegressionModel;nobs=length(dm.y))

    function grad_log_posterior(beta)::Array{Float64,1}
        d=dm.d
        x=dm.x
      y=dm.y
        priorPrec=dm.priorPrec
        chosen_indices=vec(1:nobs)
        log_prior_gradient= -priorPrec* beta

        #first step: compute y*grad_log_logit(y*beta*x)
        weights = y[chosen_indices].*grad_log_logit((x[chosen_indices,:] * beta ).*y[chosen_indices])
        #second ztep: compute y*grad_log_logit(y*beta*x)*x

        return log_prior_gradient+(reshape(weights,1,nobs)* x[chosen_indices,:])[:]
    end
    return grad_log_posterior
end
"""
`DataModel.getfullhess(dm::LogisticRegressionModel;nobs=length(dm.y))`
returns the full hessian function. The optional `nobs` keyword argument restrict the dataset to the first `nobs` data points.
"""
function DataModel.getfullhess(dm::LogisticRegressionModel;nobs=length(dm.y))

    d=dm.d
    x=dm.x
    y=dm.y
    Cinv=dm.priorPrec
    function hessfun(beta_vec)
      hess=zeros(d,d)
      hess+=Cinv;
      for i=1:(nobs)
        hess+=ddlogit(y[i]*dot(x[i,:],beta_vec)[1])*y[i]^2*x[i,:]*x[i,:]';
      end
      hess+=Cinv;
      return hess
    end
    return hessfun;
end

"""
`artificial_logreg_dm(nObs::Int64,dim::Int64;seed::Int64 =1, diagonal_covariance::Bool = false, intercept::Bool = false)` creates an artificial dataset which samples `nObs` observations from the logistic regression model. The generating weights are saved in `dm.true_weights`.
"""
function artificial_logreg_dm(nObs::Int64,dim::Int64;seed::Int64 =1, diagonal_covariance::Bool = false, intercept::Bool = false)
    srand(seed)
    # output: (X,Y) [(NxD), (N,1)], observations and labels
    # .
    # . covariance (size dxd)
    P   = diagonal_covariance ? diagm(-1+2*rand(dim)):(-1+2*rand(dim,dim))
    cov = P * P'
    # . mean (size dx1), uniformly distributed on [0,1]
    mu = rand(dim)
    # . observations (size NxD)
    X = (repmat(mu,1,nObs)+P*randn(dim,nObs))'
    if intercept
        X = hcat(ones(nObs),X)
    end
    # . true parameters (size dx1)
    priorCov = 10*eye(dim)
    priorPrec = 0.1*eye(dim)
    C = copy(priorPrec)
    if intercept
        w = vcat(randn(1),chol(priorCov)'*randn(dim))
        dim += 1
        priorPrec= eye(dim)
        priorPrec[2:end,2:end] = C
    else
        w = chol(priorCov)'*randn(dim)
    end
    # . generate labels (size Nx1)
    y01      = rand(nObs).<logit(X*w)
    y        = y01*2.-1.
    dm = LogisticRegressionModel(X,y,priorPrec)
    dm.true_weights = Nullable(w)
    srand()
    dm
end
end


text:
module LyapunovExponentsTests

printstyled(let
    message = " Warnings below (if any) are fine. "
    margin = (displaysize(stdout)[2] - length(message)) ÷ 2
    ("=" ^ margin) * message * ("=" ^ margin)
end, color=:blue)
println()
flush(stdout)
import Plots
import ForwardDiff
import DifferentialEquations
import OnlineStats
flush(stdout)
flush(stderr)
printstyled("=" ^ displaysize(stdout)[2], color=:blue)
println()

using LyapunovExponents
using Test

@time begin
include("test_testtools.jl")
include("test_utils.jl")
include("test_online_stats.jl")
include("test_smoke.jl")
include("test_ui.jl")
include("test_examples.jl")
include("test_integrators.jl")
include("test_terminators.jl")
include("test_clv.jl")
include("test_null_clv.jl")
end

end


text:
module JuliaImageBuilder

export build

using Pkg
using PackageCompiler

include("type.jl")
include("config.jl")
include("builder.jl")

end


text:
#!/usr/bin/env julia

using Colors
using Compose
using Base.MathConstants

function golden_rect(n::Int)
    poly_points = [(0, 0), (1, 0), (1, 1), (0, 1)]
    if n == 0 return context() end
    c = compose(context(), polygon(poly_points), fill(LCHab(90, 80, 70-20n)), stroke("black"))
    compose(c, (context(0,0,1/φ,1/φ, rotation=Rotation(-π/2,1,0)),  golden_rect(n-1)))
end

draw(SVG("golden_rect.svg", φ*3inch, 3inch),
    compose(golden_rect(8), linewidth(0.2mm)))


text:
struct WeekdayHolidayObsModel{T,M,A,P,E,W,H,V} <: CovidSurvey.ObservationsModel
    ps::T
    μ::M
    α::A
    ϕ::P
	expecteds::E
    weekdayeffect::W
    holidayeffect::H
	vectorize::V
end

struct SimpleObsModel{T,M,A,P,E,V} <: CovidSurvey.ObservationsModel
    ps::T
    μ::M
    α::A
    ϕ::P
	expecteds::E
	vectorize::V
end

struct SimpleObsModel2{T,M,S,A,P,E,V} <: CovidSurvey.ObservationsModel
    ps::T
    μ::M
	σ::S
    α::A
    ϕ::P
	expecteds::E
	vectorize::V
end

rescale(xs, σ, μs) = @. xs * σ + μs'
# ============================================================================
# random walk model
function random_walks!(Rts, θ, predict, latent_Rts::AbstractVector, R0s, σ_rt)
	@unpack num_regions, rt_step_indices, lockdown_indices, num_observations, link, include_early_dynamic = θ #include_early_dynamic

	for m in 1:num_regions
		rt_step_index = rt_step_indices[m]
		Rt            = Rts[m]
		lockdown      = lockdown_indices[m]
		num_obs       = num_observations[m]
		latent_Rt     = latent_Rts[m]#latent_Rts[:,m] #include_early_dynamic ? latent_Rts[m] : latent_Rts[:,m]
		R0            = R0s[m]

		if include_early_dynamic
			Rt[1:num_obs] = link.(latent_Rt[rt_step_index])
		else
			Rt[1:lockdown] .= R0
			Rt[lockdown+1:num_obs] = link.(latent_Rt[rt_step_index])
		end
		if predict
			Rt[num_obs+1:end] .= Rt[num_obs]
		end
	end
end

function random_walks!(Rts, θ, predict, latent_Rts::AbstractMatrix, R0s, σ_rt)
	@unpack num_regions, rt_step_indices, lockdown_indices, num_observations, link, include_early_dynamic = θ

	for m in 1:num_regions
		rt_step_index = rt_step_indices[m]
		Rt            = Rts[m]
		lockdown      = lockdown_indices[m]
		num_obs       = num_observations[m]
		latent_Rt     = latent_Rts[:,m] #include_early_dynamic ? latent_Rts[m] : latent_Rts[:,m]
		R0            = R0s[m]

		if include_early_dynamic
			Rt[1:num_obs] = link.(latent_Rt[rt_step_index])
		else
			Rt[1:lockdown] .= R0
			Rt[lockdown+1:num_obs] = link.(latent_Rt[rt_step_index])
		end
		if predict
			Rt[num_obs+1:end] .= Rt[num_obs]
		end
	end
end

# function random_walks!(Rts, θ, predict, latent_Rts, R0s)
# 	@unpack num_regions, rt_step_indices, lockdown_indices, num_observations, link = θ
#
# 	for m in 1:num_regions
# 		rt_step_index = rt_step_indices[m]
# 		Rt            = Rts[m]
# 		lockdown      = lockdown_indices[m]
# 		num_obs       = num_observations[m]
# 		latent_Rt     = latent_Rts[:,m]
# 		R0            = R0s[m]
#
# 		Rt[1:lockdown] .= R0
# 		Rt[lockdown+1:num_obs] = link.(latent_Rt[rt_step_index])
# 		if predict
# 			Rt[num_obs+1:end] .= Rt[num_obs]
# 		end
# 	end
# end

# ============================================================================
# semi-parametric model
function semiparametric!(Rts, θ, predict, latent_Rts, R0s, σ_rt, effects)
	@unpack num_regions, lockdown_indices, num_observations, link, rt_step_indices, covariates, covariates_starts, num_total_days, include_early_dynamic= θ #rt_step_indices

	for m in 1:num_regions
		Rt            = Rts[m]
		lockdown      = lockdown_indices[m]
		rt_step_index = rt_step_indices[m]
		num_obs       = num_observations[m]
		latent_Rt     = include_early_dynamic ? latent_Rts[m] : latent_Rts[:,m] #latent_Rts[m]
		R0            = R0s[m]
		num_total     = num_total_days[m]
		covariate     = covariates[m]
		effect        = effects[m]
		covariates_start = covariates_starts[m]

		if include_early_dynamic
			Rt[1:covariates_start-1] = link.(latent_Rt[rt_step_index])
		else
			Rt[1:lockdown] .= R0
			Rt[lockdown+1:covariates_start-1] = link.(latent_Rt[rt_step_index])
		end
		Rt[covariates_start:num_obs] = link.(
			latent_Rt[end] .+ covariate * effect
		)
		if predict
			Rt[num_obs+1:num_total] .= Rt[num_obs]
			# @unpack num_total_days, invlink = θ
			# n  = num_total_days - num_observations
			# n < 1 && return nothing
			# rw = RandomWalk(n, σ_rt, invlink(Rt[num_observations]))
			# Rt[num_observations+1:end] = link.( rand( rw ) )
		end
	end
	return nothing
end

# ============================================================================
# semi-parametric model
function mixed!(Rts, θ, predict, latent_Rts, R0s, σ_rt, effects)
	@unpack (
		num_regions,
		lockdown_indices,
		num_observations,
		link,
		rt_step_indices,
		covariates,
		covariates_starts,
		num_total_days,
		include_early_dynamic
	) = θ #rt_step_indices

	for m in 1:num_regions
		Rt               = Rts[m]
		lockdown_index   = lockdown_indices[m]
		rt_step_index    = rt_step_indices[m]
		num_obs          = num_observations[m]
		latent_Rt        = include_early_dynamic ? latent_Rts[m] : latent_Rts[:,m]#latent_Rts[m]
		R0               = R0s[m]
		num_total        = num_total_days[m]
		covariate        = covariates[m]
		effect           = effects[m]
		covariates_start = covariates_starts[m]
		predictor        = covariate * effect

		if include_early_dynamic
			for i in 1:covariates_start-1
				Rt[i] = link(latent_Rt[rt_step_index[i]])
			end
			for (j,i) in enumerate(covariates_start:num_obs)
				Rt[i] = link(latent_Rt[rt_step_index[i]] + predictor[j])
			end
		else
			Rt[1:lockdown_index] .= R0
			for (j,i) in enumerate(lockdown_index+1:covariates_start-1)
				Rt[i] = link(latent_Rt[rt_step_index[j]])
			end
			n = covariates_start - lockdown_index - 1
			for (j,i) in enumerate(covariates_start:num_obs)
				Rt[i] = link(latent_Rt[rt_step_index[j+n]] + predictor[j])
			end
		end
		if predict
			Rt[num_obs+1:num_total] .= Rt[num_obs]
			# @unpack num_total_days, invlink = θ
			# n  = num_total_days - num_observations
			# n < 1 && return nothing
			# rw = RandomWalk(n, σ_rt, invlink(Rt[num_observations]))
			# Rt[num_observations+1:end] = link.( rand( rw ) )
		end
	end
	return nothing
end

# ============================================================================
# infection model
function infections!(newly_infecteds, cumulative_infecteds, θ, ys, Rts)
	@unpack num_regions, populations = θ
	for m in 1:num_regions
		regionaldata = (
			Rt                     = Rts[m],
			population             = populations[m],
			y                      = ys[m]
		)

		initepidemic!(newly_infecteds[m], cumulative_infecteds[m], θ, regionaldata)
		runepidemic!(newly_infecteds[m], cumulative_infecteds[m], θ, regionaldata)
	end
end

function initepidemic!(newly_infected, cumulative_infected, effective_Rt, θ, regionaldata)
	@unpack num_impute = θ
	@unpack y, Rt, population = regionaldata
	newly_infected[1]      = y
	cumulative_infected[1] = zero( eltype(cumulative_infected) )
	effective_Rt[1]        = Rt[1]
	for t in 2:num_impute
		newly_infected[t] = y
		cumulative_infected[t] = cumulative_infected[t-1] + y
		St = max(population - cumulative_infected[t], 0) / population
		effective_Rt[t] = St * Rt[t]
	end
end

function initepidemic!(newly_infected, cumulative_infected, θ, regionaldata)
	@unpack num_impute = θ
	@unpack y, population = regionaldata
	newly_infected[1]      = y
	cumulative_infected[1] = zero( eltype(cumulative_infected) )
	for t in 2:num_impute
		newly_infected[t] = y
		cumulative_infected[t] = cumulative_infected[t-1] + y
		# St = max(population - cumulative_infected[t], 0) / population
	end
end

function runepidemic!(newly_infected, cumulative_infected, effective_Rts, θ, regionaldata)
	@unpack num_impute, serial_interval, num_si = θ
	@unpack population, Rt = regionaldata
	num_time_step = length(newly_infected)

	for t = (num_impute + 1):num_time_step
		# Update cumulatively infected
		C = cumulative_infected[t] = cumulative_infected[t-1] + newly_infected[t - 1]
		susceptible = max(population - C, 0)

		# Adjusts for portion of population that is susceptible
		#St = max(population - cumulative_infected[t], 0) / population
		# effective number of infectious individuals
		effectively_infectious = sum(newly_infected[τ] * serial_interval[t - τ] for τ = (t - 1):-1:max(t-num_si,1))

		unadjusted_new_infections = effectively_infectious * Rt[t]

		# number of new infections (unobserved)
		# newly_infected[t] = susceptible * (1 - exp(- unadjusted_new_infections / population))
		newly_infected[t] = susceptible * (1 - exp(- unadjusted_new_infections / population))
		# newly_infected[t] = effective_Rt[t] * effectively_infectious

		effective_Rt[t] = susceptible / population * Rt[t]
	end
end

function runepidemic!(newly_infected, cumulative_infected, θ, regionaldata)
	@unpack num_impute, serial_interval, num_si = θ
	@unpack population, Rt = regionaldata
	num_time_step = length(newly_infected)

	for t = (num_impute + 1):num_time_step
		# Update cumulatively infected
		C = cumulative_infected[t] = cumulative_infected[t-1] + newly_infected[t - 1]
		susceptible = max(population - C, 0)

		# Adjusts for portion of population that is susceptible
		effectively_infectious = sum(newly_infected[τ] * serial_interval[t - τ] for τ = (t - 1):-1:max(t-num_si,1))

		unadjusted_new_infections = effectively_infectious * Rt[t]

		# number of new infections (unobserved)
		# newly_infected[t] = susceptible * (1 - exp(- unadjusted_new_infections / population))
		newly_infected[t] = susceptible * (1 - exp(- unadjusted_new_infections / population))

	end
	return nothing
end

# ============================================================================
# unbound infection model
function infections!(newly_infecteds, θ, ys, Rts)
	@unpack num_regions = θ
	for m in 1:num_regions
		regionaldata = (
			Rt                     = Rts[m],
			y                      = ys[m]
		)

		initepidemic!(newly_infecteds[m], θ, regionaldata)
		runepidemic!(newly_infecteds[m], θ, regionaldata)
	end
end

function initepidemic!(newly_infected, θ, regionaldata)
	@unpack num_impute = θ
	@unpack y, Rt = regionaldata
	newly_infected[1] = y
	for t in 2:num_impute
		newly_infected[t] = y
	end
end

function runepidemic!(newly_infected, θ, regionaldata)
	@unpack num_impute, serial_interval, num_si = θ
	@unpack Rt = regionaldata
	num_time_step = length(newly_infected)

	for t = (num_impute + 1):num_time_step
		effectively_infectious = sum(newly_infected[τ] * serial_interval[t - τ] for τ = (t - 1):-1:max(t-num_si,1))
		newly_infected[t] = effectively_infectious * Rt[t]
	end
	return nothing
end

# ============================================================================
# observations

function infection2observation(obsmodel::SimpleObsModel)
	@unpack ps, μ = obsmodel
	@unpack delay_dispersion, delay_length, delay_dist = ps
	i2o = pdf.( Ref(delay_dist(μ, delay_dispersion)), 1:delay_length )
	i2o /= sum(i2o)
	return i2o
end

function infection2observation(obsmodel::SimpleObsModel2)
	@unpack ps, μ, σ = obsmodel
	@unpack  delay_length, delay_dist = ps
	i2o = pdf.( Ref(delay_dist(μ, σ)), 1:delay_length )
	i2o /= sum(i2o)
	return i2o
end

function expected!(obsmodel::Union{SimpleObsModel, SimpleObsModel2}, newly_infecteds)
	@unpack ps, α, expecteds = obsmodel
	@unpack num_regions, delay_length = ps

	i2o = infection2observation(obsmodel)

	for m in 1:num_regions

		expected       = expecteds[m]
		newly_infected = newly_infecteds[m]
		num_time_step  = length(expected)

		expected[1] = 1e-15*newly_infected[1]
		for t = 2:num_time_step
			expected[t] = α * sum(newly_infected[τ] * i2o[t - τ] for τ = (t - 1):-1:max(t-delay_length,1))
			# expected[t] = α * sum(newly_infected[τ] * i2o[t - τ] for τ = (t - 1):-1:max(t-delay_length,1))
		end
	end
	return nothing
end

function expected!(obsmodel::WeekdayHolidayObsModel, newly_infecteds)
	@unpack ps, α, expecteds, holidayeffect, weekdayeffect = obsmodel
	@unpack num_regions, delay_length, holidays, weekdays = ps

	i2o = infection2observation(obsmodel)

	for m in 1:num_regions

		expected       = expecteds[m]
		newly_infected = newly_infecteds[m]
		num_time_step  = length(expected)
		holiday        = holidays[m]
		weekday        = weekdays[m]

		expected[1] = 1e-15*newly_infected[1]
		for t = 2:num_time_step
			weekday_holiday_effect = (1 - holidayeffect * holiday[t]) * weekdayeffect[weekday[t]]
			expected[t] = α * weekday_holiday_effect * sum(newly_infected[τ] * i2o[t - τ] for τ = (t - 1):-1:max(t-delay_length,1))
		end
	end
	return nothing
end

function _fastinference(ys, μs, ϕ)
	dist = arraydist(NegativeBinomial2.(μs, Ref(ϕ)))
	return logpdf(dist, ys)
end

function _slowinference(ys, μs, ϕ)
	ℓ = zero(typeof(ϕ))
	for (μ,y) in zip(μs,ys)
		dist = NegativeBinomial2(μ, ϕ)
		ℓ   += logpdf(dist, y)
	end
	return ℓ
end

function _prediction(ys, μs, ϕ)
	ℓ = zero(typeof(ϕ))
	for μ in μs
		dist = NegativeBinomial2(μ, ϕ)
		y    = rand(dist)
		ℓ   += logpdf(dist, y)
	end
	return ℓ
end

function Turing.logpdf(obsmodel::CovidSurvey.ObservationsModel, observed)
	@unpack ps, ϕ, expecteds, vectorize = obsmodel
	@unpack populations, starts, stops, num_regions = ps

	# _logpdf = if vectorize
	# 	_fastinference
	# else
	# 	if eltype(ϕ) <: Real
	# 		_prediction
	# 	else
	# 		_slowinference
	# 	end
	# end
	_logpdf = _fastinference

	T = typeof(ϕ)
	ℓ = zero(T)
	for m in 1:num_regions
		ts    = starts[m]:stops[m]
		μs = expecteds[m][ts]
		ys = observed[m][ts]

		#_out_of_bounds(μs, populations[m]) && (@warn "out of bounds"; return T(Inf))
		ℓ += _logpdf(ys, μs, ϕ)
	end
	return ℓ
end

_out_of_bounds(x, upper_lim) = !all(@. 0. < x < upper_lim )

# ============================================================================
# pointwise_loglikelihoods
function Turing.pointwise_loglikelihoods(obsmodel::CovidSurvey.ObservationsModel, observed)
	@unpack ps, ϕ, expecteds, vectorize = obsmodel
	@unpack populations, starts, stops, num_regions = ps

	logliks = Vector{Float64}[]
	for m in 1:num_regions
		ts    = starts[m]:stops[m]
		μs = expecteds[m][ts]
		ys = observed[m][ts]

		loglik = [ logpdf(NegativeBinomial2(μ, ϕ), y) for (μ,y) in zip(μs,ys) ]
		push!(logliks, loglik)
	end
	return logliks
end

# ============================================================================
# predictions
function prediction(obsmodel::CovidSurvey.ObservationsModel)
	@unpack ps, ϕ, expecteds = obsmodel
	@unpack starts, stops, num_regions = ps

	ys = Vector{Float64}[]
	for m in 1:num_regions
		μs = expecteds[m]
		μs[μs .< 1e-2] .= 1e-2
		y = [ rand(NegativeBinomial2(μ, ϕ)) for μ in μs ]
		push!(ys, y)
	end
	return ys
end

# ============================================================================
# weekday effect

function weekdayeffect!(weekdayeffect, weekdayeffect_simplex)
	num_days = 7
	for i in eachindex(weekdayeffect)
		@inbounds weekdayeffect[i] = weekdayeffect_simplex[i] * num_days
	end
end

# ============================================================================
# observe sero positive

struct SimpleSeroObsModel{T,E}# <: CovidSurvey.ObservationsModel
    ps::T
	expecteds::E
end

_unrealistic(ŷ,y) = any( ŷ .> 20*y )
# _unrealistic(y,ŷ) = any( ŷ .> 1.5*y )
function Turing.logpdf(obsmodel::SimpleSeroObsModel, observeds)
	@unpack ps, expecteds = obsmodel
	@unpack indexs, delay, populations, std, num_regions = ps

	T = eltype(expecteds[1])
	ℓ = zero(T)
	for m in 1:num_regions
		index      = indexs[m]
		expected   = expecteds[m]
		population = populations[m]
		observed   = observeds[m]

		μ    = expected[index .- delay] / population
		dist = MvNormal(μ, std) # select only early observations
		#_unrealistic(μ, observed) && (return T(Inf))
		ℓ += logpdf(dist, observed)
	end
	return ℓ
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
# ? FIXME: Can this be merged with CorrelationFunctions.jl somehow?

# 2D
get_slice(a         :: AbstractArray{T, N},
          periodic  :: Bool,
          idx       :: NTuple{N, Int},
          direction :: Symbol) where {T, N} =
              get_slice(a, periodic, idx, Val(direction))

# Helper function for diagonals
function get_slice(a        :: AbstractArray{T, N},
                   periodic :: Bool,
                   idx      :: NTuple{N, Int},
                   forward  :: NTuple{N, Bool}) where {T, N}
    s = size(a)
    directions = (f ? 1 : -1 for f in forward)
    a = CircularArray(a)

    if periodic
        m = (length(a) == 3) ? idx[3] : idx[2]
        start = (f ? x - m + 1 : x + m - 1 for (f, x) in zip(forward, idx))
        indices = take(zip((countfrom(s, inc) for (s, inc) in zip(start, directions))...), s[1])
    else
        m = minimum(f ? x : s - x + 1 for (f, x, s) in zip(forward, idx, s))
        start = (f ? x - m + 1 : x + m - 1 for (f, x) in zip(forward, idx))
        indices = takewhile(x -> checkbounds(Bool, a, x...),
                            zip((countfrom(s, inc) for (s, inc) in zip(start, directions))...))
    end
    return [@inbounds a[index...] for index in indices]
end

macro def_slicer(ndims, direction, expr)
    :(get_slice(a        :: AbstractArray{T, $ndims},
                periodic :: Bool,
                idx      :: NTuple{$ndims, Int},
                _        :: Val{$direction}) where T = $expr)
end

# 2D
@def_slicer 2 :x let (_, y) = idx; a[:, y] end
@def_slicer 2 :y let (x, _) = idx; a[x, :] end
@def_slicer 2 :xy get_slice(a, periodic, idx, (true, true))
@def_slicer 2 :yx get_slice(a, periodic, idx, (false, true))

# 3D
@def_slicer 3 :x let (_, y, z) = idx; a[:, y, z] end
@def_slicer 3 :y let (x, _, z) = idx; a[x, :, z] end
@def_slicer 3 :z let (x, y, _) = idx; a[x, y, :] end

@def_slicer 3 :xy let (x, y, z) = idx; get_slice(a[:,:,z], periodic, (x, y), :xy) end
@def_slicer 3 :yx let (x, y, z) = idx; get_slice(a[:,:,z], periodic, (x, y), :yx) end
@def_slicer 3 :xz let (x, y, z) = idx; get_slice(a[:,y,:], periodic, (x, z), :xy) end
@def_slicer 3 :zx let (x, y, z) = idx; get_slice(a[:,y,:], periodic, (x, z), :yx) end
@def_slicer 3 :yz let (x, y, z) = idx; get_slice(a[x,:,:], periodic, (y, z), :xy) end
@def_slicer 3 :zy let (x, y, z) = idx; get_slice(a[x,:,:], periodic, (y, z), :yx) end

@def_slicer 3 :xyz get_slice(a, periodic, idx, (true,  true,  true))
@def_slicer 3 :yxz get_slice(a, periodic, idx, (false, true,  true))
@def_slicer 3 :xzy get_slice(a, periodic, idx, (true,  false, true))
@def_slicer 3 :zyx get_slice(a, periodic, idx, (true,  true,  false))

macro def_same_line(ndims, direction, expr)
    :(same_line_p(idx1 :: CartesianIndex{$ndims},
                  idx2 :: CartesianIndex{$ndims},
                  _    :: Val{$direction}) = $expr)
end

same_line_p(idx1      :: CartesianIndex{N},
            idx2      :: CartesianIndex{N},
            direction :: Symbol) where N =
                same_line_p(idx1, idx2, Val(direction))

@def_same_line 2 :x  idx1[2] == idx2[2]
@def_same_line 2 :y  idx1[1] == idx2[1]
@def_same_line 2 :xy (idx2[2] - idx1[2]) == (idx2[1] - idx1[1])
@def_same_line 2 :yx (idx2[2] - idx1[2]) == (idx1[1] - idx2[1])

@def_same_line 3 :x  idx1[2] == idx2[2] && idx1[3] == idx2[3]
@def_same_line 3 :y  idx1[1] == idx2[1] && idx1[3] == idx2[3]
@def_same_line 3 :z  idx1[1] == idx2[1] && idx1[2] == idx2[2]
@def_same_line 3 :xy (idx2[2] - idx1[2]) == (idx2[1] - idx1[1]) && idx1[3] == idx2[3]
@def_same_line 3 :yx (idx2[2] - idx1[2]) == (idx1[1] - idx2[1]) && idx1[3] == idx2[3]
@def_same_line 3 :xz (idx2[3] - idx1[3]) == (idx2[1] - idx1[1]) && idx1[2] == idx2[2]
@def_same_line 3 :zx (idx2[3] - idx1[3]) == (idx1[1] - idx2[1]) && idx1[2] == idx2[2]
@def_same_line 3 :yz (idx2[3] - idx1[3]) == (idx2[2] - idx1[2]) && idx1[1] == idx2[1]
@def_same_line 3 :zy (idx2[3] - idx1[3]) == (idx1[2] - idx2[2]) && idx1[1] == idx2[1]
@def_same_line 3 :xyz (idx2[3] - idx1[3]) == (idx2[2] - idx1[2]) == (idx2[1] - idx1[1])
@def_same_line 3 :yxz (idx2[3] - idx1[3]) == (idx2[2] - idx1[2]) == (idx1[1] - idx2[1])
@def_same_line 3 :xzy (idx2[3] - idx1[3]) == (idx1[2] - idx2[2]) == (idx2[1] - idx1[1])
@def_same_line 3 :zyx (idx1[3] - idx2[3]) == (idx2[2] - idx1[2]) == (idx2[1] - idx1[1])


text:
######### CmdStan optimize program example  ###########

using CmdStan, Test, Statistics

ProjDir = dirname(@__FILE__)
cd(ProjDir) #do

  bernoulli = "
  data {
    int<lower=0> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
      y ~ bernoulli(theta);
  }
  "

  bernoullidata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])
  global stanmodel, rc, optim, cnames
  stanmodel = Stanmodel(Optimize(), name="bernoulli",
    output_format=:array, model=bernoulli);

  rc, optim, cnames = stan(stanmodel, bernoullidata, ProjDir, CmdStanDir=CMDSTAN_HOME);

  if rc == 0
    println()
    display(optim)
    println()
    println("Test round.(mean(optim[1][\"theta\"]), digits=1) ≈ 0.3")
    @test round.(mean(optim["theta"]), digits=1) ≈ 0.3
  end

  # Same with saved iterations
  stanmodel = Stanmodel(Optimize(save_iterations=true), name="bernoulli",
      nchains=1,output_format=:array, model=bernoulli);

  rc, optim, cnames = stan(stanmodel, bernoullidata, ProjDir,
    CmdStanDir=CMDSTAN_HOME);

  if rc == 0
    println()
    display(optim)
    println()
    println("""Test optim["theta"][end] ≈ 0.3 rtol=0.1""")
    @test optim["theta"][end] ≈ 0.3 rtol=0.1
  end

  #end # cd


text:
using Flux
using Flux.Tracker

f(x) = 3x^2 + 2x + 1
df(x) = Tracker.gradient(f, x; nest = true)[1]
df(3)

d2f(x) = Tracker.gradient(df, x; nest = true)[1]
d2f(3)

f(W, b, x) = W*x + b
Tracker.gradient(f, 2, 3, 4)

W = param(2) # 2.0 (tracked)
b = param(3) # 3.0 (tracked)

f(x) = W * x + b

grads = Tracker.gradient(() -> f(4), params(W, b))

(grads[W] , grads[b])

W = rand(2, 5)
b = rand(2)

predict(x) = W*x .+ b

function loss(x, y)
  ŷ = predict(x)
  sum((y .- ŷ).^2)
end

x, y = rand(5), rand(2) # Dummy data
loss(x, y) # ~ 3

W = param(W)
b = param(b)

gs = Tracker.gradient(() -> loss(x, y), params(W, b))

gs[W]

gs[b]

using Flux.Tracker: update!

D = gs[W]
update!(W, -0.1D)

loss(x, y)

using Flux, Flux.Data.MNIST, Statistics
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated

# Classify MNIST digits with a simple multi-layer-perceptron

imgs = MNIST.images()

# Stack images into one large batch
X = hcat(float.(reshape.(imgs, :))...)

labels = MNIST.labels()

# One-hot-encode the labels
Y = onehotbatch(labels, 0:9)

m = Chain(
  Dense(28^2, 32, relu),
  Dense(32, 10),
  softmax)

loss(x, y) = crossentropy(m(x), y)
accuracy(x, y) = mean(onecold(m(x)) .== onecold(y))

dataset = repeated((X, Y), 200)
evalcb = () -> @show(loss(X, Y))
opt = ADAM()

Flux.train!(loss, params(m), dataset, opt, cb = throttle(evalcb, 10))
accuracy(X, Y)

# Test set accuracy
tX = hcat(float.(reshape.(MNIST.images(:test), :))...) |> gpu
tY = onehotbatch(MNIST.labels(:test), 0:9) |> gpu

accuracy(tX, tY)

using Metalhead
using Metalhead: classify
using Images

# Switch to dataset images and classify the known (1000) images
cd(homedir()*"/PacktPub/Datasets/Images")
vgg = VGG19()

img01 = load("elephant.jpg")
classify(vgg, img01)

using Unicode

img02 = load("harry.jpg")
animal = uppercase(classify(vgg, img02))
println("Meow, I am a $animal perhaps ???'")

animal_types = Metalhead.labels(vgg)
animal_types[1:5]

# How sure are we?
#
probs = Metalhead.forward(vgg, img02)
sort(probs)[996:1000]

perm = sortperm(probs)
animal_types[perm[1000]]

animal_types[perm[999]]   # <= This is a breed of DOG. sorry Harry.

animal_types[287:288]    # If not a cat -- why not one of these ?

pwd(); include("KnetJL/mnist-mlp/mlp.jl")

MLP.main("--epochs 15")

using Knet, Images
cd(Knet.dir("examples"))
pwd()

include("mnist-mlp/mlp.jl")
MLP.main("--epochs 15")

xtrn,ytrn,xtst,ytst = MLP.mnist()
map(summary,(xtrn,ytrn,xtst,ytst))

using Images
mnistview(x,i) = Images.colorview(Gray,permutedims(x[:,:,1,i],(2,1)))

hcat([mnistview(xtst,i) for i=1:5]...)

ytst[1:5]'

include("fashion-mnist/fashion-mnist.jl")

FashionMNIST.main("--epochs 15")

include("../data/fashion-mnist.jl")

xtrn,ytrn,xtst,ytst = fmnist()
map(maximum,(xtrn,ytrn,xtst,ytst))

fmnistview(x,i) = Images.colorview(Gray,permutedims(x[:,:,1,i],(2,1)))
hcat([fmnistview(xtst,i) for i=1:10]...)

ytst[1:10]'

xtst


text:
module PSFs
using FourierTools: center_pos, FourierJoin
using FourierTools, NDTools, IndexFunArrays, SpecialFunctions, FFTW
using ZernikePolynomials
export PSFParams, sinc_r, jinc_r_2d, pupil_xyz, apsf, psf, k0, kxy, aplanatic_factor
export ModeWidefield, ModeConfocal, Mode4Pi

include("aplanatic.jl")
include("PSF_types.jl")
include("util.jl")
include("pupil_pol.jl")
include("pupil.jl")
include("aPSFs.jl")

"""
    psf(sz::NTuple, pp::PSFParams; sampling=get_sampling(sz, pp))

calculates the point spread function (psf), i.e. the image of a single (very small) emitter. Most of the parameters
(such as refractive index, numerical aperture, vacuum wavelength, aberrations etc.) are hidden in the parameter structure argument `pp`,
which should be generated via the `PSFParams()` constructor. See ``PSFParams()` for details.

See also:
+ apsf():  calculates the underlying aplitude point spread function (apsf)

Example:
```jdoctest
julia> pp = PSFParams(500.0,1.4,1.52);
julia> p = psf((128,128,128),pp; sampling=(50,50,100)); #;
```
"""
function psf(sz::NTuple, pp::PSFParams; sampling=nothing, use_resampling=true, return_amp=false) # unclear why the resampling seems to be so bad
    if use_resampling == false
        amp = apsf(sz, pp, sampling=sampling)
        if return_amp
            return amp_to_int(amp), amp
        else
            return amp_to_int(amp)
        end
    end
    extra_layers = 2
    small_sz=ceil.(Int,sz./2) .+ extra_layers
    big_sz = small_sz .* 2  # size after upsampling
    if ~isnothing(sampling)
        amp_sampling = sampling .* big_sz ./ small_sz
    else
        amp_sampling = nothing
    end
    amp = apsf(small_sz, pp, sampling=amp_sampling, center_kz=true)
    if pp.FFTPlan != nothing
        P1d = plan_fft(amp,(1,2), flags=pp.FFTPlan)
        P1id = plan_ifft(amp,(1,2), flags=pp.FFTPlan)
    end
    border_in = (0,0,ceil.(Int,sz[3]./2) ./ small_sz[3],0)
    mywin = collect(window_hanning((1,1,small_sz[3],1), border_in=border_in, border_out=1)) # for speed reasons the collect is faster
    amp = upsample2(amp .* mywin,dims=(1,2,3))
    res = sum(abs2.(amp),dims=4)[:,:,:,1]
    if true # any(isodd.(sz))
        if return_amp
            return select_region!(res,new_size=sz), amp
        else
            return select_region!(res,new_size=sz)
        end
    else
        if return_amp
            return res, amp
        else
            return res
        end
    end
end

end # module


text:
###
# using LinearAlgebra:sqrt
# using GeometryTypes
###
###
"""
Return vertex coordinates fixed to the unit sphere
"""
function vertex(x::T, y::T, z::T,scale) where T
    length = sqrt(x^2 + y^2 + z^2)
    return [(i * scale) / length for i in (x,y,z)]
end

function vertex(a::Array{T,1},scale) where T
    @assert length(a) == 3
    x,y,z = a[1],a[2],a[3]
    return vertex(x, y, z,scale)
end

"""
Find a middle point and project to the unit sphere
"""
function middle_point!(point_1,point_2,middle_point_cache,verts,scale)
    # We check if we have already cut this edge first
    # to avoid duplicated verts
    smaller_index = min(point_1, point_2)
    greater_index = max(point_1, point_2)
    key = "$smaller_index, $greater_index"
    if haskey(middle_point_cache,key)
        return middle_point_cache[key]
    end
    vert_1 = verts[point_1]
    vert_2 = verts[point_2]
    middle = [sum(i)/2 for i in zip(vert_1, vert_2)]
    append!(verts,[vertex(middle,scale)])
    index = length(verts)
    middle_point_cache[key] = index
    return index
end # function

function middle_point!(point_1::T,point_2::T,middle_point_cache,verts) where {T<:OffsetInteger}
    # to avoid duplicated verts
    smaller_index = min(point_1.i, point_2.i)
    greater_index = max(point_1.i, point_2.i)
    key = "$smaller_index, $greater_index"
    if haskey(middle_point_cache,key)
        return middle_point_cache[key]
    end
    vert_1 = verts[point_1]
    vert_2 = verts[point_2]
    middle = [sum(i)/2 for i in zip(vert_1, vert_2)]
    append!(verts,[vertex(middle,scale)])
    index = length(verts)
    middle_point_cache[key] = index
    return index
end # function

"""
Create icosphere
"""
function icosphere(subdiv=3,scale = 1.0)
    PHI = (1. + sqrt(5.)) / 2.
    verts = [(vertex(-1., PHI, 0.,scale)),
        (vertex(1., PHI, 0.,scale)),
        (vertex(-1., -PHI, 0.,scale)),
        (vertex(1., -PHI, 0.,scale)),
        (vertex(0., -1., PHI,scale)),
        (vertex(0., 1., PHI,scale)),
        (vertex(0., -1., -PHI,scale)),
        (vertex(0., 1., -PHI,scale)),
        (vertex(PHI, 0., -1.,scale)),
        (vertex(PHI, 0., 1.,scale)),
        (vertex(-PHI, 0., -1.,scale)),
        (vertex(-PHI, 0., 1.,scale))]
    faces =[[1, 12, 6], # 5 faces around point 0
            [1, 6, 2],
            [1, 2, 8],
            [1, 8, 11],
            [1, 11, 12],
            [2, 6, 10], # Adjacent faces
            [6, 12, 5],
            [12, 11, 3],
            [11, 8, 7],
            [8, 2, 9],
            [4, 10, 5],  # 5 faces around 3
            [4, 5, 3],
            [4, 3, 7],
            [4, 7, 9],
            [4, 9, 10],
            [5, 10, 6],# Adjacent faces
            [3, 5, 12],
            [7, 3, 11],
            [9, 7, 8],
            [10, 9, 2]]

    middle_point_cache = Dict{String,Int}()
    for i in 1:subdiv
        faces_subdiv = Array{Int64,1}[]
        for tri in faces
            v1 = middle_point!(tri[1], tri[2],middle_point_cache,verts,scale)
            v2 = middle_point!(tri[2], tri[3],middle_point_cache,verts,scale)
            v3 = middle_point!(tri[3], tri[1],middle_point_cache,verts,scale)
            append!(faces_subdiv,[[tri[1], v1, v3]])
            append!(faces_subdiv,[[tri[2], v2, v1]])
            append!(faces_subdiv,[[tri[3], v3, v2]])
            append!(faces_subdiv,[[v1, v2, v3]])
            faces = faces_subdiv
        end
    end
    v = array_to_vertex(verts)
    f = array_to_faces(faces)
    f= convert(Vector{TriangleFace{OffsetInteger{-1,UInt32}}},f)
    sphere = Mesh(v,f)
    return sphere
end


text:
using Documenter

makedocs(
    sitename = "BAWJ",
    format = Documenter.HTML(),
    pages = [
        "Business Applications with Julia" => "index.md",
        "Blog" => "blog.md",
        "1 - The Application" => "chapter1.md",
        "Section 1 - Experimenting with Modules" => [
            "2 - Setup the Development Environment" => "chapter2.md",
            "3 - Create and Test Domain.jl" => "chapter3.md",
            "4 - Create and Test API.jl" => "chapter4.md",
            "5 - Create and Test Infrastructure.jl" => "chapter5.md",
            "6 - The Accounts Package from a User Point of View" => "chapter6.md",
            "7 - Documenting your module" => "chapter7.md",
        ],
        "Section 2 -  The Accounts Receivable Module" => [
            "8 - The Design" => "chapter8.md",
            "9 - The Sub-Module Domain" => "chapter9.md",
            "10 - The Sub-Module API" => "chapter10.md",
            "11 - The Sub-Module Infrastructure" => "chapter11.md",
            "12 - The Main-Module AppliAR.jl" => "chapter12.md",
        ],
        "Section 3 - Using containers" => [
            "13 - Creating SSH Enabled Containers" => "chapter13.md",
            "14 - Run the Application from a Notebook" => "chapter14.md",
            "15 - Running the Website from a Container" => "chapter15.md",
        ],
        "Section 4 - Miscellaneous" => [
            "16 - Logging" => "chapter16.md",
            "17 - Implementing AppliAR (Accounts Receivable) Package" => "chapter17.md",
            "18 - PosgreSQL" => "chapter18.md",
            "19 - Jenkins-slave" => "chapter19.md"
        ],
        "Installation Instruction" => "appendix.md",
        "IJulia Notebook" => "notebook.md",
        "Overview Julia" => "overview.md",
        "References" => "references.md"
    ]
)

max_stars_repo_name: rbontekoe/BAWJ

max_issues_repo_name: rbontekoe/BAWJ

max_forks_repo_name: rbontekoe/BAWJ


text:
# TODO: move to ACEatoms
# # -> rand_config?
# function rand_nhd(Nat, J::ScalarBasis, species = :X)
#    zlist = ZList(species)
#    Rs = [ rand_vec(J) for _ = 1:Nat ]
#    Zs = [ rand(zlist.list) for _ = 1:Nat ]
#    z0 = rand(zlist.list)
#    return Rs, Zs, z0
# end

# rand_config(species; kwargs...) =
#       rand_config(ZList(species); kwargs...)

# rand_config(V::AbstractCalculator; kwargs...) =
#       rand_config(zlist(V); kwargs...)

# function rand_config(zlist::Union{ZList, SZList};
#                      absrattle = 0.0, relrattle = 0.2, repeat = 3,
#                      kwargs...)
#    # start with the longest rnn
#    rnns = rnn.(zlist.list)
#    sym = chemical_symbol( zlist.list[findmax(rnns)[2]] )
#    at = bulk(sym; kwargs...) * repeat
#    for n = 1:length(at)
#       at.Z[n] = rand(zlist.list)
#    end
#    rattle!(at, maximum(rnns) * relrattle + absrattle)
#    return at
# end

# @info("Check several properties of PIBasis")
# for species in (:X, :Si, [:C, :O, :H]), N = 1:5
#    local AA, AAdag, dAA, dAAdag, dagbasis, basis, Rs, Zs, z0
#    maxdeg = 7
#    Nat = 15
#    P1 = ACE.RnYlm1pBasis(Pr; species = species)
#    dagbasis = ACE.PIBasis(P1, N, D, maxdeg)
#    basis = standardevaluator(dagbasis)
#    @info("species = $species; N = $N; length = $(length(basis))")
#    @info("test (de-)serialisation")
#    # only require dagbasis to deserialize correctly
#    println(@test all(JuLIP.Testing.test_fio(dagbasis)))
#    @info("Check Permutation invariance")
#    for ntest = 1:20
#       Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
#       p = randperm(length(Rs))
#       print_tf(@test(evaluate(basis, Rs, Zs, z0) ≈
#                      evaluate(basis, Rs[p], Zs[p], z0)))
#    end
#    println()
#    @info("Check gradients")
#    for ntest = 1:20
#       Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
#       AA = evaluate(basis, Rs, Zs, z0)
#       dAA = evaluate_d(basis, Rs, Zs, z0)
#       Us = [ rand(eltype(Rs)) .- 0.5 for _=1:length(Rs) ]
#       dAA_dUs = transpose.(dAA) * Us
#       errs = []
#       for p = 2:12
#          h = 0.1^p
#          AA_h = evaluate(basis, Rs + h * Us, Zs, z0)
#          dAA_h = (AA_h - AA) / h
#          # @show norm(dAA_h - dAA_dUs, Inf)
#          push!(errs, norm(dAA_h - dAA_dUs, Inf))
#       end
#       success = (/(extrema(errs)...) < 1e-3) || (minimum(errs) < 1e-10)
#       print_tf(@test success)
#    end
#
#    println()
#    @info("Check Standard=DAG Evaluator")
#    for ntest = 1:20
#       Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
#       AA = evaluate(basis, Rs, Zs, z0)
#       AAdag = evaluate(dagbasis, Rs, Zs, z0)
#       print_tf(@test AA ≈ AAdag)
#
#       dAA = evaluate_d(basis, Rs, Zs, z0)
#       dAAdag = evaluate_d(dagbasis, Rs, Zs, z0)
#       print_tf(@test dAA ≈ dAAdag)
#    end
#    println()
# end
# println()

#---

# #---
# @info("Basis construction and evaluation checks")
# @info("check single species")
# Nat = 15
# Rs, Zs, z0 = rand_nhd(Nat, Pr, :X)
# B = evaluate(rpibasis, Rs, Zs, z0)
# println(@test(length(rpibasis) == length(B)))
# dB = evaluate_d(rpibasis, Rs, Zs, z0)
# println(@test(size(dB) == (length(rpibasis), length(Rs))))
# B_, dB_ = evaluate_ed(rpibasis, Rs, Zs, z0)
# println(@test (B_ ≈ B) && (dB_ ≈ dB))
#
# #---
# @info("check multi-species")
# maxdeg = 5
# Pr = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)
# species = [:C, :O, :H]
# P1 = ACE.RnYlm1pBasis(Pr; species = species, D = D)
# basis = ACE.RPIBasis(P1, N, D, maxdeg)
# Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
# B = evaluate(basis, Rs, Zs, z0)
# println(@test(length(basis) == length(B)))
# dB = evaluate_d(basis, Rs, Zs, z0)
# println(@test(size(dB) == (length(basis), length(Rs))))
# B_, dB_ = evaluate_ed(basis, Rs, Zs, z0)
# println(@test (B_ ≈ B) && (dB_ ≈ dB))
#
# #---
#
# degrees = [ 12, 10, 8, 8, 8, 8 ]
#
# @info("Check a few basis properties ")
# # for species in (:X, :Si) # , [:C, :O, :H])
# for species in (:X, :Si, [:C, :O, :H]), N = 1:length(degrees)
#    local Rs, Zs, z0, B, dB, basis, D, P1, Nat
#    Nat = 15
#    D = SparsePSHDegree()
#    P1 = ACE.RnYlm1pBasis(Pr; species = species)
#    basis = ACE.RPIBasis(P1, N, D, degrees[N])
#    @info("species = $species; N = $N; deg = $(degrees[N]); len = $(length(basis))")
#    @info("   check (de-)serialization")
#    println(@test(all(JuLIP.Testing.test_fio(basis))))
#    @info("   isometry and permutation invariance")
#    for ntest = 1:30
#       Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
#       Rsp, Zsp = ACE.rand_sym(Rs, Zs)
#       print_tf(@test(evaluate(basis, Rs, Zs, z0) ≈
#                      evaluate(basis, Rsp, Zsp, z0)))
#    end
#    println()
#    @info("   check derivatives")
#    for ntest = 1:30
#       Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
#       B = evaluate(basis, Rs, Zs, z0)
#       dB = evaluate_d(basis, Rs, Zs, z0)
#       Us = [ rand(eltype(Rs)) .- 0.5 for _=1:length(Rs) ]
#       dB_dUs = transpose.(dB) * Us
#       errs = []
#       for p = 2:12
#          h = 0.1^p
#          B_h = evaluate(basis, Rs + h * Us, Zs, z0)
#          dB_h = (B_h - B) / h
#          # @show norm(dAA_h - dAA_dUs, Inf)
#          push!(errs, norm(dB_h - dB_dUs, Inf))
#       end
#       success = (/(extrema(errs)...) < 1e-3) || (minimum(errs) < 1e-10)
#       print_tf(@test success)
#    end
#    println()
#    @info("   check combine")
#    coeffs = randcoeffs(basis)
#    V = combine(basis, coeffs)
#    Vst = standardevaluator(V)
#    for ntest = 1:30
#       Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
#       v = evaluate(V, Rs, Zs, z0)
#       vst = evaluate(Vst, Rs, Zs, z0)
#       cdotB = dot(coeffs, evaluate(basis, Rs, Zs, z0))
#       print_tf(@test v ≈ cdotB ≈ vst)
#    end
#    println()
#    @info("   check graph evaluator")
#    basisst = standardevaluator(basis)
#    for ntest = 1:30
#       env = ACE.rand_nhd(Nat, Pr, species)
#       print_tf(@test evaluate(basisst, env...) ≈ evaluate(basis, env...))
#       print_tf(@test evaluate_d(basisst, env...) ≈ evaluate_d(basis, env...))
#    end
#    println()
# end
#


text:
__precompile__()
module S3

using AWSCore
using AWSS3
using Dates
using FilePathsBase
using Memento

import FilePathsBase: sync

export S3Path, sync, upload

const logger = getlogger(@__MODULE__)

# Register the module level logger at runtime so that folks can access the logger via
# `getlogger(MyModule)`. If this line is not included then the precompiled
# `MyModule.logger` won't be registered at runtime.
function __init__()
    Memento.register(logger)
end

# Couple extra methods that should probably be included in AWSS3 at some point.
"""
    upload(src::AbstractPath, dest::S3Path)

Uploads a local file to the s3 path specified by `dest`.
"""
upload

const upload = Base.cp

# Note: naming copied from Go SDK:
# https://docs.aws.amazon.com/sdk-for-go/v1/developer-guide/s3-example-presigned-urls.html
function presign(
    path::S3Path,
    duration::Period=Hour(1);
    config::AWSConfig=aws_config(),
)
    AWSS3.s3_sign_url(config, path.bucket, path.key, Dates.value(Second(duration)))
end

# TODO: Remove the sync methods below as they are now pirating.
@deprecate sync(src::AbstractString, dest::AbstractString; delete::Bool=false, config::AWSConfig=aws_config()) sync(Path(src), Path(dest); delete=delete)

end  # S3


text:
# test 5 lines with scenario specific tau
using Distributed;
addprocs(30);
@everywhere include("loadMod.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();

caseList = [13,33,123];
T = 24;
Δt = 0.25;
NN = 1000;
N = 5;

dataList = Dict();
for ci in 1:length(caseList)
    for j in procs()
        remotecall_fetch(readInData_tau,j,ci,caseList,T);
        #remotecall_fetch(readInData_old,j,T,ωSet0,10000,0);
    end

    # train the stochastic strategy
    # cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(τ, T, Δt, N, true, false, 20, 20, cutDictPG);
    cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(T, Δt, N, false, false,
        max(Int64(round(500/N)),20), max(Int64(round(500/N)),20),Dict());

    pathListData = pmap(i -> simuPath(T,pDistr), 1:NN);
    pathDict = Dict();
    for i in 1:NN
        pathDict[i] = pathListData[i];
    end
    solSDDP, LBSDDP, costSDDP = exeForward(T, Δt, NN, false, pathDict);
    listSDDP = [costSDDP[i] for i in 1:NN];
    meanSDDP = mean(listSDDP);
    sigmaSDDP = std(listSDDP);

    for j in procs()
        remotecall_fetch(readInData_tau,j,ci,caseList,T,1e4,0,"Same");
    end

    cutDicts,LBHists,UBHists,UBuHists,UBlHists,timeHists = solveMain(T, Δt, N, false, false,
        max(Int64(round(500/N)),20), max(Int64(round(500/N)),20),Dict());

    solSDDPs, LBSDDPs, costSDDPs = exeForward(T, Δt, NN, false, pathDict);
    listSDDPs = [costSDDPs[i] for i in 1:NN];
    meanSDDPs = mean(listSDDPs);
    sigmaSDDPs = std(listSDDPs);

    dataList[ci] = [LBHist,LBSDDP,listSDDP,meanSDDP,sigmaSDDP,
                    LBHists,LBSDDPs,listSDDPs,meanSDDPs,sigmaSDDPs];

    save("tauResults_specific.jld","data",dataList);
end


text:
immutable Aesthetic
    x::Symbol
    y::Symbol
end

function Aesthetic(kwargs::Dict)
    return Aesthetic(kwargs[:x], kwargs[:y])
end

function aes(;kwargs...)
    return Aesthetic(Dict(kwargs))
end


text:
# adiabatic processing test
n=4; T = 3
# ins,ans = generate(n)
# pH = pHamiltonian(ins,n)
# aqc = AQC(n)
# aqc|>AdiabaticProcessing(pH,QuExpmV();dt=1e-2,maxtime=T)
# @test run!(aqc) == T
# @test round(norm(coeff(aqc.state)[ans[1]+1])^2,3) == 0.130

# cooling-assisted processing test
ins = generate(n)
pH = pHamiltonian(ins,n)
eigen = eigfact(full(pH))
val,idx = findmin(eigen[:values])
maxEigen = maximum(eigen[:values])

aqc = AQC(n)
cooled_aqc = AQC(n)

ap1 = AdiabaticProcessing(pH,QuExpmV();dt=1e-2,maxtime=T,tlist=(0.0,1/3))
ap2 = AdiabaticProcessing(pH,QuExpmV();dt=1e-2,maxtime=T,tlist=(1/3,2/3))
ap3 = AdiabaticProcessing(pH,QuExpmV();dt=1e-2,maxtime=T,tlist=(2/3,1.0))
cooled_aqc|>ap1|>DemonCooling(hamiltonian(ap1,ap1.tlist[2]); k=10, tol=3, optimize=:off)
cooled_aqc|>ap2|>DemonCooling(hamiltonian(ap2,ap2.tlist[2]); k=10, tol=3, optimize=:off)
cooled_aqc|>ap3

aqc|>ap1|>ap2|>ap3
@show run!(cooled_aqc)
@show abs(dot(coeff(cooled_aqc.state),eigen[:vectors][:,idx]))
# @show run!(aqc)
# data = @parallel (+) for i=1:1
#   cooled_aqc.state = BellState(n)
#   run_time = run!(cooled_aqc)
#   suc_prob = abs(dot(coeff(cooled_aqc.state),eigen[:vectors][:,idx]))
#   Float64[run_time,suc_prob]
# end

# @show data

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


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

# xgboost_wrapper.h

"Calls an xgboost API function and correctly reports errors."
macro xgboost_ccall(f, argTypes, args...)
    argTypes = eval(argTypes)
    return quote
        err = ccall(($f, _xgboost), Int64, ($(argTypes...),), $(args...))
        if err != 0
            errMsg = unsafe_string(ccall((:XGBGetLastError, _xgboost), Ptr{UInt8}, ()))
            error("Call to XGBoost C function "*string($f)*" failed: $errMsg")
        end
    end
end

function XGDMatrixCreateFromFile(fname::Compat.ASCIIString, slient::Int32)
    handle = Ref{Ptr{Void}}()
    @xgboost_ccall(
        :XGDMatrixCreateFromFile,
        (Ptr{UInt8}, Int32, Ref{Ptr{Void}}),
        fname, slient, handle
    )
    return handle[]
end

function XGDMatrixCreateFromCSC(data::SparseMatrixCSC)
    handle = Ref{Ptr{Void}}()
    @xgboost_ccall(
        :XGDMatrixCreateFromCSC,
        (Ptr{UInt64}, Ptr{UInt32}, Ptr{Float32}, UInt64, UInt64, Ref{Ptr{Void}}),
        convert(Array{UInt64, 1}, data.colptr - 1),
        convert(Array{UInt32, 1}, data.rowval - 1), convert(Array{Float32, 1}, data.nzval),
        convert(UInt64, size(data.colptr)[1]),
        convert(UInt64, nnz(data)),
        handle
    )
    return handle[]
end

function XGDMatrixCreateFromCSCT(data::SparseMatrixCSC)
    handle = Ref{Ptr{Void}}()
    @xgboost_ccall(
        :XGDMatrixCreateFromCSR,
        (Ptr{UInt64}, Ptr{UInt32}, Ptr{Float32}, UInt64, UInt64, Ref{Ptr{Void}}),
        convert(Array{UInt64, 1}, data.colptr - 1),
        convert(Array{UInt32, 1}, data.rowval - 1), convert(Array{Float32, 1}, data.nzval),
        convert(UInt64, size(data.colptr)[1]),
        convert(UInt64, nnz(data)),
        handle
    )
    return handle[]
end

function XGDMatrixCreateFromMat(data::Array{Float32, 2}, missing::Float32)
    XGDMatrixCreateFromMatT(transpose(data), missing)
end

function XGDMatrixCreateFromMatT(data::Array{Float32, 2}, missing::Float32)
    nrow = size(data)[2]
    ncol = size(data)[1]
    handle = Ref{Ptr{Void}}()
    @xgboost_ccall(
        :XGDMatrixCreateFromMat,
        (Ptr{Float32}, UInt64, UInt64, Float32, Ref{Ptr{Void}}),
        data, nrow, ncol, missing, handle
    )
    return handle[]
end

function XGDMatrixSliceDMatrix(handle::Ptr{Void}, idxset::Array{Int32, 1}, len::UInt64)
    ret = Ref{Ptr{Void}}()
    @xgboost_ccall(
        :XGDMatrixSliceDMatrix,
        (Ptr{Void}, Ptr{Int32}, UInt64, Ref{Ptr{Void}}),
        handle, idxset, len, ret
    )
    return ret[]
end

function XGDMatrixFree(handle::Ptr{Void})
    @xgboost_ccall(
        :XGDMatrixFree,
        (Ptr{Void},),
        handle
    )
end

function XGDMatrixSaveBinary(handle::Ptr{Void}, fname::Compat.ASCIIString, slient::Int32)
    @xgboost_ccall(
        :XGDMatrixSaveBinary,
        (Ptr{Void}, Ptr{UInt8}, Int32),
        handle, fname, slient
    )
end

function XGDMatrixSetFloatInfo(handle::Ptr{Void}, field::Compat.ASCIIString,
                               array::Array{Float32, 1}, len::UInt64)
    @xgboost_ccall(
        :XGDMatrixSetFloatInfo,
        (Ptr{Void}, Ptr{UInt8}, Ptr{Float32}, UInt64),
        handle, field, array, len
    )
end

function XGDMatrixSetUIntInfo(handle::Ptr{Void}, field::Compat.ASCIIString,
                              array::Array{UInt32, 1}, len::UInt64)
    @xgboost_ccall(
        :XGDMatrixSetUIntInfo,
        (Ptr{Void}, Ptr{UInt8}, Ptr{UInt32}, UInt64),
        handle, field, array, len
    )
end

function XGDMatrixSetGroup(handle::Ptr{Void}, array::Array{UInt32, 1}, len::UInt64)
    @xgboost_ccall(
        :XGDMatrixSetGroup,
        (Ptr{Void}, Ptr{UInt32}, UInt64),
         handle, array, len
    )
end

function XGDMatrixGetFloatInfo(handle::Ptr{Void}, field::Compat.ASCIIString, outlen::Array{UInt64, 1})
    ret = Ref{Ptr{Float32}}()
    @xgboost_ccall(
        :XGDMatrixGetFloatInfo,
        (Ptr{Void}, Ptr{UInt8}, Ptr{UInt64}, Ref{Ptr{Float32}}),
         handle, field, outlen, ret
    )
    return ret[]
end

function XGDMatrixGetUIntInfo(handle::Ptr{Void}, field::Compat.ASCIIString, outlen::Array{UInt64, 1})
    ret = Ref{Ptr{UInt32}}()
    @xgboost_ccall(
        :XGDMatrixGetUIntInfo,
        (Ptr{Void}, Ptr{UInt8}, Ptr{UInt64}, Ref{Ptr{UInt32}}),
         handle, field, outlen, ret
    )
    return ret[]
end

function XGDMatrixNumRow(handle::Ptr{Void})
    ret = Ref{UInt64}()
    @xgboost_ccall(
        :XGDMatrixNumRow,
        (Ptr{Void}, Ref{UInt64}),
         handle, ret
    )
    return ret[]
end

function JLGetFloatInfo(handle::Ptr{Void}, field::Compat.ASCIIString)
    len = UInt64[1]
    ptr = XGDMatrixGetFloatInfo(handle, field, len)
    return unsafe_wrap(Array, ptr, len[1])
end

function JLGetUintInfo(handle::Ptr{Void}, field::Compat.ASCIIString)
    len = UInt64[1]
    ptr = XGDMatrixGetUIntInfo(handle, field, len)
    return unsafe_wrap(Array, ptr, len[1])
end

function XGBoosterCreate(cachelist::Array{Ptr{Void}, 1}, len::Int64)
    handle = Ref{Ptr{Void}}()
    @xgboost_ccall(
        :XGBoosterCreate,
        (Ptr{Ptr{Void}}, UInt64, Ref{Ptr{Void}}),
        cachelist, len, handle
    )
    return handle[]
end

function XGBoosterFree(handle::Ptr{Void})
    @xgboost_ccall(
        :XGBoosterFree,
        (Ptr{Void}, ),
        handle
    )
end

function XGBoosterSetParam(handle::Ptr{Void}, key::Compat.ASCIIString, value::Compat.ASCIIString)
    @xgboost_ccall(
        :XGBoosterSetParam,
        (Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}),
        handle, key, value
    )
end

function XGBoosterUpdateOneIter(handle::Ptr{Void}, iter::Int32, dtrain::Ptr{Void})
    @xgboost_ccall(
        :XGBoosterUpdateOneIter,
        (Ptr{Void}, Int32, Ptr{Void}),
        handle, iter, dtrain
    )
end

function XGBoosterBoostOneIter(handle::Ptr{Void}, dtrain::Ptr{Void},
                               grad::Array{Float32, 1},
                               hess::Array{Float32, 1},
                               len::UInt64)
    @xgboost_ccall(
        :XGBoosterBoostOneIter,
        (Ptr{Void}, Ptr{Void}, Ptr{Float32}, Ptr{Float32}, UInt64),
        handle, dtrain, grad, hess, len
    )
end

function XGBoosterEvalOneIter(handle::Ptr{Void}, iter::Int32,
                              dmats::Array{Ptr{Void}, 1},
                              evnames::Array{Compat.ASCIIString, 1}, len::UInt64)
    msg = Ref{Ptr{UInt8}}()
    @xgboost_ccall(
        :XGBoosterEvalOneIter,
        (Ptr{Void}, Int32, Ptr{Ptr{Void}}, Ptr{Ptr{UInt8}}, UInt64, Ref{Ptr{UInt8}}),
        handle, iter, dmats, evnames, len, msg
    )
    return unsafe_string(msg[])
end

function XGBoosterPredict(handle::Ptr{Void}, dmat::Ptr{Void}, output_margin::Int32,
                          ntree_limit::UInt32, len::Array{UInt64, 1})
    ret = Ref{Ptr{Float32}}()
    @xgboost_ccall(
        :XGBoosterPredict,
        (Ptr{Void}, Ptr{Void}, Int32, UInt32, Ptr{UInt64}, Ref{Ptr{Float32}}),
        handle, dmat, output_margin, ntree_limit, len, ret
    )
    return ret[]
end

function XGBoosterLoadModel(handle::Ptr{Void}, fname::Compat.ASCIIString)
    @xgboost_ccall(
        :XGBoosterLoadModel,
        (Ptr{Void}, Ptr{UInt8}),
        handle, fname
    )
end

function XGBoosterSaveModel(handle::Ptr{Void}, fname::Compat.ASCIIString)
    @xgboost_ccall(
        :XGBoosterSaveModel,
        (Ptr{Void}, Ptr{UInt8}),
        handle, fname
    )
end

function XGBoosterDumpModel(handle::Ptr{Void}, fmap::Compat.ASCIIString, with_stats::Int64)
    data = Ref{Ptr{Ptr{UInt8}}}()
    out_len = Ref{UInt64}(0)
    @xgboost_ccall(
        :XGBoosterDumpModel,
        (Ptr{Void}, Ptr{UInt8}, Int64, Ref{UInt64}, Ref{Ptr{Ptr{UInt8}}}),
        handle, fmap, with_stats, out_len, data
    )
    return unsafe_wrap(Array, data[], out_len[])
end


text:
using Utils

using Chain

import IterTools

const file = find_input(@__FILE__)
const input = file |> slurp

struct Packet
    version::UInt8
    type::UInt8
    content::Union{Int,Vector{Packet}}
end

function bitvector(input::AbstractString)
    @assert length(input) % 2 == 0 "input has length $(length(input))"
    bytes = @chain input IterTools.partition(_, 2) map(c -> parse(UInt8, join(c), base=16), _)
    re = BitVector((byte & (1 << i)) != 0 for byte=bytes for i=7:-1:0)
    # Make sure there's enought to read packets of 4.
    push!(re, false, false, false)
end

decimal(bits::AbstractArray{Bool}) = foldl((l,r) -> (l << 1) | r, bits, init=0)

function decode(bv::BitVector, limit, pos=1)::Tuple{Int,Vector{Packet}}
    found = Packet[]
    while pos <= length(bv) && limit != 0
        version = decimal(bv[pos:pos+2])
        pos += 3
        type = decimal(bv[pos:pos+2])
        pos += 3
        if type == 4
            num = 0
            while true
                chunk = decimal(bv[pos:pos+4])
                pos += 5
                num = (num << 4) | (chunk & 0b1111)
                ((chunk & 0b10000) == 0) && break
            end
            push!(found, Packet(version, type, num))
        else
            if bv[pos]
                len = decimal(bv[pos+1:pos+11])
                pos += 1 + 11
                inc, children = decode(bv[pos:end], len)
            else
                len = decimal(bv[pos+1:pos+15])
                pos += 1 + 15
                inc, children = decode(bv[pos:pos+len-1], -1)
            end
            pos += inc - 1
            push!(found, Packet(version, type, children))
        end
        limit > 0 && (limit -= 1)
    end
    pos, found
end

decode(string::AbstractString) = @chain string bitvector decode(1)

function sumversions(paket::Packet)::Int
    re = paket.version
    if !isa(paket.content, Number)
        re += map(sumversions, paket.content) |> sum
    end
    re
end

function part1(input)::Int
    @chain input begin
        decode
        _[2]
        only
        sumversions
        sum
    end
end

assertequal(part1(input), 965)

function part2(pak::Packet)::Int
    pak.type == 4 && return pak.content
    evald = part2.(pak.content)
    if pak.type == 0
        sum(evald)
    elseif pak.type == 1
        prod(evald)
    elseif pak.type == 2
        minimum(evald)
    elseif pak.type == 3
        maximum(evald)
    else
        @assert length(evald) == 2
        left, right = evald
        bool = if pak.type == 5
            left > right
        elseif pak.type == 6
            left < right
        elseif pak.type == 7
            left == right
        end
        convert(Int, bool)
    end
end

part2(input) = decode(input)[2] |> only |> part2

assertequal(part2(input), 116672213160)

using Test

@testset "utilities" begin
    @test bitvector("D2")[1:8] == [true, true, false, true, false, false, true, false]
    @test bitvector("FE")[1:8] == [true, true, true, true, true, true, true, false]
    @test bitvector("FED2")[1:16] == [true, true, true, true, true, true, true, false, true, true, false, true, false, false, true, false]
    @test decimal([true,true,false]) == 6
    @test part1("8A004A801A8002F478") == 16
    @test part1("620080001611562C8802118E34") == 12
    @test part1("C0015000016115A2E0802F182340") == 23
    @test part1("A0016C880162017C3686B18A3D4780") == 31
end

max_stars_repo_name: CmdQ/AoC

max_issues_repo_name: CmdQ/AoC

max_forks_repo_name: CmdQ/AoC


text:
using BinDeps
using Conda

@BinDeps.setup

libeccodes = library_dependency("libeccodes")
Conda.add_channel("conda-forge")
Conda.update()
provides(Conda.Manager, Dict("eccodes"=>libeccodes))

@BinDeps.install Dict(:libeccodes => :libeccodes)


text:
module BouncingEllipsoidSimulation

using ModiaLang
import Modia3D
using  Modia3D.ModiaInterface

BouncingEllipsoid = Model(
    boxHeigth = 0.1,
    gravField = UniformGravityField(g=9.81, n=[0, -1, 0]),
    world = Object3D(feature=Scene(gravityField=:gravField,
                                   visualizeFrames=false,
                                   defaultFrameLength=0.2,
                                   visualizeBoundingBox = true,
                                   enableContactDetection=true,
                                   visualizeContactPoints=false)),
    worldFrame = Object3D(parent=:world, feature=Visual(shape=CoordinateSystem(length=0.5))),
    ground = Object3D(parent=:world,
                      translation=:[0.0,-boxHeigth/2,0.0],
                      feature=Solid(shape=Box(lengthX=4.0, lengthY=:boxHeigth, lengthZ=3.0),
                                    visualMaterial=VisualMaterial(color="DarkGreen", transparency=0.5),
                                    solidMaterial="Steel",
                                    collision=true)),
    ellipsoid = Object3D(feature=Solid(shape=Ellipsoid(lengthX=0.1, lengthY=0.2, lengthZ=0.3),
                                       visualMaterial=VisualMaterial(color="Blue"),
                                       solidMaterial="Steel",
                                       collision=true)),
    free = FreeMotion(obj1=:world, obj2=:ellipsoid, r=Var(init=[0.0, 1.0, 0.0]), w=Var(init=[5.0, 0.0, -2.0]))
)

bouncingEllipsoid = @instantiateModel(buildModia3D(BouncingEllipsoid), unitless=true, log=false, logStateSelection=false, logCode=false)

#@show bouncingEllipsoid.parameterExpressions
#@show bouncingEllipsoid.parameters

stopTime = 2.5
tolerance = 1e-8
requiredFinalStates = [-0.8720220901260356, -2.706225861204731, 2.1494619380602984, -0.6654182338218089, -7.480644002905684, 0.8158622722948691, 19.6901637595772, -0.5127178856804271, 0.5844275371325878, 7.417393481447869, -0.15143390608966664, -0.19610822193495095]
simulate!(bouncingEllipsoid, stopTime=stopTime, tolerance=tolerance, log=true, logStates=true, logEvents=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(bouncingEllipsoid, ["free.r" "free.rot"; "free.v" "free.w"], figure=1)

end


text:
module TableOfElements
using JSON
using AtomicLevels
using Unicode
using Printf

struct Element
    name::String
    symbol::String
    Z::Int
    mass::Float64 # amu
    conf::Configuration
    Ip::Float64
    Ip_eV::Float64
    χ::Float64 # Pauling units
end

struct ElementsTable
    elements::Vector{Element}
    Z::Dict{String,Int}
end

function _load_table_of_elements()
    mp = r"([0-9\.]+)\([0-9]+\)"
    Z = Dict{String,Int}()
    elements = map(JSON.parsefile(joinpath(dirname(@__FILE__), "../data/data.json"))) do element
        mstr = element["atomicMass"]
        mass = if typeof(mstr) <: String
            m = match(mp, mstr)
            if m != nothing
                parse(Float64, m[1])
            else
                Inf
            end
        elseif typeof(mstr) <: Vector
            mstr[1]
        else
            Inf
        end

        getvalue(name, default=NaN) = (name in keys(element) && element[name] ≠ "" ?
                                       element[name] : default)

        Ip = getvalue("ionizationEnergy")
        χ = getvalue("electronegativity")

        conf = parse(Configuration{Orbital},
                     replace(element["electronicConfiguration"], "." => " "))
        symbol = element["symbol"]
        ZZ = element["atomicNumber"]
        Z[symbol] = ZZ
        Element(element["name"], symbol, ZZ,
                mass, conf, Ip, 0.01036410Ip, χ)
    end
    ElementsTable(elements, Z)
end

const data = _load_table_of_elements()

get_element(Z::Integer) = data.elements[Z]
get_element(symbol::String) = get_element(data.Z[titlecase(symbol)])
get_element(symbol::Symbol) = get_element("$(symbol)")

import Base: show
function show(stream::IO, element::Element)
    write(stream,
          @sprintf("%s [%s], Z = %i, mass = %.2f",
                   element.name, element.symbol,
                   element.Z, element.mass))
    !isnan(element.Ip) &&
        write(stream,
              @sprintf(", Ip = %.2f kJ/mol = %.2f eV",
                       element.Ip, element.Ip_eV))
    !isnan(element.χ) &&
        write(stream, @sprintf(", χ = %.2f", element.χ))
end

export get_element, show

end # module


text:
# using FunHPC
using Base.Test

unshift!(LOAD_PATH, "../src")

include("FunsTest.jl")
include("FoldableTest.jl")
include("FunctorTest.jl")
include("StencilFunctorTest.jl")
include("MaybesTest.jl")
include("MonadTest.jl")
include("MultiDictsTest.jl")

include("CommTest.jl")
include("GIDsTest.jl")
include("FunRefsTest.jl")
include("ParTest.jl")

using Comm
run_main() do
    FunsTest.main()
    FoldableTest.main()
    FunctorTest.main()
    StencilFunctorTest.main()
    MaybesTest.main()
    MonadTest.main()
    MultiDictsTest.main()

    CommTest.main()
    GIDsTest.main()
    FunRefsTest.main()
    ParTest.main()
end


text:
# Deprecations in FITSIO module. Additional deprecations in FITSIO.Libcfitsio
# are in cfitsio.jl.

import Base: @deprecate

# Deprecated in v0.6
@deprecate readkey(hdu::HDU, key) read_key(hdu, key)
@deprecate readheader(hdu::HDU) read_header(hdu::HDU)
@deprecate getcomment(hdr::FITSHeader, key) get_comment(hdr, key)
@deprecate setcomment!(hdr::FITSHeader, key, comment) set_comment!(hdr, key, comment)
@deprecate getindex(hdu::ImageHDU, I...) read(hdu, I...)