text:
struct StrictChunker
    w_max::Int
end

function pack_stripe(A::SparseMatrixCSC{Tv, Ti}, method::StrictChunker, args...; kwargs...) where {Tv, Ti}
    @inbounds begin
        # matrix notation...
        # i = 1:m rows, j = 1:n columns
        m, n = size(A)

        w_max = method.w_max

        A_pos = A.colptr
        A_idx = A.rowval

        #hst = zeros(Int, m)

        spl = Vector{Int}(undef, n + 1) # Column split locations

        c = A_pos[2] - A_pos[1] #The cardinality of the first column in the part
        j = 1
        K = 0
        spl[1] = 1
        for j′ = 2:n
            c′ = A_pos[j′ + 1] - A_pos[j′] #The cardinality of the candidate column
            w = j′ - j #Current block size
            d = true
            if c == c′ && w != w_max
                l′ = A_pos[j′]
                for l = A_pos[j]:(A_pos[j + 1] - 1)
                    if A_idx[l] != A_idx[l′]
                        d = false
                        break
                    end
                    l′ += 1
                end
            else
                d = false
            end
            if !d
                K += 1
                spl[K + 1] = j′
                j = j′
                c = c′
            end
        end
        K += 1
        spl[K + 1] = n + 1

        resize!(spl, K + 1)

        return SplitPartition{Ti}(K, spl)
    end
end


text:
using .NetworkResilience
using Statistics

function main(t)

	ind = parse(Int64, ARGS[1])

	configs = []
	n = 67
	for j=1:n
		for i=1:n-j
			ne = n - i - j
			push!(configs, [i j ne])
		end
	end
	break1 = length(configs)
	# println(break1)

	#Make the file name
	directory = "austria/"
	indexdict = Dict("67" => ind, "100" => ind-break1)
	subdirectorydict = Dict("67" => "67tri/", "100" => "100tri/")
	ntot = sum(configs[ind])
	thisid = indexdict[string(ntot)]
	thissubdir = subdirectorydict[string(ntot)]
	fname = "data/" * directory * thissubdir * string(thisid) * ".txt"

	#Compute the cascade
	ns = configs[ind][1]
	nd = configs[ind][2]
	ne = configs[ind][3]
	ensemblesize = 300
	tol = 5e-4
	n = ns + nd + ne
	q = 0.0
	I = 1.0
	D = 1.0
	k = 4
	Ptot = 1.0
	κ = 5.0
	ac, fmax, iter = swingcascadebisection(ensemblesize, tol, n, ns, nd, q, k, κ, I, D)

	open(fname, "w") do f
		write(f, "$ns $nd $ne $fmax \n")
		write(f, "$ac $fmax $iter $iter \n")
	end

end

main(ARGS)

ext: jl

max_stars_repo_name: o-smith/GridResilience

max_issues_repo_name: o-smith/GridResilience

max_forks_repo_name: o-smith/GridResilience


text:
using Revise
using DataFrames, Distributions, LinearAlgebra, Finometrics, Dates

struct Universe
  df::DataFrame #source dataframe must be properly formatted with a :date column but can include extra observations
  datebegin::Date
  dateend::Date
  returns::Matrix{MFloat64} #should contain the numerical returns data in the DataFrame
  #NOTE for performance reasons, this is M by N!!!

  μ::Vector{Float64} #we can calc this in advance
  N::Int #number of observations in the sample period
  M::Int #number of return series

  missingmap::Vector{Vector{Int}} # a collection of all non-missing indices in returns
end

#creates an investment universe which will be used to generate portfolio returns
function Universe(df::AbstractDataFrame; datebegin::Date = minimum(df[:date]),
  dateend::Date = maximum(df[:date]))::Universe
  subdf::SubDataFrame = view(df,(df[:date] .≥ datebegin) .& (df[:date] .≤ dateend),:)

  N::Int = size(subdf,1) #get the size in number of observations
  M::Int = size(subdf, 2) - 1 #exclude date column

  xcols::Vector{Symbol} = setdiff(names(df), [:date])
  returns::Matrix{MFloat64} = Matrix{MFloat64}(subdf[xcols])'

  #compute the mean for use later
  μ::Vector{Float64} = (r::Int->mean(skipmissing(returns[r,:]))).(1:M)

  #map the missing values
  missingmap::Vector{Vector{Int}} = Vector{Vector{Float64}}()
  for c ∈ 1:N
    push!(missingmap, collect(1:M)[(!ismissing).(returns[:, c])])
  end

  return Universe(df,datebegin, dateend, returns, μ, N, M, missingmap)
end

#lightweight structure for portfolio math
struct Portfolio
  U::Universe
  wₜ::Vector{Float64}
  R::Vector{Float64}
  pμ::Base.RefValue{Float64}

  #the below vector contains all valid weights, and does not include invalid weights
  #thus the shape of wvec is a ragged array, with each array of length equal to the
  #number of valid values
  wvec::Vector{Vector{Float64}} #store this so we don't have to reallocate
  wvecsize::Vector{Int} #holds the lengths so we don't ahve to call length
end

#update the portfolio object
#innermost loop and likely a performance bottleneck, so highly optimized
function Portfolio!(P, wₜ)::Portfolio

  P.wₜ .= wₜ
  for c ∈ 1:P.U.N #for each column
    for r ∈ 1:P.wvecsize[c] #for each valid value in the column
      P.wvec[c][r] = P.wₜ[P.U.missingmap[c][r]] # get the valud value from the missing map
    end
    adjfactor::Float64 = sum(P.wvec[c])
    P.wvec[c] ./= adjfactor #rescale to account for missing values
  end

  #build the returns (Note: Threads.@threads provides signficant performance gains)
  for c ∈ 1:P.U.N
    P.R[c] = sum(P.wvec[c] .* P.U.returns[P.U.missingmap[c],c])
  end

  #update the pointer
  P.pμ[] = mean(P.R)
  return P
end

#construct a brand new portfolio object
function Portfolio(U, wₜ)::Portfolio
  wvecsize::Vector{Int} = (v->length(v)).(U.missingmap)
  wvec::Vector{Vector{Float64}} = (n->Vector{Float64}(undef, n)).(wvecsize)

  P::Portfolio = Portfolio( #pre-allocate if needed, reduces code reuse
    U,
    Vector{Float64}(undef, U.M),
    Vector{Float64}(undef, U.N),
    Ref(-1.0),
    #Matrix{Float64}(undef, U.M, U.N))
    wvec,
    wvecsize)

  return Portfolio!(P, wₜ)
end

function testFMCov(N=560, M=5000, iter::Int=100)
  cols::Vector{Symbol} = (n->Symbol("X",n)).(1:M)
  crosscor::MFloat64 = 0.25
  stdev::MFloat64 = 0.1
  allmean = .02

  ρ::Matrix{MFloat64} = ones(M,M) .* crosscor
  ρ[diagind(ρ)] .= 1.0 #fix the diagonals
  σ::Vector{MFloat64} = (n->stdev).(1:M)
  Σ::Matrix{MFloat64} = ρ .* stdev^2
  L = cholesky(Σ).L

  #make the test matrix
  testmatrix::Matrix{MFloat64} = rand(Normal(), N, M)

  for r ∈ 1:N
    testmatrix[r,:] .= L * testmatrix[r,:] .+ allmean
  end

  #make missing values
  for c ∈ 1:M
    for r ∈ 1:N
      rand([true, false]) && (testmatrix[r,c] = missing)
    end
  end

  #make the dataframe
  df = DataFrame(testmatrix)
  names!(df, cols)
  df[:date] = Vector{MDate}(missing, N)

  datebegin::Date = Date(2011,11,11)

  df[:date] .= datebegin .+ (Day).(1:N)

  U = Universe(df)
  wₜ = ones(M) ./ M

  P = Portfolio(U,wₜ)

  @time for i ∈ 1:iter
    wₜ .= rand(M)
    Portfolio!(P, wₜ)
  end

  return nothing
end

testFMCov()

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#!/usr/bin/env julia

using PastaQ
using LinearAlgebra
using ITensors
using Plots

function PastaQ.gate(::GateName"R"; theta::Real, phi::Real)
    [
        cos(theta/2)    (-im * exp(-im * phi) * sin(theta/2))
        (-im * exp(im * phi) * sin(theta/2))     cos(theta/2)
    ]
end

function PastaQ.gate(::GateName"M"; Theta::Real)
    [
        cos(Theta)    0    0    (-im * sin(Theta))
        0    cos(Theta)    (-im * sin(Theta))    0
        0    (-im * sin(Theta))    cos(Theta)    0
        (-im * sin(Theta))    0    0    cos(Theta)
    ]
end

function run(N, depth)
    # Random circuit.
    gates = Vector{Tuple}[]

    rdms = []
    for i in 1:N
        push!(rdms, i)
    end

    rdm = rand((rdms))

    for i in 1:depth
        one_qubit_layer = Tuple[]
        two_qubit_layer = Tuple[]
        err_qubit_layer = Tuple[]

        begin (rdm < N)
        idx_err = rdm
        end

        for j in idx_err
            gate = ("E", j)
            push!(err_qubit_layer, gate)
        end

        for j in 1:N
            if j != rdm
                gate = ("R", j, (theta=2pi*rand(), phi=2pi*rand()))
                push!(one_qubit_layer, gate)
            end
        end

        # Alternate start qubit for pairs.
        idx_first = i % 2 + 1

        for j in idx_first:2:(N-1)
            if j != rdm
                gate = ("M", (j, j+1), (Theta=2pi*rand(),))
                push!(two_qubit_layer, gate)
            end
        end

        push!(gates, one_qubit_layer)
        push!(gates, two_qubit_layer)
        push!(gates, err_qubit_layer)

    end

    psi = runcircuit(N, gates)
end


text:
module Showing

using AWS
using AWS.S3

import Base.show, Base.Meta.quot

print_indented(io::IO,s,indent) = print(io," "^indent,s)
show_indented(io::IO,s::String,indent) = print(io,s*"\n")
show_indented(io::IO,::Void,indent) = print(io,"Not set\n")
show_indented(io::IO,i,indent) = println(io,i)

function show_indented(io::IO, s::Array, indent)
    println(io,eltype(s)," Array")
    for x in s
        print_indented(io,"",indent)
        show_indented(io,x,indent+4)
    end
end

function show_indented(io::IO, s::Dict, indent)
    println(io,eltype(s)," Dict")
    for (k,v) in collect(s)
        print_indented(io,"",indent)
        show_indented(io,string(k," => ",v),indent+4)
    end
end

function Base.show(io::IO, awsEnv::AWS.AWSEnv)
    print(io, """AWSEnv:
    aws_id: $("*******$(awsEnv.aws_id[end-2:end])")
    aws_seckey: "*********$(awsEnv.aws_seckey[end-2:end])"
    region: $(awsEnv.region)
    scheme: $(awsEnv.ep_scheme)
    host: $(awsEnv.ep_host)
    path: $(awsEnv.ep_path)
    ver: $(awsEnv.sig_ver)
    timeout: $(awsEnv.timeout)
    dry_run: $(awsEnv.dry_run)
    debug: $(awsEnv.dbg)
""")
end

macro show_func(n, t)
    #n,t,id_sym = eval(n),eval(t),eval(id_sym)
    block = Expr(:block)
    for x in fieldnames(t)
        push!(block.args,quote
            $(esc(:print_indented))(io,$(string(x)*": "),indent)
            $(esc(:show_indented))(io,getfield(r, ($(quot(x)))),indent+4)
        end)
    end
    ret = quote
        function $(esc(:show_indented))(io::IO,r::$(t),indent)
            $(esc(:print))(io,$(string(n)))
            println(io)
            $(block)
        end
        $(esc(:show))(io::IO,r::$(t)) = $(esc(:show_indented))(io,r,0)
    end
    ret
end

# Create default show method for all types
for m in [AWS, AWS.S3]
    for n in names(m)
        t = getfield((m), (n))
        if isa(t,Type) && (t !== AWS.AWSEnv)
            @eval @show_func $n $t
        end
    end
end

end


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

##########################################################################################92

@doc raw"""
    sphere_point_to_uv(point::Point) :: Vec2d

Convert a 3D `point` ``P = (P_x, P_y, P_z)`` on the surface of the unit sphere
into a 2D `Vec2d` using the following spherical coordinates:

```math
u = \frac{\phi}{2\pi} = \frac{\arctan (P_y / P_x)}{2\pi},
    \quad
v = \frac{\theta}{\pi} = \frac{\arccos (P_z)}{\pi}
```

See also: [`Point`](@ref), [`Vec2d`](@ref), [`Sphere`](@ref)
"""
function sphere_point_to_uv(point::Point)
    v = acos(point.z) / π
    u = atan(point.y, point.x) / (2.0 * π)
    u = u>=0 ? u : u + 1.0
    return Vec2d(u,v)
end

@doc raw"""
    plane_point_to_uv(point::Point) :: Vec2d

Convert a 3D `point` ``P = (P_x, P_y, P_z)`` on the surface of the unit plane
into a 2D `Vec2d` using the following periodical coordinates:

```math
u = P_x - \lfloor P_x \rfloor,
    \quad
v = P_y - \lfloor P_y \rfloor,
```

where ``\lfloor \cdot \rfloor`` indicates the rounding down approximation,
in order to guarantee that ``u, v \in [0, 1)``.

See also: [`Point`](@ref), [`Vec2d`](@ref), [`Plane`](@ref)
"""
function plane_point_to_uv(point::Point)
    u = point.x - floor(point.x)
    v = point.y - floor(point.y)
    return Vec2d(u,v)
end

@doc raw"""
    cube_point_to_uv(point::Point) :: Vec2d

Convert a 3D `point` ``P = (P_x, P_y, P_z)`` on the surface of the unit cube
into a 2D `Vec2d` using the following  coordinates:

```math
P_x = \frac{1}{2} \lor P_x = -\frac{1}{2}
\quad \Rightarrow \quad
u = P_y +  \frac{1}{2} \; , \;   v = P_z +  \frac{1}{2}
```
```math
P_y = \frac{1}{2} \lor P_y = -\frac{1}{2}
\quad \Rightarrow \quad
u = P_x +  \frac{1}{2} \; , \;   v = P_z +  \frac{1}{2}
```
```math
P_z = \frac{1}{2} \lor P_z = -\frac{1}{2}
\quad \Rightarrow \quad
u = P_x +  \frac{1}{2} \; , \;   v = P_y +  \frac{1}{2}
```
```math
P_x \neq \frac{1}{2},  -\frac{1}{2} \land
P_y \neq \frac{1}{2},  -\frac{1}{2} \land
P_z \neq \frac{1}{2},  -\frac{1}{2}
\quad \Rightarrow \quad
\mathrm{throw Exception}
```

See also: [`Point`](@ref), [`Vec2d`](@ref), [`Cube`](@ref)
"""
function cube_point_to_uv(point::Point)
    if (point.x ≈ 0.5 || point.x ≈ -0.5)
        u, v  = point.y + 0.5, point.z + 0.5
    elseif (point.y ≈ 0.5 || point.y ≈ -0.5)
        u, v  = point.x + 0.5, point.z + 0.5
    elseif (point.z ≈ 0.5 || point.z ≈ -0.5)
        u, v  = point.x + 0.5, point.y + 0.5
    else
        throw(ArgumentError("the given point do not belong to the unit cube."))
    end

    return Vec2d(u,v)
end

function torus_point_to_uv(point::Point)
    len_point = norm(point)
    u = asin(point.y/len_point) / π
    v = atan(point.z, point.x) / (2.0 * π)
    v>=0 ? nothing : v+= 1.0
    u>=0 ? nothing : u+= 1.0
    return Vec2d(u,v)
end

@doc raw"""
    triangle_point_to_uv(triangle::Triangle, point::Point) :: Vec2d

Return the barycentic coordinates of the given `point` for the input
`triangle`.

If the triangle is made of the vertexes ``(A,B,C)`` (memorized in this order),
then the point ``P`` has coordinates ``(u,v) = (\beta, \gamma)`` such that:
```math
    P(\beta, \gamma) = A + \beta \,(B - A) + \gamma \,(C-A)
```
The analitic resolution of this linear system is:
```math
\begin{aligned}
&\beta = \frac{
            (P_x - A_x)(C_y - A_y) - (P_y - A_y)(C_x - A_x)
        }{
            (B_x - A_x)(C_y - A_y) - (B_y - A_y)(C_x - A_x)
        } \\
&\gamma = \frac{
            (P_x - A_x)(B_y - A_y) - (P_y - A_y)(B_x - A_x)
        }{
            (C_x - A_x)(B_y - A_y) - (C_y - A_y)(B_x - A_x)
        }
\end{aligned}
```

**NOTE**: this function do not check if ``P`` is on the plane defined by ``(A,B,C)``,
neither if ``P`` is inside the triangle made of them!

See also: [`Triangle`](@ref), [`Vec2d`](@ref), [`Point`](@ref)
"""
function triangle_point_to_uv(triangle::Triangle, point::Point)
    A, B, C = Tuple(P for P in triangle.vertexes)
    P = point

    β_num = (P.x - A.x)*(C.y - A.y) - (P.y - A.y)*(C.x - A.x)
    β_den = (B.x - A.x)*(C.y - A.y) - (B.y - A.y)*(C.x - A.x)
    β = β_num/β_den

    γ_num = (P.x - A.x)*(B.y - A.y) - (P.y - A.y)*(B.x - A.x)
    γ_den = (C.x - A.x)*(B.y - A.y) - (C.y - A.y)*(B.x - A.x)
    γ = γ_num/γ_den

    Vec2d(β, γ)
end

##########################################################################################92

@doc raw"""
    sphere_normal(point::Point, ray_dir::Vec) :: Normal

Compute the `Normal` of a unit sphere.

The normal is computed for the given `Point` ``point = (P_x, P_y, P_z)``
(with ``\sqrt{P_x^2 + P_y^2 + P_z^2}=1``) on the
surface of the sphere, and it is chosen so that it is always in the opposite
direction with respect to the given `Vec` `ray_dir`.

See also: [`Point`](@ref), [`Ray`](@ref), [`Normal`](@ref), [`Sphere`](@ref)
"""
function sphere_normal(point::Point, ray_dir::Vec)
    result = Normal(point.x, point.y, point.z)
    Vec(point) ⋅ ray_dir < 0.0 ? nothing : result = -result
    return result
end

@doc raw"""
    plane_normal(point::Point, ray_dir::Vec) :: Normal

Compute the `Normal` of a unit plane.

The normal is computed for the given `point` ``P = (P_x, P_y, 0)`` on the
surface of the plane, and it is chosen so that it is always in the opposite
direction with respect to the given `Vec` `ray_dir`.

See also: [`Point`](@ref), [`Ray`](@ref), [`Normal`](@ref), [`Plane`](@ref)
"""
function plane_normal(point::Point, ray_dir::Vec)
    result = Normal(0., 0., 1.)
    Vec(0., 0., 1.) ⋅ ray_dir < 0.0 ? nothing : result = -result
    return result
end

@doc raw"""
    cube_normal(point::Point, ray_dir::Vec) :: Normal

Compute the `Normal` of a unit cube.

The normal is computed for the given `point` on the
surface of the cube, and it is chosen so that it is always in the opposite
direction with respect to the given `Vec` `ray_dir`.

See also: [`Point`](@ref), [`Ray`](@ref), [`Normal`](@ref), [`Cube`](@ref)
"""
function cube_normal(point::Point, ray_dir::Vec)
    if (point.x ≈ 0.5 || point.x ≈ -0.5)
        result = Normal(1., 0., 0.)
    elseif (point.y ≈ 0.5 || point.y ≈ -0.5)
        result = Normal(0., 1., 0.)
    elseif (point.z ≈ 0.5 || point.z ≈ -0.5)
        result = Normal(0., 0., 1.)
    else
        throw(ArgumentError("the given point do not belong to the unit cube."))
    end

    result ⋅ ray_dir < 0.0 ? nothing : result = -result
    return result
end

"""
    torus_normal(p::Point, ray_dir::Vec, R::Float64) -> Normal

Compite the [`Normal`](@ref) of a torus

The normal is computed for [`Point`](@ref) (a point on the surface of the
torus), and it is chosen so that it is always in the opposite
direction with respect to `ray_dir` ([`Vec`](@ref)).
"""
function torus_normal(p::Point, ray_dir::Vec, R::Float64)
    R_z = copysign(R / √(1+(p.x/p.z)^2), p.z)
    R_x = copysign(p.x / p.z * R_z, p.x)
    R_p = Vec(R_x, 0, R_z)
    result = Normal(Vec(p - R_p))
    result ⋅ ray_dir < 0.0 ? nothing : result = -result
    return result
end

@doc raw"""
    triangle_normal(triangle::Triangle, ray_dir::Vec) :: Normal

Compute the `Normal` of a given triangle.

The normal for a triangle with vertexes ``(A, B, C)`` is computed as follows:
```math
    n = \pm (B-A) \times (C-A)
```
where the sign is chosen so that it is always in the opposite
direction with respect to the given `ray_dir`.

See also: [`Point`](@ref), [`Ray`](@ref), [`Normal`](@ref), [`Triangle`](@ref)
"""
function triangle_normal(triangle::Triangle, ray_dir::Vec)
    result = (triangle.vertexes[2] -  triangle.vertexes[1]) ×
                (triangle.vertexes[3] -  triangle.vertexes[1])
    result ⋅ ray_dir < 0.0 ? nothing : result = -result
    return Normal(result)
end

@doc raw"""
    triangle_barycenter(triangle::Triangle) :: Point

Return the barycenter of the given `triangle`.

For a triangle with vertexes ``(A, B, C)``, the barycenter
is ``M``:
```math
    M = \frac{A + B + C}{3}
```

See also: [`Triangle`](@ref), [`Point`](@ref)
"""
function triangle_barycenter(triangle::Triangle)
    A, B, C = Tuple(P for P in triangle.vertexes)
    result = Point(A.x+B.x+C.x, A.y+B.y+C.y, A.z+B.z+C.z)*1/3
    return result
end

##########################################################################################92

"""
    ray_intersection(shape::Shape, ray::Ray) :: Union{HitRecord, Nothing}

Compute the intersection between a `Ray` and a `Shape`.

See also: [`Ray`](@ref), [`Shape`](@ref)
"""
function ray_intersection(shape::Shape, ray::Ray)
    return ErrorException(
            "ray_intersection is an abstract method"*
            "and cannot be called directly"
            )
end

"""
    ray_intersection(AABB::AABB, ray::Ray) :: Bool

Check if the `ray` intersects the `AABB`.
Return `true` if intersection occurs, `false` otherwise.

See also: [`Ray`](@ref), [`AABB`](@ref), [`HitRecord`](@ref)
"""
function ray_intersection(AABB::AABB, ray::Ray)
    (tmin, tmax) = Tuple( sort( [
                        (AABB.m.x - ray.origin.x) / ray.dir.x,
                        (AABB.M.x - ray.origin.x) / ray.dir.x
                    ]) )
    (tymin, tymax) = Tuple( sort( [
                        (AABB.m.y - ray.origin.y) / ray.dir.y,
                        (AABB.M.y - ray.origin.y) / ray.dir.y
                    ]) )

    ((tmin > tymax) || (tymin > tmax)) && (return false)

    (tymin > tmin) && (tmin = tymin)
    (tymax < tmax) && (tmax = tymax)

    (tzmin, tzmax) = Tuple( sort( [
                        (AABB.m.z - ray.origin.z) / ray.dir.z,
                        (AABB.M.z - ray.origin.z) / ray.dir.z
                    ]) )

    ((tmin > tzmax) || (tzmin > tmax)) && (return false)

    (tzmin > tmin) && (tmin = tzmin)
    (tzmax < tmax) && (tmax = tzmax)

    if (ray.tmin ≤ tmin ≤ ray.tmax) || ( ray.tmin ≤ tmax ≤ ray.tmax)
        return true
    else
        return false
    end
end

"""
    ray_intersection(sphere::Sphere, ray::Ray) :: Union{HitRecord, Nothing}

Check if the `ray` intersects the `sphere`.
Return a `HitRecord`, or `nothing` if no intersection is found.

See also: [`Ray`](@ref), [`Sphere`](@ref), [`HitRecord`](@ref)
"""
function ray_intersection(sphere::Sphere, ray::Ray)

    (ray_intersection(sphere.AABB, ray) == true) || (return nothing)

    inv_ray = inverse(sphere.T) * ray
    origin_vec = Vec(inv_ray.origin)

    a = squared_norm(inv_ray.dir)
    b = 2.0 * origin_vec ⋅ inv_ray.dir
    c = squared_norm(origin_vec) - 1.0
    Δ = b * b - 4.0 * a * c

    (Δ > 0.0) || (return nothing)

    tmin = (-b - √Δ) / (2.0 * a)
    tmax = (-b + √Δ) / (2.0 * a)

    if (tmin > inv_ray.tmin) && (tmin < inv_ray.tmax)
        first_hit_t = tmin
    elseif (tmax > inv_ray.tmin) && (tmax < inv_ray.tmax)
        first_hit_t = tmax
    else
        return nothing
    end

    hit_point = at(inv_ray, first_hit_t)

    return HitRecord(
        sphere.T * hit_point,
        sphere.T * sphere_normal(hit_point, inv_ray.dir),
        sphere_point_to_uv(hit_point),
        first_hit_t,
        ray,
        sphere
    )
end

"""
    ray_intersection(plane::Plane, ray::Ray) :: Union{HitRecord, Nothing}

Check if the `ray` intersects the `plane`.
Return a `HitRecord`, or `nothing` if no intersection is found.

See also: [`Ray`](@ref), [`Plane`](@ref), [`HitRecord`](@ref)
"""
function ray_intersection(plane::Plane, ray::Ray)
    inv_ray = inverse(plane.T) * ray

    !(inv_ray.dir.z ≈ 0.) || (return nothing)

    hit_t = - inv_ray.origin.z / inv_ray.dir.z

    ( (hit_t > inv_ray.tmin) && (hit_t < inv_ray.tmax) ) || (return nothing)

    hit_point = at(inv_ray, hit_t)

    return HitRecord(
        plane.T * hit_point,
        plane.T * plane_normal(hit_point, inv_ray.dir),
        plane_point_to_uv(hit_point),
        hit_t,
        ray,
        plane
    )
end

"""
    ray_intersection(cube::Cube, ray::Ray) :: Union{HitRecord, Nothing}

Check if the `ray` intersects the `cube`.
Return a `HitRecord`, or `nothing` if no intersection is found.

The implementation is only a long boring list of `if`-`else` block statements,
and may have to be optimized.

See also: [`Ray`](@ref), [`Cube`](@ref), [`HitRecord`](@ref)
"""
function ray_intersection(cube::Cube, ray::Ray)
    (ray_intersection(cube.AABB, ray) == true) || (return nothing)

    inv_ray = inverse(cube.T) * ray
    d = inv_ray.dir
    O = inv_ray.origin

    (tmin, tmax) = Tuple( sort( [ (-0.5 - O.x) / d.x, (0.5 - O.x) / d.x ]) )
    (tymin, tymax) = Tuple( sort( [ (-0.5 - O.y) / d.y, (0.5 - O.y) / d.y ]) )

    ((tmin > tymax) || (tymin > tmax)) && (return false)

    (tymin > tmin) && (tmin = tymin)
    (tymax < tmax) && (tmax = tymax)

    (tzmin, tzmax) = Tuple( sort( [ (-0.5 - O.z) / d.z, (0.5 - O.z) / d.z ]) )

    ((tmin > tzmax) || (tzmin > tmax)) && (return false)

    (tzmin > tmin) && (tmin = tzmin)
    (tzmax < tmax) && (tmax = tzmax)

    if (inv_ray.tmin ≤ tmin ≤ inv_ray.tmax)
        hit_point = at(inv_ray, tmin)
        return HitRecord(
            cube.T * hit_point,
            cube.T * cube_normal(hit_point, inv_ray.dir),
            cube_point_to_uv(hit_point),
            tmin,
            ray,
            cube
        )
    elseif ( inv_ray.tmin ≤ tmax ≤ inv_ray.tmax)
        hit_point = at(inv_ray, tmax)
        return HitRecord(
            cube.T * hit_point,
            cube.T * cube_normal(hit_point, inv_ray.dir),
            cube_point_to_uv(hit_point),
            tmax,
            ray,
            cube
        )
    else
        return nothing
    end
end

function ray_intersection(torus::Torus, ray::Ray)
    inv_ray = inverse(torus.T) * ray

    d = normalize(inv_ray.dir)
    o = inv_ray.origin
    norm²_d = squared_norm(d)
    norm²_o = squared_norm(Vec(o))
    r = torus.r
    R = torus.R
    c4 = norm²_d^2
    c3 = 4 * norm²_d * (Vec(o) ⋅ d)
    c2 = 2 * norm²_d * (norm²_o - r^2 - R^2) + 4 * (Vec(o) ⋅ d)^2 + 4 * R^2 * (d.y)^2
    c1 = 4 * (norm²_o - r^2 - R^2) *  (Vec(o) ⋅ d) + 8 * R^2 * o.y * d.y
    c0 = (norm²_o - r^2 - R^2)^2 - 4 * R^2 * (r^2 - (o.y)^2)

    t_ints = roots(Polynomial([c0, c1, c2, c3, c4]))

    hit_t = Union{Float64, Nothing}
    hit_t = nothing
    for i in t_ints
        if typeof(i) == Float64
            if (i > inv_ray.tmin) && (i < inv_ray.tmax)
                hit_t = i
            end
        end
    end

    if hit_t == nothing
        return nothing
    end

    hit_point = at(inv_ray, hit_t)

    return HitRecord(
        torus.T * hit_point,
        torus.T * torus_normal(hit_point, inv_ray.dir, torus.R),
        torus_point_to_uv(hit_point), # manca la funzione
        hit_t,
        ray,
        torus
    )
end

@doc raw"""
    ray_intersection(triangle::Triangle, ray::Ray) :: Union{HitRecord, Nothing}

Check if the `ray` intersects the `triangle`.
Return a `HitRecord`, or `nothing` if no intersection is found.

For a triangle with vertexes ``(A, B, C)`` and a ray defined with the
simple equation ``r(t) = O + t \, \vec{d}``, the coordinates ``(u,v) = (\beta, \gamma)``
and the `t` value of intersection are obtained solving this linear system:
```math

\begin{bmatrix}
    B_x-A_x & C_x-A_x & -d_x \\
    B_y-A_y & C_y-A_y & -d_y \\
    B_z-A_z & C_z-A_z & -d_z
\end{bmatrix}

\begin{bmatrix}
u \\
v \\
t
\end{bmatrix}
=
\begin{bmatrix}
O_x - A_x\\
O_z - A_z\\
O_z - A_z
\end{bmatrix}
```

See also: [`Ray`](@ref), [`Triangle`](@ref), [`HitRecord`](@ref)
"""
function ray_intersection(triangle::Triangle, ray::Ray)

    A, B, C = Tuple(P for P in triangle.vertexes)
    m = [ray.origin.x-A.x;  ray.origin.y-A.y; ray.origin.z-A.z]
    M = [
        B.x-A.x C.x-A.x -ray.dir.x ;
        B.y-A.y C.y-A.y -ray.dir.y ;
        B.z-A.z C.z-A.z -ray.dir.z ;
    ]

    try
        w = transpose(m) / transpose(M)
        u, v, hit_t = Tuple(x for x in w)
        ( ray.tmin < hit_t < ray.tmax ) || (return nothing)
        ( (u>0.) && (v>0.) && (1-u-v>0.) ) || (return nothing)
        hit_point = at(ray, hit_t)
        return HitRecord(
            hit_point,
            triangle_normal(triangle, ray.dir),
            Vec2d(u,v),
            hit_t,
            ray,
            triangle
        )
    catch Excep
        return nothing
    end
end

"""
    ray_intersection(world::World, ray::Ray) :: Union{HitRecord, Nothing}

Determine whether the `ray` intersects any of the objects of the given `world`.
Return a `HitRecord`, or `nothing` if no intersection is found.

See also: [`Ray`](@ref), [`World`](@ref), [`HitRecord`](@ref)
"""
function ray_intersection(world::World, ray::Ray)
    closest = nothing

    for shape in world.shapes
        intersection = ray_intersection(shape, ray)

        # The ray missed this shape, skip to the next one
        !(isnothing(intersection)) || continue

        # There was a hit, and it was closer than any other hit found before
        ( isnothing(closest) || (intersection.t < closest.t) ) &&  (closest = intersection)
    end

    return closest
end

##########################################################################################92

"""
    add_shape!(world::World, shape::Shape)

Append a new `shape` to the given `world`.

See also: [`Shape`](@ref), [`World`](@ref)
"""
function add_shape!(world::World, S::Shape)
    push!(world.shapes, S)
    return nothing
end

"""
    add_light!(world::World, pointlight::PointLight)

Append a new `pointlight` to the given `world`.

See also: [`PointLight`](@ref), [`World`](@ref)
"""
function add_light!(world::World, pointlight::PointLight)
    push!(world.point_lights, pointlight)
    return nothing
end

##########################################################################################92

"""
    quick_ray_intersection(shape::Shape, ray::Ray) :: Bool

Quickly determine whether the `ray` hits the `shape` or not.

See also: [`Shape`](@ref), [`Ray`](@ref)
"""
function quick_ray_intersection(shape::Shape, ray::Ray)
    return ErrorException(
            "quick_ray_intersection is an abstract method"*
            "and cannot be called directlly")
end

"""
    quick_ray_intersection(sphere::Sphere, ray::Ray) :: Bool

Quickly checks if the `ray` intersects the `sphere` or not.

See also: [`Sphere`](@ref), [`Ray`](@ref)
"""
function quick_ray_intersection(sphere::Sphere, ray::Ray)
    inv_ray = inverse(sphere.T) * ray
    origin_vec = Vec(inv_ray.origin)

    a = squared_norm(inv_ray.dir)
    b = 2.0 * origin_vec ⋅ inv_ray.dir
    c = squared_norm(origin_vec) - 1.0
    Δ = b * b - 4.0 * a * c

    (Δ > 0.0) || (return false)

    tmin = (-b - √Δ) / (2.0 * a)
    tmax = (-b + √Δ) / (2.0 * a)

    return ((inv_ray.tmin < tmin < inv_ray.tmax) || (inv_ray.tmin < tmax < inv_ray.tmax))
end

"""
    quick_ray_intersection(plane::Plane, ray::Ray) :: Bool

Quickly checks if the `ray` intersects the `plane` or not.

See also: [`Plane`](@ref), [`Ray`](@ref)
"""
function quick_ray_intersection(plane::Plane, ray::Ray)
    inv_ray = inverse(plane.T) * ray
    !(inv_ray.dir.z ≈ 0.) || (return false)

    t = -inv_ray.origin.z / inv_ray.dir.z
    return (inv_ray.tmin < t < inv_ray.tmax)
end

"""
    is_point_visible(
            world::World,
            point::Point,
            observer_pos::Point
            ) :: Bool

Return `true` if the straight line connecting `observer_pos` to `point`
do not intersect any of the shapes of `world` between the two points,
otherwise return `false`.

See also: [`World`](@ref), [`Point`](@ref)
"""
function is_point_visible(world::World, point::Point, observer_pos::Point)
    direction = point - observer_pos
    dir_norm = norm(direction)

    ray = Ray(observer_pos, direction, 1e-2 / dir_norm, 1.0)
    for shape in world.shapes
        if (quick_ray_intersection(shape, ray) == true) && (shape.flag_pointlight==false)
            return false
        end
    end

    return true
end


text:
"""
    A = JopLeakyIntegration(spc; α=0.99, dim=1)

`A` is Leaky integration (i.e. a running sum) along dimension `dim`
of `spc::JetSpace{T,N}`.  The result is scaled by `α`.
"""

#=
From stanford: http://sep.stanford.edu/sep/prof/geelib/leakint.lop
    tt = 0.0

    if (adj)
        do i =  n, 1, -1 {
            tt = rho*tt + yy(i)
            xx(i) += tt
        }
    else
        do i =  1, n {
            tt = rho*tt + xx(i)
            yy(i) += tt
        }
    }
=#

function JopLeakyIntegration(spc::JetSpace{T}; α=0.99, dim=1) where {T}
    JopLn(dom = spc, rng = spc, df! = JopLeakyIntegration_df!, df′! = JopLeakyIntegration_df′!,
        s = (dim=dim, α=T(α)))
end

export JopLeakyIntegration

function JopLeakyIntegration_df!(d::AbstractArray{T}, m::AbstractArray{T}; dim, α, kwargs...) where {T}
    n1 = size(d,1)
    n2 = div(length(d),n1)
    _d = reshape(d,n1,n2)
    _m = reshape(m,n1,n2)
    for k2 = 1:n2
        tmp = T(0)
        for k1 = 1:n1
            tmp = α * tmp + _m[k1,k2]
            _d[k1,k2] = tmp
        end
    end
    d
end

function JopLeakyIntegration_df′!(m::AbstractArray{T}, d::AbstractArray{T}; dim, α, kwargs...) where {T}
    n1 = size(d,1)
    n2 = div(length(d),n1)
    _d = reshape(d,n1,n2)
    _m = reshape(m,n1,n2)
    for k2 = 1:n2
        tmp = T(0)
        for k1 = n1:-1:1
            tmp = α * tmp + _d[k1,k2]
            _m[k1,k2] = tmp
        end
    end
    m
end


text:
using FractionalDiffEq

function ϕ(x)
    if x == 0
        return 19.00001
    else
        return 19.0
    end
end

function f(t, y, ϕ)
    return 3.5*y*(1-ϕ/19)
end
t0=0; T=2
tspan = (t0, T); h=0.5
prob = FDDEProblem(f, ϕ, 0.97, 0.8, tspan)

result = solve(prob, h, DelayPI())
using Plots
tspan = collect(0:h:T)
plot(tspan, result[:])


text:
#= -------------------------------------------------------------------------
# @file problema1_sin_disipador_sundials.jl
#
# @date 11/27/17 18:07:28
# @author Martin Noblia
# @email mnoblia@disroot.org
#
# @brief
#
# @detail
#
#  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
---------------------------------------------------------------------------=#
import Sundials
using PyPlot

const  C_j   =  0.5
const  C_c   =  6.8
const  R_jc  =  1.0
const  R_ca  =  30
const  T_a   =  30
P_d  = (150 - 30) / (1+30)
a_11 = (-1 / C_c) * ( ( 1 / R_jc ) + ( 1 / R_ca ))
a_12 = (1 / (R_jc * C_c))
a_21 = (1 / (C_j * R_jc))
a_22 = (-1/ (C_j * R_jc))

step(t) = (t >=0.0) ? P_d: 0.0
pulse(t) = (P_d * sin(2.0*π*100.0*t) + P_d) / 2.0
"""
docs
"""
function system(t, T, dT)
   dT[1] = a_11 * T[1] + a_12 * T[2] + (T_a / (R_ca * C_c))
   dT[2] = a_21 * T[1] + a_22 * T[2] + (pulse(t) / C_j)
end

T₀ = [0.0, 0.0] # initial condition
t = linspace(0, 1300, 1000000) # time vector
sol = Sundials.cvode(system, T₀, collect(t))
# plot the results
plot(t, sol[:,1], label=L"T_c")
plot(t, sol[:,2], label=L"T_j")
xlabel(L"t")
ylabel(L"T")
legend()


text:
using  Languages, TextAnalysis, Flux, PyPlot, Statistics

Arr = ["well done",
     "good work",
     "great effort",
     "nice work",
     "excellent",
     "weak",
     "poor effort",
     "not good",
     "poor work",
     "could have done better"]

# positve or negative sentiment to each 'document' string
y = [true true true true true false false false false false]

docs=[]
for i in 1:length(Arr)
    push!(docs, StringDocument(Arr[i]))
end

crps=Corpus(docs)
update_lexicon!(crps)
doc_term_matrix=DocumentTermMatrix(crps)
word_dict=doc_term_matrix.column_indices

tk_idx(s) = haskey(word_dict, s) ? i=word_dict[s] : i=0

function pad_corpus(c, pad_size)
    M=[]
    for doc in 1:length(c)
        tks = tokens(c[doc])
        if length(tks)>=pad_size
            tk_indexes=[tk_idx(w) for w in tks[1:pad_size]]
        end
        if length(tks)<pad_size
            tk_indexes=zeros(Int64,pad_size-length(tks))
            tk_indexes=vcat(tk_indexes, [tk_idx(w) for w in tks])
        end
        doc==1 ? M=tk_indexes' : M=vcat(M, tk_indexes')
    end
    return M
end

pad_size =4
doc_pad_size=4
padded_docs = pad_corpus(crps, pad_size)
x = padded_docs'
data = [(x, y)]

N = size(padded_docs,1)  #Number of documents (10)
max_features = 8
vocab_size = 20

embedding_matrix=Flux.glorot_normal(max_features, vocab_size)

m = Chain(x -> embedding_matrix * Flux.onehotbatch(reshape(x, doc_pad_size*size(x,2)), 0:vocab_size-1),
          x -> reshape(x, max_features, doc_pad_size, trunc(Int64(size(x,2)/doc_pad_size))),
          x -> mean(x, dims=2),
          x -> reshape(x, max_features, :),
          Dense(max_features, 2),
          softmax
)

loss_h=[]
accuracy_train=[]
accuracy(x, y) = mean(x .== y)

loss(x, y) = sum(Flux.binarycrossentropy.(m(x), y))
optimizer = Flux.Descent(0.01)

for epoch in 1:400
    Flux.train!(loss, Flux.params(m), data, optimizer)
    println(loss(x, y).data, " ", accuracy(m(x).>0.5,y))
    push!(loss_h, loss(x, y).data)
    push!(accuracy_train, accuracy(m(x).>0.5,y))
end

println(m(x).>0.5)
accuracy(m(x).>0.5,y)


text:
#=
Created on Tuesday 14 July 2020
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Plotting utility for OTSIN
=#

using Plots, RecipesBase, LaTeXStrings

@userplot MarginalsHeatmap

@recipe function f(h::MarginalsHeatmap; cmap=cgrad([:white, :blue]), cola="#2a9d8f", colb="#f4a261", showaxis=false)
    if length(h.args) != 1 || !(typeof(h.args[1]) <: AbstractMatrix)
        error("Marginals Heatmap should be given a matrix.  Got: $(typeof(h.args))")
    end
    P, = h.args
    a, b = marginals(P)

    # set up the subplots
    legend := false
    link := :both
    framestyle := [:none :axes :none]
    grid := false
    layout := @layout [tophist           _
                       hist2d{0.9w,0.9h} righthist]

    # main histogram2d
    @series begin
        seriescolor := cmap
        seriestype := :heatmap
        subplot := 2
        showaxis := showaxis
        P
    end

    # these are common to both marginal histograms
    fillalpha := 0.8
    linealpha := 0.3
    seriestype := :bar

    # upper histogram
    @series begin
        fillcolor := cola
        subplot := 3
        orientation := :h
        title := L"$\mathbf{a}$"
        a
    end

    # right histogram
    @series begin
        fillcolor := colb
        subplot := 1
        title := L"$\mathbf{b}$"
        b
    end
end

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
#=
Continuation of data, to e.g. fill missing values
=#
using SignalDecomposition
export sinusoidal_continuation

_numbertype(T) = _numbertype(eltype(T))
_numbertype(::Type{Union{T, Missing}}) where T = T
_numbertype(::Type{T}) where {T <: Real} = T

# TODO: perhaps it is worth considering a version that makes true time in seconds
# instead of months, which in the future would allow support of diurnal cycles
# TODO: YES! Make true time be seconds!

"""
    sinusoidal_continuation(T::ClimArray, fs = [1, 2]; Tmin = -Inf, Tmax = Inf)
Fill-in the missing values of spatiotemporal field `T`, by fitting sinusoidals
to the non-missing values, and then using the fitted sinusoidals for the missing values.

Frequencies are given per year (frequency 2 means 1/2 of a year).

`Tmin, Tmax` limits are used to clamp the result into this range (no clamping in the
default case).
"""
function sinusoidal_continuation(T, frequencies = [1.0, 2.0]; Tmin = -Inf, Tmax = Inf)
    E = _numbertype(T)
    lpv = Sinusoidal(E.(frequencies ./ DAYS_IN_ORBIT))
    fullT = copy(T)
    # TODO: this must be extended to a general "true time" function
    truetime = realtime_days(dims(T, Time).val, E)
    for i in otheridxs(T, Time)
        x = T[i...]
        any(ismissing, x) || continue # this timeseries needs no correction
        mi = findall(!ismissing, x)
        x0 = Array{Float32}(x[mi])
        t0 = truetime[mi]
        sea, res = decompose(t0, x0, lpv)
        c = cosines(truetime, lpv)
        imi = findall(ismissing, x)
        x[imi] .= c[imi]
        fullT[i...] .= x
    end
    if !isinf(Tmin) && !isinf(Tmin)
        fullT = clamp.(fullT, Tmin, Tmax)
    end
    return ClimArray(nomissing(fullT.data), T.dims; name = T.name, attrib = T.attrib)
end

function cosines(t, lpv::Sinusoidal)
    s = fill(lpv.A[1], length(t))
    for i in 2:length(lpv.A)
        @. s += lpv.A[i] * cos(2π * lpv.fs[i-1] * t + lpv.φ[i])
    end
    return s
end


text:
#LIKELIHOOD SCORING FUNCS
function IPM_likelihood(sources::AbstractVector{<:Tuple{<:AbstractMatrix{<:AbstractFloat},<:Integer}}, observations::AbstractMatrix{<:Integer}, obs_lengths::AbstractVector{<:Integer}, bg_scores::AbstractArray{<:AbstractFloat}, mix::BitMatrix, revcomp::Bool=REVCOMP, returncache::Bool=false, cache::AbstractVector{<:AbstractFloat}=zeros(0), clean::AbstractVector{<:Bool}=Vector(falses(size(observations)[2])))
    source_wmls=[size(source[1])[1] for source in sources]
    O = size(bg_scores)[2]
    source_stops=[obsl-wml+1 for wml in source_wmls, obsl in obs_lengths] #stop scannng th source across the observation as the source reaches the end
    L=maximum(obs_lengths)+1

    obs_src_idxs=mix_pull_idxs(mix) #get vectors of sources emitting in each obs

    revcomp ? (srcs=[cat(source[1],revcomp_pwm(source[1]),dims=3) for source in sources]; motif_expectations = [((MOTIF_EXPECT/2)/obsl) for obsl in obs_lengths]; mat_dim=2) : (srcs=[source[1] for source in sources]; ; motif_expectations = [(MOTIF_EXPECT/obsl) for obsl in obs_lengths]; mat_dim=1) #setup appropriate reverse complemented sources if necessary and set log_motif_expectation-nMica has 0.5 per base for including the reverse complement, 1 otherwise

    lme_vec=zeros(length(sources))

    obs_lhs=Vector{Vector{Float64}}() #setup likelihood vecs for threaded operation-reassembled on return
    nt=Threads.nthreads()
    for t in 1:nt-1
        push!(obs_lhs,zeros(Int(floor(O/nt))))
    end
    push!(obs_lhs, zeros(Int(floor(O/nt)+(O%nt))))

    Threads.@threads for t in 1:nt
        revcomp && (weavevec=zeros(3))
        revcomp ? (score_mat=zeros(maximum(source_stops),2)) : (score_mat=zeros(maximum(source_stops)))
        opt = floor(O/nt) #obs per thread
        score_matrices=Vector{typeof(score_mat)}(undef, length(sources)) #preallocate
        osi_emitting=Vector{Int64}() #preallocate
        lh_vec = zeros(L) #preallocated likelihood vector is one position (0 initialiser) longer than the longest obs

        for i in 1:Int(opt+(t==nt)*(O%nt))
            o=Int(i+(t-1)*opt)
            if clean[o]
                obs_lhs[t][i]=cache[o]
            else
                obsl = obs_lengths[o]
                oidxs=obs_src_idxs[o]
                mixwmls=source_wmls[oidxs]

                obs_cardinality = length(oidxs) #the more sources, the greater the cardinality_penalty
                if obs_cardinality > 0
                    revcomp ? score_sources_ds!(score_mat, score_matrices, view(srcs,oidxs), view(observations,:,o), view(source_stops,oidxs,o)) :
                         score_sources_ss!(score_mat, score_matrices, view(srcs,oidxs), view(observations,:,o), view(source_stops,oidxs,o)) #get scores for this observation

                    lme_vec.=motif_expectations[o]
                    penalty_sum = sum(lme_vec[1:obs_cardinality])
                    penalty_sum > 1. && (penalty_sum=1.)
                    cardinality_penalty=log(1.0-penalty_sum)
                else
                    cardinality_penalty=0.0
                end

                revcomp ? (obs_lhs[t][i]=weave_scores_ds!(weavevec, lh_vec, obsl, view(bg_scores,:,o), score_matrices, oidxs, mixwmls, log(motif_expectations[o]), cardinality_penalty, osi_emitting)) :
                    (obs_lhs[t][i]=weave_scores_ss!(lh_vec, obsl, view(bg_scores,:,o), score_matrices, oidxs, mixwmls, log(motif_expectations[o]), cardinality_penalty, osi_emitting))

                empty!(osi_emitting)

            end
        end
    end

    returncache ? (return lps([lps(obs_lhs[t]) for t in 1:nt]), vcat(obs_lhs...)) : (return lps([lps(obs_lhs[t]) for t in 1:nt]))
end
                @inline function mix_pull_idxs(A::AbstractArray{Bool})
                    n=count(A)
                    S=[Vector{Int64}() for o in 1:size(A,1)]
                    cnt=1
                    for (i,a) in pairs(A)
                        if a
                            push!(S[i[1]],i[2])
                            cnt+=1
                        end
                    end
                    return S
                end

                @inline function revcomp_pwm(pwm::AbstractMatrix{<:AbstractFloat}) #in order to find a motif on the reverse strand, we scan the forward strand with the reverse complement of the pwm, reordered 3' to 5', so that eg. an PWM for an ATG motif would become one for a CAT motif
                    return pwm[end:-1:1,end:-1:1]
                end

                @inline function score_sources_ds!(score_mat, score_matrices::AbstractVector{<:AbstractMatrix{<:AbstractFloat}}, sources, observation::AbstractVector{<:Integer}, source_stops)
                    for (s,source) in enumerate(sources)
                        for t in 1:source_stops[s]
                            for position in 1:size(source,1)
                                score_loc = t+position-1 #score_loc is the position of the obs to be scored by PWM
                                score_mat[t,1] += source[position,observation[score_loc],1] #add the appropriate log PWM value from the source to the score
                                score_mat[t,2] += source[position,observation[score_loc],2] #add the appropriate log PWM value from the source to the score

                            end
                        end
                        score_matrices[s]=score_mat[1:source_stops[s],:] #copy score matrix to vector
                        score_mat.=0. #reset score matrix
                    end
                end

                @inline function score_sources_ss!(score_mat, score_matrices::AbstractVector{<:AbstractArray{<:AbstractFloat}}, sources, observation::AbstractVector{<:Integer}, source_stops)
                    for (s,source) in enumerate(sources)
                        for t in 1:source_stops[s]
                            for position in 1:size(source,1)
                                score_loc = t+position-1 #score_loc is the position of the obs to be scored by PWM
                                score_mat[t] += source[position,observation[score_loc]] #add the appropriate log PWM value from the source to the score

                            end
                        end
                        score_matrices[s]=score_mat[1:source_stops[s]] #copy score matrix to vector
                        score_mat.=0. #reset score matrix
                    end
                end

                @inline function weave_scores_ds!(weavevec, lh_vec, obsl::Integer, bg_scores::SubArray, score_mat::AbstractVector{<:AbstractMatrix{<:AbstractFloat}}, obs_source_indices::AbstractVector{<:Integer}, source_wmls::AbstractVector{<:Integer}, log_motif_expectation::AbstractFloat, cardinality_penalty::AbstractFloat, osi_emitting)
                    for i in 2:obsl+1 #i=1 is ithe lh_vec initializing 0, i=2 is the score of the first background position (ie t=1)
                        t=i-1
                        score = lps(lh_vec[i-1], bg_scores[t], cardinality_penalty)

                        #logic: all observations are scored from t=wml to the end of the obs-therefore check at each position for new sources to add (indexed by vector position to retrieve source wml and score matrix)
                        if length(osi_emitting)<length(obs_source_indices)
                            for n in 1:length(obs_source_indices)
                                if !(n in osi_emitting)
                                    t>= source_wmls[n] && (push!(osi_emitting,n))
                                end
                            end
                        end

                        for n in osi_emitting
                            wml = source_wmls[n]
                            from_score = lh_vec[i-wml+1] #score at the first position of the PWM
                            score_array = score_mat[n] #get the source score matrix
                            score_idx = t - wml + 1 #translate t to score_array idx for emission score
                            f_emit_score = score_array[score_idx,1] #emission score at the last position of the PWM
                            r_emit_score = score_array[score_idx,2]

                            weavevec .= score, lps(from_score, f_emit_score, log_motif_expectation), lps(from_score, r_emit_score, log_motif_expectation)

                            score=logsumexp(weavevec)
                        end
                        lh_vec[i] = score
                    end
                    return lh_vec[obsl+1]
                end

                @inline function weave_scores_ss!(lh_vec, obsl::Integer, bg_scores::SubArray, score_mat::AbstractVector{<:AbstractVector{<:AbstractFloat}}, obs_source_indices::AbstractVector{<:Integer}, source_wmls::AbstractVector{<:Integer}, log_motif_expectation::AbstractFloat, cardinality_penalty::AbstractFloat, osi_emitting)
                    for i in 2:obsl+1 #i=1 is ithe lh_vec initializing 0, i=2 is the score of the first background position (ie t=1)
                        t=i-1
                        score = lps(lh_vec[i-1], bg_scores[t], cardinality_penalty)
                        #logic: all observations are scored from t=wml to the end of the obs-therefore check at each position for new sources to add (indexed by vector position to retrieve source wml and score matrix)
                        if length(osi_emitting)<length(obs_source_indices)
                            for n in 1:length(obs_source_indices)
                                if !(n in osi_emitting)
                                    t>= source_wmls[n] && (push!(osi_emitting,n))
                                end
                            end
                        end

                        for n in osi_emitting
                            wml = source_wmls[n]
                            from_score = lh_vec[i-wml+1] #score at the first position of the PWM
                            score_array = score_mat[n] #get the source score matrix
                            score_idx = t - wml + 1 #translate t to score_array idx for emission score
                            f_emit_score = score_array[score_idx] #emission score at the last position of the PWM

                            score=logaddexp(score, lps(from_score, f_emit_score, log_motif_expectation))
                        end
                        lh_vec[i] = score
                    end
                    return lh_vec[obsl+1]
                end


text:
@testset "$S" for S in [RefStdString, GCStdString, STL.SharedStdString]
    str = S()
    T = eltype(str)

    @test length(str) == 0
    @test isempty(str)

    if S isa RefStdString
        free(str)
    end

    str = S("Hello, World!")
    @test length(str) == 13

    @test str[1] == T('H')
    @test str[8] == T('W')
    @test str[13] == T('!')

    str[13] = '\0'
    @test length(str) == 13
    @test str[13] == T('\0')

    if S isa RefStdString
        free(str)
    end
    GC.gc(true)
end

# Set reproducible random number seed
Random.seed!(0)
testAbstractChar(StdChar, generator(StdChar))
Random.seed!(0)
testAbstractString(RefStdString, generator(RefStdString))
Random.seed!(0)
testAbstractString(GCStdString, generator(GCStdString))
Random.seed!(0)
testAbstractString(STL.SharedStdString, generator(STL.SharedStdString))


text:
module CodeGeneration

export generatetype, generatestruct, generateimport,
    generateblock, generatemethod, generatemodule

# Generate abstract types

generatetype(name::Symbol) = :(abstract type $name end)
generatetype(name::Symbol, supertype::Symbol) = :(abstract type $name <: $supertype end)

# Generate structs

generatestruct(name::Symbol, fields::Vararg{Tuple{Symbol, Symbol}, N}) where {N} =
    :(struct $name
        $(map(x -> :($(x[1])::$(x[2])), fields)...)
    end)

generatestruct(name::Symbol, supertype::Symbol, fields::Vararg{Tuple{Symbol, Symbol}, N}) where {N} =
    :(struct $name <: $supertype
        $(map(x -> :($(x[1])::$(x[2])), fields)...)
    end)

# Generate import

generateimport(mod::Expr, ::Vararg{Symbol,0}) = :(import $mod)
generateimport(mod::Expr, objs::Vararg{Symbol,N}) where {N} = :(import $mod: $(objs...))

# Generate blocks

generateblock(exprs::Expr...) = Expr(:block, exprs...)

# Generate methods

const SymbolOrExpr = Union{Expr, Symbol}

function generatemethod(
    name::SymbolOrExpr,
    parameters::Vector{T},
    code::Expr,
    ::Vararg{SymbolOrExpr, 0}
) where {T <: SymbolOrExpr}

    :(function $name($(parameters...))
        $code
    end)
end

function generatemethod(
    name::SymbolOrExpr,
    parameters::Vector{T},
    code::Expr,
    whereparams::Vararg{SymbolOrExpr, N}
) where {T <: SymbolOrExpr, N}

    :(function $name($(parameters...)) where {$(whereparams...)}
        $code
    end)
end

# Generate modules

generatemodule(name::Symbol, exprs::Expr...) =
    :(module $name
        $(exprs...)
    end)

end


text:
using InteractiveUtils: subtypes
using .StructArrays: StructArray, components

function _subtypes(type::Type)
  out = Any[]
  _subtypes!(out, type)
end

function _subtypes!(out, type::Type)
  if !isabstracttype(type)
    push!(out, type)
  else
    foreach(T -> _subtypes!(out, T), subtypes(type))
  end
  out
end

# Use lispy tuple programming to iterate over all fields of a `StructArray`.
# Use ugly `@eval` to resolve ambiguities arising from existing `mul!` definitions.
for Dtype in _subtypes(AbstractDerivativeOperator)
  @eval Base.@propagate_inbounds function mul!(dest::A, D::$Dtype, u::A) where {A<:StructArray}
    iterate_mul!(Tuple(components(dest)), D, Tuple(components(u)))
  end

  @eval Base.@propagate_inbounds function mul!(dest::A, D::$Dtype, u::A, α) where {A<:StructArray}
    iterate_mul!(Tuple(components(dest)), D, Tuple(components(u)), α)
  end

  @eval Base.@propagate_inbounds function mul!(dest::A, D::$Dtype, u::A, α, β) where {A<:StructArray}
    iterate_mul!(Tuple(components(dest)), D, Tuple(components(u)), α, β)
  end
end

Base.@propagate_inbounds function iterate_mul!(dest::NTuple{N}, D::AbstractDerivativeOperator, u::NTuple{N}, args...) where N
  mul!(first(dest), D, first(u), args...)
  iterate_mul!(Base.tail(dest), D, Base.tail(u), args...)
end

iterate_mul!(::Tuple{}, D::AbstractDerivativeOperator, ::Tuple{}, args...) = nothing


text:
## Resimulation MH

# MH algorithm that fully resimulates the underlying simulator
function resimulation_mh(T::Int, observations::ChoiceMap, n_iters::Int,
                         tracked_vars=[:β], obs_noise=5.0)
    trs, scores = Vector{Trace}(undef, n_iters) ,Vector{Float64}(undef, n_iters)
    data = DataFrame(fill(Float64, length(tracked_vars)), tracked_vars, n_iters)
    # Generate initial trace
    trace, _ = generate(bayesian_sir, (T, obs_noise), observations)
    trs[1], scores[1] = trace, get_score(trace)
    data[1, tracked_vars] = [trace[v] for v in tracked_vars]
    # Get all addreses for the simulator
    abm_addrs = Gen.select([:step => t => :agents for t in 1:T]...)
    println("Running resimulation MH...")
    @showprogress for i in 2:n_iters
        trace, _ = mh(trace, Gen.select(:β)) # Propose new parameters
        trace, _ = mh(trace, gaussian_drift, ()) # Propose new parameters
        trace, _ = mh(trace, abm_addrs) # Resimulate entire SIR model
        trs[i], scores[i] = trace, get_score(trace)
        data[i, tracked_vars] = [trace[v] for v in tracked_vars]
    end
    return trs, scores, data
end


text:
using Test
using PETSc, MPI, LinearAlgebra, SparseArrays

# This shows some examples of how SNES can be used.
# there are 2 ways to use it:
#    1) you will have julia vectors (but PETSc matrices) within your user-routine
#       That is fine for sequential runs.
#    2) you will receive pointers to PETSc vectors within your user routines.
#       That is important for parallel simulations, where global residual
#       and solution vectors are passed to the user-routines, but we work with
#       local vectors

# Case 1) Julia vectors

x = ones(2)
# this is an example of a case where julia vectors are passed to the routines
function fn_julia!(fx::Vector, x::Vector, args...)
  fx[1] = x[1]^2 + x[1]*x[2] - 3
  fx[2] = x[1]*x[2] + x[2]^2 - 6
end

J = zeros(2,2)
PJ = PETSc.MatSeqDense(J)
function update_jac_julia!(x::Vector, args...)
    J[1,1] = 2x[1] + x[2]
    J[1,2] = x[1]
    J[2,1] = x[2]
    J[2,2] = x[1] + 2x[2]
end

S_julia = PETSc.SNES{Float64}(MPI.COMM_SELF; ksp_rtol=1e-4, pc_type="none", ksp_monitor=false, ksp_converged_reason=true)
PETSc.setfunction!(S_julia, fn_julia!, PETSc.VecSeq(zeros(2)))
PETSc.setjacobian!(S_julia, update_jac_julia!, PJ, PJ)
x_julia = PETSc.VecSeq([2.0, 3.0]);
b_julia = PETSc.VecSeq([0.0, 0.0]);
PETSc.solve!(x_julia, S_julia, b_julia)
@test x_julia.array ≈ [1.0,2.0] rtol=1e-4

# Case 2) PETSc vectors
# this is an example of a case where pointers to PETSc vectors are passed to the routines
function fn!(cfx, cx, user_ctx)
  # We could do Global->Local here on cfx/cx, provided a pointer to the local
  #  vector is available in user_ctx
  x  = PETSc.unsafe_localarray(Float64, cx;  write=false)   # read array
  fx = PETSc.unsafe_localarray(Float64, cfx; read=false)    # write array

  fx[1] = x[1]^2 + x[1]*x[2] - 3
  fx[2] = x[1]*x[2] + x[2]^2 - 6

  Base.finalize(fx)
  Base.finalize(x)
end

function update_jac!(cx, J1::PETSc.AbstractMat{Float64}, args...)
    x  = PETSc.unsafe_localarray(Float64, cx;  write=false)
    @show J1, typeof(J1)
    J1[1,1] = 2x[1] + x[2]
    J1[1,2] = x[1]
    J1[2,1] = x[2]
    J1[2,2] = x[1] + 2x[2]

    Base.finalize(x)
    PETSc.assemble(J1)
end

# structure with which we can pass data to the user-routines above
mutable struct Data
    vec
    julia
end

julia_vec = 0;  # we want pointers to local vectors
S = PETSc.SNES{Float64}(MPI.COMM_SELF, julia_vec;
                ksp_rtol=1e-4,
                pc_type="none",
                ksp_monitor=true,
                ksp_converged_reason=true)

data        = Data([100;2], 1)
S.user_ctx  = data;      # we can pack anything we need in this struct

PETSc.setfunction!(S, fn!, PETSc.VecSeq(zeros(2)))
PJ = PETSc.MatSeqDense(zeros(2,2))
PETSc.setjacobian!(S, update_jac!, PJ, PJ)
x  = PETSc.VecSeq([2.0, 3.0]);
b  = PETSc.VecSeq([0.0, 0.0]);
PETSc.solve!(x, S, b)
@test x.array ≈ [1.0,2.0] rtol=1e-4

# cleanup
#finalize(x)
#finalize(J);
#PETSc.destroy(x0);
#PETSc.destroy(x);
#PETSc.destroy(b);
#PETSc.destroy(PJ);
#PETSc.destroy(S);

#include("test_dmstag.jl")

#PETSc.finalize()


text:
using ReactorApp

push!(ARGS, "arg")
ReactorApp.julia_main()


text:
mutable struct Impact{T} <: Contact{T}
    Nx::Adjoint{T,SVector{3,T}}
    p::SVector{3,T}
    offset::SVector{6,T}

    function Impact(body::Body{T}, normal::AbstractVector; p::AbstractVector = zeros(3),  offset::AbstractVector = zeros(3)) where T
        # Derived from plane equation a*v1 + b*v2 + distance*v3 = p - offset
        V1, V2, V3 = orthogonalcols(normal) # gives two plane vectors and the original normal axis
        A = [V1 V2 V3]
        Ainv = inv(A)
        ainv3 = Ainv[3,SA[1; 2; 3]]
        Nx = ainv3'
        offset = [offset;0.0;0.0;0.0]

        new{T}(Nx, p, offset), body.id
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, constraint::Impact{T}) where {T}
    summary(io, constraint)
    println(io,"")
    println(io, " Nx:     "*string(constraint.Nx))
    println(io, " offset: "*string(constraint.offset))
end


text:
module SpeedDate

export read_fasta, count_mutations, estimate_time, lower, middle, upper

using BioSequences, FASTX, Distributions, Roots, PairwiseListMatrices, NamedArrays

# For bit-parallel mutation counting
import Twiddle: enumerate_nibbles, nibble_mask, count_nonzero_nibbles, count_1111_nibbles

# Needed for using BioSequences' bit-wise operations code generator
import BioSequences: bitindex, encoded_data, index, offset, bitmask

include("read_sequences.jl")
include("count_mutations.jl")
include("estimation_algorithm.jl")
#include("plotting/visualize.jl")

end # Module


text:
module MClassification
    using DataFrames, MLJBase, Distances
    using EasyStream
    export euclidean_distance

    euclidean_distance = Distances.Euclidean()
    include("./MCluster.jl")

    mutable struct MClassifier <: EasyStream.Deterministic
        r_limit::Float64
        n_features::Int64
        metric

        function MClassifier(;r_limit = 0.1, metric = euclidean_distance)
            new_MClassifier = new()
            new_MClassifier.r_limit = r_limit
            new_MClassifier.metric = metric
            return new_MClassifier
        end
    end

    function fit(model::MClassifier, verbosity::Int, X::Array{T, 2}, Y::CategoricalArray) where {T<:Number}
        fitresult = Set{MCluster.MC}()
        for i in 1:length(Y)
            mc = MCluster.MC(X[i, :], Y[i])
            push!(fitresult, mc)
        end
        cache = nothing
        report = nothing
        return fitresult, cache, report
    end

    function fit(model::MClassifier, verbosity::Int, X, Y::CategoricalArray)
        x = Matrix(X)
        fitresult = Set{MCluster.MC}()
        for i in 1:length(Y)
            label = Y[i]
            mc = MCluster.MC(x[i, :], label)
            push!(fitresult, mc)
        end
        cache = nothing
        report = nothing
        return fitresult, cache, report
    end

    function predict(classifier::MClassifier, fitresult, instance::Array{T, 1}) where {T<:Number}
        distances = Array{Any, 1}()
        for micro_cluster in fitresult
            Base.push!(distances, [classifier.metric(instance, micro_cluster.centroid), micro_cluster])
        end

        sort!(distances, by = x -> x[1])
        array_size = length(distances)
        micro_cluster = distances[1][2]
        return micro_cluster.label
    end

    function update_predict(classifier::MClassifier, fitresult, instance::Array{T, 1}) where {T<:Number}
        distances = Array{Any, 1}()
        for micro_cluster in fitresult
            Base.push!(distances, [classifier.metric(instance, micro_cluster.centroid), micro_cluster])
        end

        sort!(distances, by = x -> x[1])
        array_size = length(distances)
        micro_cluster = distances[1][2]

        if MCluster.predict_r(micro_cluster, instance) <= classifier.r_limit
            MCluster.append!(micro_cluster, instance)
        else
            push!(fitresult, MCluster.MC(instance, micro_cluster.label))
            count = 0
            mcs_farthest = Array{MCluster.MC, 1}()

            for i in 0:array_size - 1
                if distances[array_size - i][2].label == distances[1][2].label
                    count += 1
                    push!(mcs_farthest, distances[array_size - i][2])
                    if count == 2
                        merge_farthest(fitresult, mcs_farthest[1], mcs_farthest[2])
                        break
                    end
                end
            end
        end
        return micro_cluster.label
    end

    function update_predict(classifier::MClassifier, fitresult, samples::Array{T, 2}) where {T<:Number}
        return [update_predict(classifier, fitresult, samples[i, :]) for i in 1:nrows(samples)]
    end

    function update_predict(classifier::MClassifier, fitresult, samples) where {T<:Number}
        X = Matrix(samples)
        update_predict(classifier, fitresult, X)
    end

    function append!(classifier::MClassifier, micro_cluster::MCluster.MC)
        push!(classifier.micro_clusters, micro_cluster)
    end

    function merge_farthest(micro_clusters::Set, mc_a::MCluster.MC, mc_b::MCluster.MC)
        mc_a = MCluster.merge(mc_a, mc_b)
        delete!(micro_clusters, mc_b)
    end

    include("./MLJBase.jl")
    include("./EasyStream.jl")
end # module


text:
module GenderDemoInertia
using DataFrames
using Dates
using Distributions
using Statistics
using CSV
using Debugger
using Missings
using DataFramesMeta
using ProgressMeter
import Pandas
import DataValues

include("model_demographic_inertia_gender.jl")
include("../test/conftest.jl")

export run_gender_inertia_simulation, TESTDATA

function run_gender_inertia_simulation(initial_values, runs, duration)
    res = demographic_inertia_simulation(initial_values, runs, duration)
    res;
    end

function get_summaries(initial_values, runs, duration)
    simulation_data = demographic_inertia_simulation(initial_values, runs, duration)
    simulation_data, summary = calculate_simulation_summary(simulation_data)
    return simulation_data, summary ;
end

function get_summaries(initial_values)
    get_summaries(initial_values, initial_values[])
end

#@run run_gender_inertia_simulation(TESTDATA, 2, 20)
#@run run_gender_inertia_simulation(TESTDATA, 100, 6)
#@run GenderDemoInertia.get_summaries(GenderDemoInertia.TESTDATA, 100, 6)

end  # Module end


text:
#
# R2BP systems in our solar system
#

"""
Constant `R2BPParameters` for our sun!
"""
const Sun = R2BPParameters(1.327124400419393e11u"km^3/s^2"; name=:Sun)

"""
Constant `R2BPParameters` for Mercury.
"""
const Mercury = R2BPParameters(22031.78000000002u"km^3/s^2"; name=:Mercury)

"""
Constant `R2BPParameters` for Venus.
"""
const Venus = R2BPParameters(324858.592u"km^3/s^2"; name=:Venus)

"""
Constant `R2BPParameters` for your home planet!
"""
const Earth = R2BPParameters(398600.4354360959u"km^3/s^2"; name=:Earth)

"""
Constant `R2BPParameters` for our moon.
"""
const Moon = R2BPParameters(4902.800066163796u"km^3/s^2"; name=:Moon)

"""
Constant `R2BPParameters` (alias for our mooon).
"""
const Luna = Moon

"""
Constant `R2BPParameters` for Mars.
"""
const Mars = R2BPParameters(42828.37362069909u"km^3/s^2"; name=:Mars)

"""
Constant `R2BPParameters` for Jupiter.
"""
const Jupiter = R2BPParameters(1.2668653492180079e8u"km^3/s^2"; name=:Jupiter)

"""
Constant `R2BPParameters` for Saturn.
"""
const Saturn = R2BPParameters(3.793120749865224e7u"km^3/s^2", name=:Saturn)

"""
Constant `R2BPParameters` for Uranus.
"""
const Uranus = R2BPParameters(5.793951322279009e6u"km^3/s^2"; name=:Uranus)

"""
Constant `R2BPParameters` for Neptune.
"""
const Neptune = R2BPParameters(6.835099502439672e6u"km^3/s^2"; name=:Neptune)

"""
Constant `R2BPParameters` for Pluto. We couldn't leave you out again!
"""
const Pluto = R2BPParameters(869.6138177608749u"km^3/s^2"; name=:Pluto)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing

@testset "test interpolation of discrete constant time-variant field" begin
    f = DCTV(0.0 => 0.0, 1.0 => 1.0)
    # time interpolation of time-variant fields results it's
    # time-invariant counterpart
    @test isapprox(f(-1.0), DCTI(0.0))
    @test isapprox(f( 0.0), DCTI(0.0))
    @test isapprox(f( 0.3), DCTI(0.3))
    @test isapprox(f( 0.5), DCTI(0.5))
    @test isapprox(f( 0.9), DCTI(0.9))
    @test isapprox(f( 1.0), DCTI(1.0))
    @test isapprox(f( 1.5), DCTI(1.0))

    f2 = DCTV(0.0 => 0.0, 0.25 => -0.1, 0.50 => -0.1)
    @test isapprox(f2(0.0), DCTI(0.0))
    @test isapprox(f2(0.25), DCTI(-0.1))
    @test isapprox(f2(0.50), DCTI(-0.1))
    @test isapprox(f2(0.35), DCTI(-0.1))
end


text:
module Highbury

export majority_element, distributed_majority_element, parallel_majority_element

include("majority_voting.jl")

end


text:
module BenchMergewith

using BangBang: mergewith!!, right
using BenchmarkTools: @benchmarkable, BenchmarkGroup

const suite = BenchmarkGroup()

src1 = Dict{Int,Int}(rand(Int, 1000) .=> 1)

for (clabel, cfun) in [:right => right, :+ => +]
    s0 = suite[:combine=>clabel] = BenchmarkGroup()

    s1 = s0[:desttype=>:compatible] = BenchmarkGroup()
    s1[:destkeys=>:empty] =
        @benchmarkable(mergewith!!($cfun, dest, $src1); setup = (dest = Dict{Int,Int}()))
    s1[:destkeys=>:equal] =
        @benchmarkable(mergewith!!($cfun, dest, $src1); setup = (dest = copy($src1)))
    s1[:destkeys=>:half] = @benchmarkable(
        mergewith!!($cfun, dest, $src1);
        setup = (dest = Dict{Int,Int}(collect(keys($src1))[1:end÷2] .=> 1))
    )

    s2 = s0[:desttype=>:widen] = BenchmarkGroup()
    s2[:destkeys=>:empty] =
        @benchmarkable(mergewith!!($cfun, dest, $src1); setup = (dest = Dict{Bool,Int}()))
end

end  # module
BenchMergewith.suite


text:
using Statistics

#=
Grading policy

Your overall score in this class will be a weighted average of your scores for the different components, with the following weights:

16% for the lecture exercises (divided equally among the 16 out of 19 lectures)
18% for the homeworks (divided equally among 6 (out of 7) homeworks)
2% for the Project 0
36% for the Projects (divided equally among 4 (out of 5)
12% for the Midterm exam (timed)
16% for the final exam (timed)

To earn a verified certificate for this course, you will need to obtain an overall score of 60% or more of the maximum possible overall score.

Lecture Exercises, Problem Sets, and Projects

    The lowest 3 scores among the 19 lectures will be dropped, so only 16 out of 19 lectures will count .

    The lowest 1 scores among the 7 homeworks will be dropped, so only 6 out of 7 homeworks will count .

    The lowest 1 score among the 5 projects (excluding project 0) will be dropped, so only 4 out of 5 projects will count .
=#

lectures_grades  = [0/11,0/18,0/12,0/11,0/9,0/13,0/9,0/14,0/5,0/18,0/18,0/11,0/24,0/10,0/24,0/20,0/20,0/9,0/8]
homeworks_grades = [0/61,0/28,0/12,0/25,0/26,0/32,0/36]
projects_grades  = [0/39,0/72,0/39,0/14,0/14] # exclude Project 0 here
project0_grade   = 0
midterm_grade    = 0/38
final_grade      = 0
lectures_grades_retained=sort(lectures_grades)[4:end]
homeworks_grades_retained=sort(homeworks_grades)[2:end]
projects_grades_retained=sort(projects_grades)[2:end]

grades = [mean(lectures_grades_retained),mean(homeworks_grades_retained),project0_grade,mean(projects_grades_retained),midterm_grade,final_grade]

w_lectures  = 0.16
w_homeworks = 0.18
w_project0  = 0.02
w_projects  = 0.36
w_midterm   = 0.12
w_final     = 0.16
weigths     = [w_lectures,w_homeworks,w_project0,w_projects,w_midterm,w_final]

overallGrade = grades' * weigths


text:
#=##############################################################################
# # DESCRIPTION
#     Examples of VTK formatting and usage of VTKtools
# # AUTHORSHIP
#   * Author    : Eduardo J Alvarez
#   * Email     : Edo.AlvarezR@gmail.com
#   * Created   : Nov 2017
#   * License   : MIT License
=###############################################################################

function airfoil_example()
  x,y = readcontour(data_path*"S809.txt"; header_len=2)
  plot_airfoil(x,y; style="--.k")

  println("Close figure and press ENTER")
  readline()

  # Separate upper and lower surfaces to make the contour injective in x
  upper, lower = splitcontour(x,y)

  # Parameterize both surfaces independently
  fun_upper = gt.parameterize(upper[1], upper[2], zeros(Float64, size(upper[1])); inj_var=1)
  fun_lower = gt.parameterize(lower[1], lower[2], zeros(Float64, size(lower[1])); inj_var=1)

  # New discretization for both surfaces
  upper_points = gt.discretize(fun_upper, 0, 1, 50, 1.0)
  lower_points = gt.discretize(fun_lower, 0, 1, 50, 1.0)

  # Put both surfaces back together from TE over the top and from LE over the bottom.
  reverse!(upper_points)                           # Trailing edge over the top
  new_x = [point[1] for point in upper_points]
  new_y = [point[2] for point in upper_points]     # Leading edge over the bottom
  new_x = vcat(new_x, [point[1] for point in lower_points])
  new_y = vcat(new_y, [point[2] for point in lower_points])

  plot_airfoil(new_x, new_y; style="^r", label="Parameterized")
  plot_airfoil(x,y; style="--.k", label="Original")
  legend(loc="best")

  println("Close figure and press ENTER")
  readline()

  # Upper surface sections
  sec1 = (0.35, 20, 3.0, true) # 35% of the line has 20 sections in ratio 3.0 around center
  sec2 = (0.65, 33, 3.0, true) # 65% of the line has 33 sections in ratio 3.0 around center

  # New discretization for both surfaces
  upper_points = gt.multidiscretize(fun_upper, 0, 1, [sec1,sec2])
  lower_points = gt.discretize(fun_lower, 0, 1, 50, 8.0; central=true)

  # Put both surfaces back together from TE over the top and from LE over the bottom.
  reverse!(upper_points)                           # Trailing edge over the top
  new_x = [point[1] for point in upper_points]
  new_y = [point[2] for point in upper_points]     # Leading edge over the bottom
  new_x = vcat(new_x, [point[1] for point in lower_points])
  new_y = vcat(new_y, [point[2] for point in lower_points])

  plot_airfoil(new_x, new_y; style="--^r", label="Parameterized")
end

### INTERNAL FUNCTIONS #########################################################
"Plots the contour of an airfoil given in x,y"
function plot_airfoil(x::Array{Float64,1}, y::Array{Float64,1};
                      label="Airfoil", style="-k", figfactor=1.0,
                      title_str="Airfoil geometry")
  # Sizes the figure
  figsize = [7*1.5,5*0.5]*figfactor
  xmin, xmax = -0.05, 1.05
  yrange = (xmax-xmin)/figsize[1] * figsize[2]
  ymin, ymax = -yrange/2, yrange/2
  fig1 = figure(title_str, figsize=(figsize[1], figsize[2]))
  xlim([xmin, xmax])
  ylim([ymin, ymax])

  PyPlot.plot(x,y, style, label=label)
  xlabel("x")
  ylabel("y")
  grid(true, color="0.8", linestyle="--")
  title(title_str)

  if label!="Airfoil"; legend(loc="best"); end;
end

"Receives a .dat file as pulled from airfoiltools.com containing the x and y
contour coordinates of an airfoil, and returns arrays x and y."
function readcontour(file_name; header_len=1)
  x, y = Float64[], Float64[]

  open(file_name) do f
    for (i,line) in enumerate(eachline(f))

      # Ignores header
      if i<=header_len
        nothing
      # Parses each line
      else
        this_x, this_y = split(line)
        push!(x, parse(Float64, this_x))
        push!(y, parse(Float64, this_y))
      end

    end
  end
  return x,y
end

"""
  Receives an airfoil contour and splits it up in upper and lower surfaces as
  divided by the chord line. It returns `(upper, lower)` with `upper=(x,y)` the
  points of the upper surface, ditto for `lower`. Both `upper` and `lower` are
  given in increasing order in x (i.e., from leading to trailing edge).
"""
function splitcontour(x,y)
  # ERROR CASES
  if !(x[1] in [0.0, 1.0])
    error("Invalid contour. x[1] must be either 0.0 or 1.0, got $(x[1]).")
  end

  # Flag indicating whether the contour start at the trailing or leading edge
  start_TE = x[1]==1.0

  # Find the opposite end of the contour
  end_i = -1
  for (i, xi) in enumerate(x)
    if i==1
      nothing
    # Case of starting from the trailing edge
    elseif start_TE && xi > x[i-1]
      end_i = i-1
      break
    # Case of leading edge
    elseif !start_TE  && xi < x[i-1]
      end_i = i-1
      break
    end
  end

  # ERROR CASE
  if end_i==-1
    error("Logic error! End of airfoil not found!")
  end

  # Splits them up
  x_sec1, y_sec1 = x[1:end_i], y[1:end_i]
  x_sec2, y_sec2 = x[end_i:end], y[end_i:end]

  # Sorts them from LE to TE
  if x_sec1[1] > 0.5; reverse!(x_sec1); reverse!(y_sec1); end;
  if x_sec2[1] > 0.5; reverse!(x_sec2); reverse!(y_sec2); end;

  # Determines upper and lower surfaces
  if Statistics.mean(y_sec1) > Statistics.mean(y_sec2)
    upper = [x_sec1, y_sec1]
    lower = [x_sec2, y_sec2]
  else
    upper = [x_sec2, y_sec2]
    lower = [x_sec1, y_sec1]
  end

  return upper, lower
end
################################################################################


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiMetricsV1beta1PodMetrics <: SwaggerModel
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    timestamp::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: timestamp
    window::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Duration } # spec name: window
    containers::Any # spec type: Union{ Nothing, Vector{IoK8sApiMetricsV1beta1ContainerMetrics} } # spec name: containers

    function IoK8sApiMetricsV1beta1PodMetrics(;metadata=nothing, timestamp=nothing, window=nothing, containers=nothing)
        o = new()
        validate_property(IoK8sApiMetricsV1beta1PodMetrics, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiMetricsV1beta1PodMetrics, Symbol("timestamp"), timestamp)
        setfield!(o, Symbol("timestamp"), timestamp)
        validate_property(IoK8sApiMetricsV1beta1PodMetrics, Symbol("window"), window)
        setfield!(o, Symbol("window"), window)
        validate_property(IoK8sApiMetricsV1beta1PodMetrics, Symbol("containers"), containers)
        setfield!(o, Symbol("containers"), containers)
        o
    end
end # type IoK8sApiMetricsV1beta1PodMetrics

const _property_map_IoK8sApiMetricsV1beta1PodMetrics = Dict{Symbol,Symbol}(Symbol("metadata")=>Symbol("metadata"), Symbol("timestamp")=>Symbol("timestamp"), Symbol("window")=>Symbol("window"), Symbol("containers")=>Symbol("containers"))
const _property_types_IoK8sApiMetricsV1beta1PodMetrics = Dict{Symbol,String}(Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("timestamp")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("window")=>"IoK8sApimachineryPkgApisMetaV1Duration", Symbol("containers")=>"Vector{IoK8sApiMetricsV1beta1ContainerMetrics}")
Base.propertynames(::Type{ IoK8sApiMetricsV1beta1PodMetrics }) = collect(keys(_property_map_IoK8sApiMetricsV1beta1PodMetrics))
Swagger.property_type(::Type{ IoK8sApiMetricsV1beta1PodMetrics }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiMetricsV1beta1PodMetrics[name]))}
Swagger.field_name(::Type{ IoK8sApiMetricsV1beta1PodMetrics }, property_name::Symbol) =  _property_map_IoK8sApiMetricsV1beta1PodMetrics[property_name]

function check_required(o::IoK8sApiMetricsV1beta1PodMetrics)
    true
end

function validate_property(::Type{ IoK8sApiMetricsV1beta1PodMetrics }, name::Symbol, val)
end


text:
using Comonicon

# pkg activate <string> --shared=false

module PkgCmd

using Comonicon

"""
activate the environment at `path`.

# Arguments

- `path`: the path of the environment

# Flags

- `--shared`: whether activate the shared environment
"""
@cast function activate(path; shared::Bool = false)
    println("activating $s (shared=$shared)")
end

@main

end

PkgCmd.command_main()


text:
using Spectra
using Test

#### Dummy x
x = collect(0.1:0.1:100.)

#### gaussienne test
p = [10.,50.,5.]
y = p[1] .* exp.(-log(2) .* ((x.-p[2])./p[3]).^2)

y_calc, y_peaks = gaussiennes([p[1]],[p[2]],[p[3]],x)

@test isapprox(y,y_calc)

#### lorentzienne test
p = [10.,50.,5.]
y = p[1] ./ (1.0 .+ ((x.-p[2])./p[3]).^2)

y_calc, y_peaks = lorentziennes([p[1]],[p[2]],[p[3]],x)

@test isapprox(y,y_calc)

#### pearson7 test
a = [10.,50.,5., 5.]
y = a[1] ./ (1.0 .+ ((x.-a[2])./a[3]).^2 .* (2.0.^(1.0./a[4]) .- 1.0))

y_calc, y_peaks = pearson7([a[1]],[p[2]],[p[3]],[a[4]],x)

@test isapprox(y,y_calc)

#### pseudovoigt test
p = [10.,50.,5.,0.5]

y_lor = p[4] .* (p[1] ./ (1 .+ ((x.-p[2])./p[3]).^2))
y_gauss = (1.0 .-p[4]) .* (p[1] .* exp.(-log(2) .* ((x.-p[2])./p[3]).^2))
y = y_gauss + y_lor

y_calc, y_peaks = pseudovoigts([p[1]],[p[2]],[p[3]],[p[4]],x)
@test isapprox(y,y_calc)

#### normal distribution test
p = [1.,50.,5.]
y = p[1]./(p[3].*sqrt(2.0 .*pi)) .* exp.(-0.5.*((x.-p[2])./p[3]).^2)

y_calc, y_peaks = normal_dist([p[1]],[p[2]],[p[3]],x)
@test isapprox(y,y_calc)


text:
# # Script to wrap PETSc header files with Julia code
# #=
# using Clang_orig.cindex
# using Clang_orig.wrap_c
# using Compat

# import Clang_orig.wrap_c.repr_jl
# =#

# using Clang.cindex
# using Clang.wrap_c
# using Compat

# import Clang.wrap_c.repr_jl

# include("rewriter.jl")
# PETSC_INCLUDE = "../../deps/RealDouble/petsc-3.6.0/include"
# ARCH_INCLUDE = "../../deps/RealDouble/petsc-3.6.0/arch-linux2-c-debug/include"  # TEMPORARY
# petsc_header = [joinpath(PETSC_INCLUDE, "petsc.h")]
# h1 = "/usr/lib/gcc/x86_64-linux-gnu/4.8/include"  # get some C datatype definitions like size_t
# h2 = joinpath(PETSC_INCLUDE, "petscsys.h")
# h3 = joinpath(ARCH_INCLUDE, "petscconf.h")
# # Set up include paths
# clang_includes = ASCIIString[]
# push!(clang_includes, PETSC_INCLUDE)
# push!(clang_includes, ARCH_INCLUDE)
# push!(clang_includes, h1)
# #push!(clang_includes, h2)
# #push!(clang_includes, h3)

# println("clang_includes = ", clang_includes)

# # Clang arguments
# clang_extraargs = [ "-std=c99", "-D", "__STDC_LIMIT_MACROS", "-D", "__STDC_CONSTANT_MACROS", "-D", "PETSC_USE_REAL_SINGLE" ]

# # Callback to test if a header should actually be wrapped (for exclusion)
# function wrap_header(top_hdr::ASCIIString, cursor_header::ASCIIString)
#   return startswith(dirname(cursor_header), PETSC_INCLUDE)
# end

# lib_file(hdr::ASCIIString) = "petsc"
# output_file(hdr::ASCIIString) = "PETSc.jl"

# function wrap_cursor(name::ASCIIString, cursor)
#   println("\nwrapping name = ", name)
#   exc = true
# #= debug: wrap everything
#   if length(name) > 0
#     if name[1] == '_'
#       exc = false
#     end
#   end
# =#
#   return exc
# end

# const wc = wrap_c.init(;
#                         headers = petsc_header,
#                         output_file = "libPETSc_h.jl",
#                         common_file = "libPETSc_common.jl",
#                         clang_includes      = clang_includes,
#                         clang_args          = clang_extraargs,
#                         header_wrapped      = wrap_header,
#                         header_library      = lib_file,
#                         header_outputfile   = output_file,
#                         cursor_wrapped      = wrap_cursor,
#                         clang_diagnostics = true,
#                         rewriter = petsc_rewriter)

# run(wc)


text:
using Test
using BioXP

include("parser.jl")
include("randomize.jl")
include("BioXP.jl")
include("dg.jl")

max_stars_repo_name: hbsmith/BioXP

max_issues_repo_name: hbsmith/BioXP

max_forks_repo_name: hbsmith/BioXP


text:
g = DAG(:A => :M, :A => :D, :M => :D)

@test node(g, :A) == 1
@test_throws LabelNotFound(:U) node(g, :U)

@test nodes_indices(g, [:A, :M]) == [1, 3]
@test_throws LabelNotFound(:X) nodes_indices(g, [:A, :X])

@test nodes_labels(g, [1, 3]) == [:A, :M]

@test edges_indices(g, [(:A => :M), (:M => :D)]) == [1 => 3, 3 => 2]


text:
module Tableaus

    import Polynomials
    import Polynomials: Polynomial

    using ..RungeKutta
    using ..RungeKutta: big, @big, _legendre, _shifted_legendre

    include("tableaus/dirk.jl")

    export TableauCrankNicolson,
           TableauCrouzeix,
           TableauKraaijevangerSpijker,
           TableauQinZhang

    include("tableaus/erk.jl")

    export TableauExplicitEuler, TableauForwardEuler,
           TableauExplicitMidpoint,
           TableauHeun2,
           TableauHeun3,
           TableauRalston2,
           TableauRalston3,
           TableauRunge, TableauRunge2,
           TableauKutta, TableauKutta3,
           TableauRK416, TableauRK4,
           TableauRK438,
           TableauSSPRK2,
           TableauSSPRK3

    include("tableaus/firk.jl")

    export TableauImplicitEuler, TableauBackwardEuler,
           TableauImplicitMidpoint,
           TableauSRK3

    include("tableaus/gauss.jl")

    export TableauGauss

    include("tableaus/lobatto.jl")

    export TableauLobattoIIIA, TableauLobattoIIIB, TableauLobattoIIIC, TableauLobattoIIIC̄,
           TableauLobattoIIID, TableauLobattoIIIE, TableauLobattoIIIF, TableauLobattoIIIG
    export get_lobatto_nullvector

    include("tableaus/radau.jl")

    export TableauRadauIA, TableauRadauIIA,
           TableauRadauIB, TableauRadauIIB

    TableauList = (
       # explicit
       forward_euler         = TableauForwardEuler,
       explicit_euler        = TableauExplicitEuler,
       explicit_midpoint     = TableauExplicitMidpoint,
       heun2                 = TableauHeun2,
       heun3                 = TableauHeun3,
       ralston2              = TableauRalston2,
       ralston3              = TableauRalston3,
       runge                 = TableauRunge,
       runge2                = TableauRunge2,
       kutta                 = TableauKutta,
       kutta3                = TableauKutta3,
       rk4                   = TableauRK4,
       rk416                 = TableauRK416,
       rk438                 = TableauRK438,
       ssprk2                = TableauSSPRK2,
       ssprk3                = TableauSSPRK3,
       # dirk
       crank_nicolson        = TableauCrankNicolson,
       crouzeix              = TableauCrouzeix,
       kraaijevanger_spijker = TableauKraaijevangerSpijker,
       qin_zhang             = TableauQinZhang,
       # firk
       backward_euler        = TableauBackwardEuler,
       implicit_euler        = TableauImplicitEuler,
       implicit_midpoint     = TableauImplicitMidpoint,
       srk3                  = TableauSRK3,
    )

    export TableauList

end


text:
# import Particles package
using Particles

# import HDF5 package for storing solution
using HDF5

# parameters
Δt = 0.1    # time step size
nt = 200    # number of time steps
np = 10000  # number of particles
nx = 16     # number of grid points

# output file
h5file = "vlasov_poisson.hdf5"

# random initial conditions
x₀ = randn(np)
v₀ = randn(np)

# shift x₀ to the interval [0,1]
xmax = ceil(maximum(abs.(x₀)))
x₀ .+= xmax
x₀ ./= 2*xmax

# concatenate x₀ and v₀
z₀ = [ collect(z) for z in zip(x₀, v₀) ]

# vector field
function lorentz_force!(ż, z, p::PoissonSolver)
    for i in eachindex(ż, z)
        ż[i][1] = z[i][2]
        ż[i][2] = eval_field(p, z[i][1])
    end
end

# solution storage
function copy_to_hdf5(h5z, z, n)
    h5z[:,:,n+1] = reduce(hcat, z)
end

# create Poisson solver
p = PoissonSolver{eltype(x₀)}(nx)
solve!(p, x₀)

# create an Equation instance with callbacks
equ = Equation((ż, z) -> lorentz_force!(ż, z, p), z₀; f_preproc = z -> solve!(p, reduce(hcat, z)[1,:]))

# create HDF5 file and copy initial conditions
h5  = h5open(h5file, "w")
h5z = d_create(h5, "z", eltype(x₀), ((2, np, nt+1), (2, np, -1)), "chunk", (2,np,1))
copy_to_hdf5(h5z, z₀, 0)

# create Simulation instance
sim = Simulation(equ, ExplicitEuler, Δt, nt; f_postproc = (z, n) -> copy_to_hdf5(h5z, z, n))

# run simulation
try
    z = sim()
finally
    close(h5)
end

# load Plots package
using Plots

# read array from HDF5 file
z = h5read(h5file, "z")

# compute plot ranges
vmax = ceil(maximum(abs.(v₀)))
xlim = (0, 1)
vlim = (-vmax, +vmax)

# plot initial condition
scatter(mod.(z[1,:,1], 1), z[2,:,1],
        marker = 3,
        xlim = xlim,
        ylim = vlim,
        title = "Vlasov-Poisson",
        legend = false,
        size = (800, 600)
)

# save figure to file
savefig("vlasov_poisson_z₀.png")

# create animation
anim = @animate for n in 0:nt
    scatter(mod.(z[1,:,n+1], 1), z[2,:,n+1],
        marker = 3,
        xlim = xlim,
        ylim = vlim,
        title = "Vlasov-Poisson",
        legend = false,
        size = (800, 600)
    )
end

# save animation to file
gif(anim, "vlasov_poisson_anim.gif", fps=10)

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
# ---
# cover: ../logo.png
# ---

# Configure cover in front matter


text:
__precompile__()

module JuliaSQL

    #includes
    include("io.jl")
    include("structs.jl")

    function connect(con_url::String)
        driver_name, con_info = split(con_url, "://")[1], nothing
        if io.driver_support(driver_name)
            if io.driver_class(driver_name) == "mysql-family"
                #url = "mysql://$uid:$psswd@$host:$port/$schema_name?$args"
                url_splits = split(replace(con_url, "$driver_name://", ""), "/") ; user_info, host_info = split(url_splits[1], "@")
                uid, psswd = split(user_info, ":")
                host, port = split(host_info, ":")
                schema_name, args = split(url_splits[2], "?")[1], split(url_splits[2], "?")[2]
                con_info = Dict(
                    "driver" => io.driver_class(driver_name),
                    "uid" => uid,
                    "psswd" => psswd,
                    "host" => host,
                    "port" => port,
                    "schema" => schema_name,
                    "args" => args != "" ? string("&", args) : args
                )
            elseif io.driver_class(driver_name) == "sqlite-family"
                con_info = Dict(
                    "driver" => io.driver_class(driver_name),
                    "schema" => con_url[length(driver_name) + 4:end],
                    "args" => ""
                )
            end

            result = structs.con(Dict("driver" => con_info["driver"], "url" => io.create_jdbc_url(con_info), "args" => con_info["args"]))
        else
            print("$driver_name is unsupported driver")
            result = false
        end

        return result
    end

end


text:
"""
    CrossEntropy(NeuralNetwork)

Cross entropy loss function.

    Usage:
```julia-repl
CrossEntropy(NeuralNetwork)
```
Arguments:
- `NeuralNetwork` : The NeuralNetwork structure.
"""
function CrossEntropy(NeuralNetwork)
    ΣScore = 0.0

	for i in 1:size(NeuralNetwork.y, 1)
        ΣScore += NeuralNetwork.y[i] * log(1e-15 + NeuralNetwork.output[i])
    end

    MeanΣScore = 1.0 / size(NeuralNetwork.y, 1) * ΣScore
    loss = -MeanΣScore
	return loss
end


text:
"""
$(TYPEDEF)

One sample test for proportion

Constructors
------------
* `OneSampleProp(p1::Real, p0::Real)`

Arguments
---------
* `p1`: True proportion

* `p0`: Proportion under null hypothesis

Fields
------
$(FIELDS)
"""
type OneSampleProp <: TrialTest
    p1::Real
    p0::Real

    # Validator
    function OneSampleProp(p1, p0)

        if !((0 < p1 < 1) & (0 < p0 < 1))

            error("Proportion values must be in (0, 1)")

        end # end if
        new(p1, p0)

    end # end function

end # function

# One sample test for proportion
function hypotheses{T <: OneSampleProp}(test::T, n::Real, std::Void, alpha::Real, side::String)

    diff = test.p1 - test.p0
    se = sqrt(1 / n * test.p1 * (1 - test.p1))
    z = diff / se
    if side == "two"

        p = cdf(ZDIST, z - quantile(ZDIST, 1 - alpha / 2)) + cdf(ZDIST, -z - quantile(ZDIST, 1 - alpha / 2))

    else

        p = cdf(ZDIST, abs(z) - sqrt(test.p0 * (1 - test.p0) / (test.p1 * (1 - test.p1))) * quantile(ZDIST, 1 - alpha))

    end # end if
    return p

end # end function


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

using Markdown
using InteractiveUtils

# ╔═╡ 1663336c-3902-11eb-0158-f54434ea3080
using DrWatson

# ╔═╡ ac5c6510-3901-11eb-0c47-bb3143471434
begin
    @quickactivate

    using Colors
	using PlaymatSimulator

	import PlaymatSimulator.Actors.Image

	US = deserialize("$pd/tmp/user_selection.jls")
	GS = ingredients("$pd/Base/notebooks/game_settings.jl")

	SCREEN_WIDTH = GS.SCREEN_WIDTH
	SCREEN_HEIGHT = GS.SCREEN_HEIGHT
	SCREEN_BORDER = GS.SCREEN_BORDER

	md"""
	## GAME STAGE
	This notebook should define the `stage` of the game. The `stage` should consist of elements to be drawn below the cards and other ui elements such as background images/ANations, some ui elements, game zone/area markers, etc. Provide those objects to the game engine by defining a single dictionary object, `stage = Dict{Symbol, Any}(...)`.
	"""
end

# ╔═╡ fdf59b16-4da2-11eb-204f-0d0e1a522f6c
const SHADE_PATH = "$pd/Base/ui/zones/area_blk.png"

# ╔═╡ a2195ae4-4da2-11eb-2112-e54a53878565
const zone_shade = load(SHADE_PATH)

# ╔═╡ 235edb4c-38fe-11eb-3112-bde291f6f5b5
begin
	dice_faces = [ load("$pd/Base/ui/dice/$fn") for fn in
		readdir("$pd/Base/ui/dice") if occursin("gif", fn) ]

	serialize("$pd/Base/ui/dice/dice_faces.jls", dice_faces)

	dfs = deserialize("$pd/Base/ui/dice/dice_faces.jls")

	function create_die(dfs; id="die_$(randstring(5))", x=0, y=0)
		d = Dice(id,
			[ GIF("dieface_$(randstring(5))", f) for f in dfs ],
			length(dfs),
			)

		for f in d.faces
			f.data[:parent_id] = id
			f.x = x
			f.y = y
		end

		return d
	end

	function create_glass_counter(id="ctr_$(randstring(5))"; x=0, y=0)
		c = Counter("ctr_$(randstring(5))",
			[ GIF("ctrface_$(randstring(5))", load("$pd/Base/ui/counters/glass.gif")) ],
			nothing,
			)
		c.faces[begin].data[:parent_id] = id
		c.faces[begin].x = x
		c.faces[begin].y = y
		return c
	end

	STAGE = OrderedDict(
		:background => Image("$(GS.BKG_NAME)",
			load(GS.BKG_PATH),
			w=SCREEN_WIDTH, h=SCREEN_HEIGHT
			),
		:six_sided_die => create_die(dfs,
			x=ceil(Int32, 0.85SCREEN_WIDTH),
			y=ceil(Int32, 0.9SCREEN_HEIGHT)
			),
		:glass_counter => create_glass_counter(
			x=ceil(Int32, 0.8SCREEN_WIDTH),
			y=ceil(Int32, 0.9SCREEN_HEIGHT)
			),
		"Library" => Image("Library",
			load(SHADE_PATH),
			x=SCREEN_BORDER,
			y=ceil(Int32, SCREEN_HEIGHT * 0.6) + SCREEN_BORDER,
			w=ceil(Int32, SCREEN_WIDTH * 0.15) - SCREEN_BORDER,
			h=ceil(Int32, SCREEN_HEIGHT * 0.4 - 2SCREEN_BORDER),
			alpha=50,
			),
		"Battlefield" => Image("Battlefield",
			load(SHADE_PATH),
			x=ceil(Int32, SCREEN_WIDTH * 0.15) + SCREEN_BORDER,
			y=SCREEN_BORDER,
			w=ceil(Int32, SCREEN_WIDTH * 0.7) - SCREEN_BORDER,
			h=ceil(Int32, SCREEN_HEIGHT - 2SCREEN_BORDER),
			alpha=50,
			),
		"Command" => Image("Command",
			load(SHADE_PATH),
			x=ceil(Int32, SCREEN_WIDTH * 0.85 + SCREEN_BORDER),
			y=SCREEN_BORDER,
			w=ceil(Int32, SCREEN_WIDTH * 0.15 - 2SCREEN_BORDER),
			h=ceil(Int32, SCREEN_HEIGHT * 0.4 - SCREEN_BORDER),
			alpha=50,
			),
		"Graveyard" => Image("Graveyard",
			load(SHADE_PATH),
			x=ceil(Int32, SCREEN_WIDTH * 0.85 + SCREEN_BORDER),
			y=ceil(Int32, SCREEN_HEIGHT * 0.4 + SCREEN_BORDER),
			w=ceil(Int32, SCREEN_WIDTH * 0.15 - 2SCREEN_BORDER),
			h=ceil(Int32, SCREEN_HEIGHT * 0.6 - 2SCREEN_BORDER),
			alpha=50,
			),
		"Hand" => Image("Hand",
			load(SHADE_PATH),
			x=SCREEN_BORDER,
			y=SCREEN_BORDER,
			w=ceil(Int32, SCREEN_WIDTH * 0.15 - SCREEN_BORDER),
			h=ceil(Int32, SCREEN_HEIGHT * 0.6 - SCREEN_BORDER),
			alpha=50,
			),
		)
end

# ╔═╡ Cell order:
# ╟─1663336c-3902-11eb-0158-f54434ea3080
# ╟─ac5c6510-3901-11eb-0c47-bb3143471434
# ╟─fdf59b16-4da2-11eb-204f-0d0e1a522f6c
# ╠═a2195ae4-4da2-11eb-2112-e54a53878565
# ╟─235edb4c-38fe-11eb-3112-bde291f6f5b5


text:
type TarSource <: AbstractSource
  filename::String
  url::String
end
TarSource(url::String) = TarSource(basename(url), dirname(url))

# Downloads a tarfile and untars it
# It expects archive should contain a single directory
function download(package::AbstractPackage, source::TarSource, workdir::String; kwargs...)
  const src = sourcedir(package, workdir)
  const presource = dirname(src)
  const tar = joinpath(presource, source.filename)

  mkpath(presource)
  isfile(tar) || download(source.url * "/" * source.filename, tar)
  cd(presource) do
    if !isdir(src)
      dir = split(readchomp(`tar -tf $tar`), "\n")[1]
      run(`tar -xvf $tar`)
      mv(dir, basename(src))
    end
  end
end

type FileSource <: AbstractSource
  filename::String
  url::String
end

FileSource(url::String) = FileSource(basename(url), dirname(url))

# Downloads a tarfile and untars it
# It expects archive should contain a single directory
function download(package::AbstractPackage, source::FileSource, workdir::String; kwargs...)
  const src = sourcedir(package, workdir)
  const tar = joinpath(src, source.filename)

  mkpath(src)
  isfile(tar) || download(source.url * "/" * source.filename, tar)
end

max_stars_repo_name: UCL/DebbyPacker

max_issues_repo_name: UCL/DebbyPacker

max_forks_repo_name: UCL/DebbyPacker


text:
"""
RationalQuadraticKernel([ρ=1.0[,α=2.0]])
The rational-quadratic kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y)=(1+ρ²||x−y||²/α)^(-α)
```
where `α` is a shape parameter of the Euclidean distance. Check [`GammaRationalQuadraticKernel`](@ref) for a generalization.
"""
struct RationalQuadraticKernel{Tr,Tα<:Real} <: Kernel{Tr}
    transform::Tr
    α::Tα
    function RationalQuadraticKernel{Tr, Tα}(t::Tr, α::Tα) where {Tr, Tα}
        @check_args(RationalQuadraticKernel, α, α > zero(Tα), "α > 1")
        return new{Tr, Tα}(t, α)
    end
end

function RationalQuadraticKernel(ρ::Real=1.0, α::Real=2.0)
    RationalQuadraticKernel(ScaleTransform(ρ),α)
end

function RationalQuadraticKernel(ρ::AbstractVector{<:Real}, α::Real=2.0)
    RationalQuadraticKernel(ARDTransform(ρ), α)
end

function RationalQuadraticKernel(t::Tr, α::Tα=2.0) where {Tr<:Transform, Tα<:Real}
    return RationalQuadraticKernel{Tr, Tα}(t, α)
end

params(k::RationalQuadraticKernel) = (params(transform(k)),k.α)
opt_params(k::RationalQuadraticKernel) = (opt_params(transform(k)),k.α)

@inline kappa(κ::RationalQuadraticKernel, d²::T) where {T<:Real} = (one(T)+d²/κ.α)^(-κ.α)

metric(::RationalQuadraticKernel) = SqEuclidean()

"""
`GammaRationalQuadraticKernel([ρ=1.0[,α=2.0[,γ=2.0]]])`
The Gamma-rational-quadratic kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y)=(1+ρ^(2γ)||x−y||^(2γ)/α)^(-α)
```
where `α` is a shape parameter of the Euclidean distance and `γ` is another shape parameter.
"""
struct GammaRationalQuadraticKernel{Tr, Tα<:Real, Tγ<:Real} <: Kernel{Tr}
    transform::Tr
    α::Tα
    γ::Tγ
    function GammaRationalQuadraticKernel{Tr,Tα,Tγ}(t::Tr, α::Tα, γ::Tγ) where {Tr, Tα<:Real, Tγ<:Real}
        @check_args(GammaRationalQuadraticKernel, α, α > one(Tα), "α > 1")
        @check_args(GammaRationalQuadraticKernel, γ, γ >= one(Tγ), "γ >= 1")
        return new{Tr, Tα, Tγ}(t, α, γ)
    end
end

function GammaRationalQuadraticKernel(ρ::Real=1.0, α::Real=2.0, γ::Real=2.0)
    GammaRationalQuadraticKernel(ScaleTransform(ρ), α, γ)
end

function GammaRationalQuadraticKernel(ρ::AbstractVector{<:Real}, α::Real=2.0, γ::Real=2.0)
    GammaRationalQuadraticKernel(ARDTransform(ρ),α,γ)
end

function GammaRationalQuadraticKernel(t::Tr,α::Tα=2.0,γ::Tγ=2.0) where {Tr<:Transform, Tα<:Real, Tγ<:Real}
    GammaRationalQuadraticKernel{Tr, Tα, Tγ}(t, α, γ)
end

params(k::GammaRationalQuadraticKernel) = (params(k.transform),k.α,k.γ)
opt_params(k::GammaRationalQuadraticKernel) = (opt_params(k.transform),k.α,k.γ)

@inline kappa(κ::GammaRationalQuadraticKernel, d²::T) where {T<:Real} = (one(T)+d²^κ.γ/κ.α)^(-κ.α)

metric(::GammaRationalQuadraticKernel) = SqEuclidean()