text:
#=

  Euler #26 in Julia.

  Problem 26:
  """
  A unit fraction contains 1 in the numerator. The decimal representation of the
  unit fractions with denominators 2 to 10 are given:

      1/2	= 	0.5
      1/3	= 	0.(3)
      1/4	= 	0.25
      1/5	= 	0.2
      1/6	= 	0.1(6)
      1/7	= 	0.(142857)
      1/8	= 	0.125
      1/9	= 	0.(1)
      1/10	= 	0.1

  Where 0.1(6) means 0.166666..., and has a 1-digit recurring cycle. It can be
  seen that 1/7 has a 6-digit recurring cycle.

  Find the value of d < 1000 for which 1/d contains the longest recurring cycle in
  its decimal fraction part.
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

function get_rep_len(n)
    foundRemainders = zeros(n)
    value = 1;
    position = 1;
    while foundRemainders[value] == 0 && value != 0
        foundRemainders[value] = position
        value = round(Int,1+((value*10) % n)) # added 1
        position += 1
    end
    return position-foundRemainders[value]
end

# 0.00195827s
function euler26a()
    maxLen = 0
    maxD = 0
    for d = 2:1000
        len = get_rep_len(d)
        if len > maxLen
            maxLen = len
            maxD = d
        end
    end
    return maxD
end

# Checks only primes
# 0.00048207s
function euler26b()
    maxLen = 0
    maxD = 0
    for d in 2:1000
        if isPrime(d)
            len = get_rep_len(d)
            if len > maxLen
                maxLen = len
                maxD = d
            end
        end
    end
    return maxD
end

# Only primes
# 0.01460937s
function euler26c()
    return sort(primes(999)
        .|>d->[get_rep_len(d),d]
        , by=x->x[1],rev=true)[1][2]
end

# filter primes first
# 0.01558240s
function euler26d()
    return sort(filter(isPrime, 2:1000)
           .|>d->[get_rep_len(d),d]
           , by=x->x[1],rev=true)[1][2]

end

# A variant. But it's still not complete chaining...
# 0.01890910s
function euler26e()
    return (sort((2:1000).|>x->[get_rep_len(x)*x*isPrime(x),x])|>last)[2]
end

# run_euler(euler26a)
run_euler(euler26b)
# run_euler(euler26c)
# run_euler(euler26d)
# run_euler(euler26e)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
for (op, name) in zip((:min, :max), (:_vminimum, :_vmaximum))
  # Reduce over all dimensions
  @eval $name(A, ::Colon) = vreduce($op, A)

  # Reduce over a single dimension
  @eval $name(A, region::Int) = region <= ndims(A) ? vreduce($op, A, dims=region) : A

  # General case: recursive
  @eval function $name(A::AbstractArray, region::Tuple)
    if length(region) == 1
      # For tuple with single element, simply unwrap
      $name(A, region[1])
    elseif length(region) == 2
      # For tuple with two elements, two evaluations suffice
      $name($name(A, region[1]), region[2])
    else
      # Otherwise recurse
      $name($name(A, region[1]), region[2:end])
    end
  end
end

"""
```julia
vminimum(A; dims)
```
Find the least value contained in `A`, optionally over dimensions specified by `dims`.
As `Base.minimum`, but vectorized

## Examples
```julia
julia> using VectorizedStatistics

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> vminimum(A, dims=1)
1×2 Matrix{Int64}:
 1  2

julia> vminimum(A, dims=2)
 2×1 Matrix{Int64}:
 1
 3
```
"""
vminimum(A; dims=:) = _vminimum(A, dims)
export vminimum

"""
```julia
vmaximum(A; dims)
```
Find the greatest value contained in `A`, optionally over dimensions specified by `dims`.
As `Base.maximum`, but vectorized.

## Examples
```julia
julia> using VectorizedStatistics

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> vmaximum(A, dims=1)
1×2 Matrix{Int64}:
 3  4

julia>  vmaximum(A, dims=2)
 2×1 Matrix{Int64}:
 2
 4
```
"""
vmaximum(A; dims=:) = _vmaximum(A, dims)
export vmaximum

"""
```julia
vextrema(A; dims)
```
Find the maximum and minimum of `A`, optionally along the dimensions specified by `dims`.
As `Base.extrema`, but vectorized.

## Examples

julia> A = reshape(Vector(1:2:16), (2,2,2))
2×2×2 Array{Int64, 3}:
 [:, :, 1] =
  1  5
  3  7

 [:, :, 2] =
   9  13
  11  15

julia> extrema(A, dims = (1,2))
1×1×2 Array{Tuple{Int64, Int64}, 3}:
 [:, :, 1] =
  (1, 7)

 [:, :, 2] =
  (9, 15)
"""
vextrema(A; dims=:) = _vextrema(A, dims)
_vextrema(A, region) = collect(zip(_vminimum(A, region), _vmaximum(A, region)))
_vextrema(A, ::Colon) = (_vminimum(A, :), _vmaximum(A, :))
export vextrema

ext: jl


text:
function get_filepath_info_no_ext(filepath_no_ext)
    pattern = r"^(?P<dir>.+\\)*(?P<file>.+)$"
    _matched = match(pattern, filepath_no_ext)
    return (typeof(_matched["dir"]) == Void ? "" : _matched["dir"]), _matched["file"]
end

function write_simple_file(vtkobject::T, filepath_or_vtmfile::Union{String, WriteVTK.MultiblockFile}) where T<:AbstractVTKSimpleData
    vtkfile = nothing
    if T <: VTKUnstructuredData
        points = vtkobject.point_coords
        cells = [MeshCell(VTK_CELL_TYPE[vtkobject.cell_types[i]], vtkobject.cell_connectivity[i])
            for i in 1:num_of_cells(vtkobject)]
        vtkfile = vtk_grid(filepath_or_vtmfile, points, cells)
    elseif T <: VTKStructuredData
        points = vtkobject.point_coords
        vtkfile = vtk_grid(filepath_or_vtmfile, points)
    elseif T <: VTKRectilinearData
        points = vtkobject.point_coords
        vtkfile = vtk_grid(filepath_or_vtmfile, points...)
    elseif T <: VTKUniformRectilinearData
        vtkfile = vtk_grid(filepath_or_vtmfile, vtkobject.extents..., origin=vtkobject.origin, spacing=vtkobject.spacing)
    end

    for v in keys(vtkobject.point_data)
        vtk_point_data(vtkfile, vtkobject.point_data[v], v)
    end
    for v in keys(vtkobject.cell_data)
        vtk_cell_data(vtkfile, vtkobject.cell_data[v], v)
    end

    return vtk_save(vtkfile)[1]
end

function write_blocked_file(multiblock_object::AbstractVTKMultiblockData, filename_noext::String)
    xvtm = XMLDocument()
    xroot = create_root(xvtm, "VTKFile")
    set_attribute(xroot, "type", "vtkMultiBlockDataSet")
    set_attribute(xroot, "version", "1.0")
    if ENDIAN_BOM == 0x04030201
        set_attribute(xroot, "byte_order", "LittleEndian")
    else
        set_attribute(xroot, "byte_order", "BigEndian")
    end
    new_child(xroot, "vtkMultiBlockDataSet")
    xMBDS = find_element(xroot, "vtkMultiBlockDataSet")

    filenames = []
    i = 1
    for block in simple_block_generator(multiblock_object)
        vtkFilename_noext = @sprintf("%s.z%02d", filename_noext, i)
        if isa(block, VTKPolyData)
            push!(filenames, splitdir(write_vtp(block, vtkFilename_noext))[2])
        else
            push!(filenames, splitdir(write_simple_file(block, vtkFilename_noext))[2])
        end
        i += 1
    end

    for i in 1:length(filenames)
        fname = filenames[i]
        xBlock = new_child(xMBDS, "Block")
        set_attribute(xBlock, "index", "$i")

        xDataSet = new_child(xBlock, "DataSet")
        set_attribute(xDataSet, "index", "0")
        set_attribute(xDataSet, "file",  fname)
    end
    vtmfilepath = string(filename_noext, ".vtm")
    save_file(xvtm, vtmfilepath)

    return vtmfilepath
end

function write_timeseries_file(timeseries_object::AbstractTimeSeriesVTKData, filename_noext::String)
    xpvd = XMLDocument()
    xroot = create_root(xpvd, "VTKFile")
    set_attribute(xroot, "type", "Collection")
    set_attribute(xroot, "version", "1.0")
    if ENDIAN_BOM == 0x04030201
        set_attribute(xroot, "byte_order", "LittleEndian")
    else
        set_attribute(xroot, "byte_order", "BigEndian")
    end
    set_attribute(xroot, "compressor", "vtkZLibDataCompressor")
    new_child(xroot, "Collection")
    xMBDS = find_element(xroot, "Collection")

    main_dir, file = get_filepath_info_no_ext(filename_noext)
    subdir = main_dir*file
    mkpath(subdir)

    _length = length(timeseries_object)
    no_of_digits = 1
    while no_of_digits*10 - 1 <= _length
        no_of_digits += 1
    end

    filenames = []
    for i in 1:length(timeseries_object)
        new_filepath_no_ext = string(subdir, "\\", file, "_", "0"^(no_of_digits - length(string(i))), i)
        T = typeof(timeseries_object[i])
        if T <: VTKPolyData
            push!(filenames, write_vtp(timeseries_object[i], new_filepath_no_ext))
        elseif T <: AbstractVTKSimpleData
            push!(filenames, write_simple_file(timeseries_object[i], new_filepath_no_ext))
        elseif T <: AbstractVTKMultiblockData
            push!(filenames, write_blocked_file(timeseries_object[i], new_filepath_no_ext))
        end
    end

    for i in 1:length(filenames)
        fname = filenames[i]
        t = timeseries_object.timemarkers[i]
        xDataSet = new_child(xMBDS, "DataSet")
        set_attribute(xDataSet, "timestep", @sprintf("%f", t))
        set_attribute(xDataSet, "part", "0")
        set_attribute(xDataSet, "file", fname)
    end

    pvdfilepath = string(filename_noext, ".pvd")
    save_file(xpvd, pvdfilepath)

    return pvdfilepath
end

function valid_to_write(vtkobject::AbstractStaticVTKData)
    T = typeof(vtkobject)
    T <: AbstractVTKMultiblockData && begin
        _out = true
        for i in vtkobject
            _T = typeof(i)
           (_out = _T.parameters[1] <: Union{Float64,Int}) || break
        end
        _out
    end || T.parameters[1] <: Union{Float64,Int} || return false, "Only float or integer data can be written."
    return true, ""
end

function valid_to_write(vtkobject::AbstractTimeSeriesVTKData)
    for b in vtkobject
       valid, _error = valid_to_write(b)
       valid ? continue : return false, _error
    end
    return true, ""
end

function write_vtk(vtkobject::T, filepath_no_ext::String; time_resolution::Int=1, validation=false, repeat_cells=false) where T<:AbstractVTKData
    if validation
        valid, _error = is_valid(vtkobject, repeat_cells=repeat_cells)
        valid || throw("Invalid data, cannot be written: $_error")
        valid, _error = valid_to_write(vtkobject)
        valid || throw("Invalid data, cannot be written: $_error")
    end

    if T <: VTKPolyData
        return write_vtp(vtkobject, filepath_no_ext)
    elseif T <: AbstractVTKSimpleData
        return write_simple_file(vtkobject, filepath_no_ext)
    elseif T <: AbstractVTKMultiblockData
        return write_blocked_file(vtkobject, filepath_no_ext)
    elseif T <: AbstractTimeSeriesVTKData
        if time_resolution == 1
            return write_timeseries_file(vtkobject, filepath_no_ext)
        else
            _vtkobject = deepcopy(vtkobject)
            increase_resolution!(_vtkobject, time_resolution)
            return write_timeseries_file(_vtkobject, filepath_no_ext)
        end
    end
    throw("Type of object input is not supported.")
end


text:
module AffineTest

using Projections, Test

@testset "affine" begin
A = reshape(collect(LinRange(0, 1, 25)), 5, 5)[:, 1:4]
b = collect(LinRange(1, 5, 5))
y = ones(size(A, 2))
@test isapprox(
    project(Affine(A, b), y),
    [15.36,  9.12,  2.88,  -3.36],
    rtol = 1e-4)
end
end


text:
import Clustering
import PyPlot as Plt
import Seaborn as Sns

function clustermap(mtx::T where T <: AbstractMatrix{Float64}, gene_names::Vector{String}; gene_ord::Union{<:AbstractVector{<:Integer}, Bool}=true,
        cell_ord::Union{<:AbstractVector{<:Integer},Bool}=true, diag_genes::Bool=false, transpose::Bool=false, kwargs...)
    if typeof(cell_ord) === Bool
        if cell_ord
            cell_dists = 1 .- cor(mtx);
            cell_ord = Clustering.hclust(cell_dists, linkage=:ward).order;
        else
            cell_ord = 1:size(mtx, 2)
        end
    end

    if typeof(gene_ord) === Bool
        if gene_ord
            if diag_genes
                gene_ord = sortperm(vec(mapslices(x -> findmax(x)[2], mtx[:, cell_ord], dims=2)), rev=true);
            else
                gene_dists = 1 .- cor(mtx');
                gene_ord = Clustering.hclust(gene_dists, linkage=:ward).order;
            end
        else
            gene_ord = 1:size(mtx, 1)
        end
    end

    mtx = mtx[gene_ord, cell_ord]
    tick_fun = Plt.yticks
    if !transpose
        mtx = copy(mtx')
        tick_fun = Plt.xticks
    end

    plt = Sns.heatmap(mtx; kwargs...)
    tick_fun(1:length(gene_ord), gene_names[gene_ord], rotation=(transpose ? 0 : 90))
    return plt, cell_ord, gene_ord
end

function label_axis!(ax, label::String; x=-0.2, y=1.05, fontsize=24)
    ax.text(x, y, label, transform=ax.transAxes, fontsize=fontsize, va="top", ha="right", fontfamily="sans-serif", fontname="Helvetica")
end

function plot_embedding(embedding::Matrix{<:Real}, groups::Vector; mark_groups::Bool=true, ax=Plt.gca(),
        rasterized::Bool=true, fontsize::Int=10, palette::Union{Nothing, Vector{String}}=nothing, kwargs...)
    annot_labels = sort(unique(groups))
    ids_per_clust = B.split_ids(denserank(groups))
    kwargs = Dict{Symbol, Any}(kwargs...)
    for (i,ids) in enumerate(ids_per_clust)
        if palette !== nothing
            kwargs[:color] = palette[i]
        end

        ax.scatter(embedding[1,ids], embedding[2,ids]; label=annot_labels[i], rasterized=rasterized, kwargs...)
    end

    if mark_groups
        clust_centers = hcat([median(embedding[:, ids], dims=2) for ids in ids_per_clust]...);
        bbox_style = Dict(:facecolor=>"white", :alpha=>0.75, :edgecolor=>"none", :pad=>0.0)
        for (t, x, y) in zip(annot_labels, clust_centers[1,:], clust_centers[2,:])
            ax.annotate(t, (x,y), horizontalalignment="center", fontsize=fontsize, bbox=bbox_style)
        end
    end

    return ax
end

max_stars_repo_name: hiraksarkar/BaysorAnalysis

max_issues_repo_name: hiraksarkar/BaysorAnalysis

max_forks_repo_name: hiraksarkar/BaysorAnalysis


text:
using Transformers.Basic
using PreSumm.Search

struct BeamSearch
    width::Integer
    steps::Integer
    length_normalization::AbstractFloat
end

function (search::BeamSearch)(
    vocabulary::Vocabulary{T},
    predict::Function;
    start_token::T="[CLS]",
    end_token::T="[SEP]"
)::AbstractVector{T} where T
    return beam_search(
        search.width,
        vocabulary.list,
        predict,
        sequence -> last(sequence) != end_token,
        search.steps,
        initial_sequence=[start_token],
        length_normalization=search.length_normalization
    )
end

function BeamSearch(
    width::Integer,
    steps::Integer;
    length_normalization::AbstractFloat=0.0,
)
    return BeamSearch(width, steps, length_normalization)
end


text:
# DESCRIPTION
# Closed-loop aircraft dynamics.
#
# INPUT
# - dx: state derivative
# - x: current state vector
# - u_curr: current control vector
# - t: time
#
# OUTPUT
# - xdot: system dynamics
#
# ------------------------------------------------------------------------------

function GenHel_DI!(dx,x,u_curr,t)
    # unpack states
    u = x[1]
    v = x[2]
    w = x[3]
    p = x[4]
    q = x[5]
    r = x[6]
    ϕ = x[7]
    # unpack command filter states
    p_cmd = x[NSTATES+1]
    intperr_cmd = x[NSTATES+2]
    q_cmd = x[NSTATES+3]
    intqerr_cmd = x[NSTATES+4]
    r_cmd = x[NSTATES+5]
    intrerr_cmd = x[NSTATES+6]
    δ_lat = x[NSTATES+7]
    δ_lon = x[NSTATES+8]
    δ_ped = x[NSTATES+9]
    # control perturbations
    Δδ_lat = u_curr[1]-δ_lat
    Δδ_lon = u_curr[2]-δ_lon
    Δδ_ped = u_curr[4]-δ_ped
    # feedforward
    pdot_cmd = -1/τ_p*p_cmd + k_lat/τ_p*Δδ_lat
    qdot_cmd = -1/τ_q*q_cmd + k_lon/τ_q*Δδ_lon
    rdot_cmd = -1/τ_r*r_cmd + k_ped/τ_r*Δδ_ped
    # feedback compensation
    νp = pdot_cmd + kp_roll*(p_cmd-p) + ki_roll*(intperr_cmd)
    νq = qdot_cmd + kp_pitch*(q_cmd-q) + ki_pitch*(intqerr_cmd)
    νr = rdot_cmd + kp_yaw*(r_cmd-r) + ki_yaw*(intrerr_cmd)
    # pseudo control
    ν=[νp, νq, νr]
    # total velocity
    V=sqrt(u^2+v^2+w^2)/1.688
    # matrix interpolation
    idx1=Int(min(max(floor(V/20)+1,1),7))
    idx2=idx1+1;
    intscale=min((V-(idx1-1)*20)/20,1.)
    CA=CATAB[:,:,idx1]+intscale*(CATAB[:,:,idx2]-CATAB[:,:,idx1])
    CBinv=CBinvTAB[:,:,idx1]+intscale*(CBinvTAB[:,:,idx2]-CBinvTAB[:,:,idx1])
    # dynamic inversion
    Δu=CBinv*(ν-CA*[p, q, r])
    # control inputs
    u=[Δu[1]+δ_lat; Δu[2]+δ_lon; u_curr[3]; Δu[3]+δ_ped]
    # state derivative
    der = GenHel!(dx[1:NSTATES],x[1:NSTATES],u,t)
    # rigid-body state derivative
    xdot=zeros(NSTATES+NDISTATES,1)
    xdot[1:NSTATES]=der
    # command filters state derivative
    xdot[NSTATES+1] = pdot_cmd
    xdot[NSTATES+2] = p_cmd-p
    # turn compensation
    if V<40
        xdot[NSTATES+3] = qdot_cmd
    elseif V<60
        xdot[NSTATES+3] = qdot_cmd+r*sin(ϕ)*((V-40)/(60-40))
    else
        xdot[NSTATES+3] = qdot_cmd+r*sin(ϕ)
    end
    xdot[NSTATES+4] = q_cmd-q
    # turn coordination
    if V<40
        xdot[NSTATES+5] = rdot_cmd
    elseif V<60
        xdot[NSTATES+5] = rdot_cmd+g/(V*1.688)*ϕ*((V-40)/(60-40))
    else
        xdot[NSTATES+5] = rdot_cmd+g/(V*1.688)*ϕ
    end
    xdot[NSTATES+6] = r_cmd-r
    xdot[NSTATES+7] = (u_curr[1]-δ_lat)/τ_u
    xdot[NSTATES+8] = (u_curr[2]-δ_lon)/τ_u
    xdot[NSTATES+9] = (u_curr[4]-δ_ped)/τ_u
    #
    return xdot
end

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
# Continuous interpolation of discrete frequency-domain data
using Interpolations
import Interpolations: interpolate

function interpolate(s::Spectrum)
    l,h = limits(s)
    fn = extrapolate(interpolate(s[:], BSpline(Cubic(Natural())), OnGrid()), Flat())
    # Interpolate.jl only supports increasing ranges, so will have to
    # use flip `h` and `l` here and invert the axis in the next line.
    scaled = scale(fn, range(l; stop=h, length=length(s)))
    δ -> scaled[l + (h - δ)]
end

function resample(s::Spectrum, shifts::AbstractArray)
    intp = interpolate(s)
    intp.(shifts)
end


text:
/* allow comparing booleans. */

int main() {
  if (true == true) { printInt(42); }
  return 0 ;

}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
import Munkres: munkres

function GPM_orthogonal{T}(C::Matrix{T}, x₁::Matrix{T}, num_iters::Int, r::Int = 3, c::Int = 3)
    α = max(0, -eigmin(C))
    C̃ = C + α*I

    x = Vector{Matrix{T}}(num_iters)
    x[1] = x₁

    for i in 2:num_iters
        D = C̃ * x[i-1]
        O = split_blocks(D, r, c)
        map!(nearest_orthogonal, O)
        x[i] = join_blocks(O)
    end

    return x
end

function GPM_permutation{T}(C::Matrix{T}, x₁::Matrix{T}, num_iters::Int, d::Int)
    α = max(0, -eigmin(C))
    C̃ = C + α*I

    x = Vector{Matrix{T}}(num_iters)
    x[1] = x₁

    for i in 2:num_iters
        D = C̃ * x[i-1]
        O = split_blocks(D, d, d)
        map!(O) do Q
            eye(Int, d)[:, munkres(-Q')]
        end
        x[i] = join_blocks(O)
    end

    return x
end


text:
module GCMMesh

include("Mesh.jl")
include("Graphs.jl")
include("Partition.jl")

end # module

max_stars_repo_name: CliMA/GCMMesh

max_issues_repo_name: CliMA/GCMMesh

max_forks_repo_name: CliMA/GCMMesh


text:
#Julia Vega Plotting

#Pkg.add("Vega")
using Vega
stock = ["chairs", "tables", "desks", "rugs", "lamps"];
quantity = [15, 10, 10, 5, 20];
piechart(x = stock, y = quantity)


text:
using Plots, Test
@testset "lens!" begin
    pl = plot(1:5)
    lens!(pl, [1,2], [1,2], inset = (1, bbox(0.0,0.0,0.2,0.2)), colorbar = false)
    @test length(pl.series_list) == 4
    @test pl[2][:colorbar] == :none
end # testset

@testset "vline, vspan" begin
    vl = vline([1], widen = false)
    @test Plots.xlims(vl) == (1,2)
    @test Plots.ylims(vl) == (1,2)
    vl = vline([1], xlims=(0,2), widen = false)
    @test Plots.xlims(vl) == (0,2)
    vl = vline([1], ylims=(-3,5), widen = false)
    @test Plots.ylims(vl) == (-3,5)

    vsp = vspan([1,3], widen = false)
    @test Plots.xlims(vsp) == (1,3)
    @test Plots.ylims(vsp) == (0,1) # TODO: might be problematic on log-scales
    vsp = vspan([1,3], xlims=(-2,5), widen = false)
    @test Plots.xlims(vsp) == (-2,5)
    vsp = vspan([1,3], ylims=(-2,5), widen = false)
    @test Plots.ylims(vsp) == (-2,5)
end # testset


text:
"""
2-D/ 3-D elastic forward modeling using staggered-grid velocity-stress finite-difference formulation.
"""
struct FdtdElastic end
"""
2-D/ 3-D acoustic forward modeling using staggered-grid velocity-stress finite-difference formulation.
"""
struct FdtdAcou end
"""
2-D/ 3-D viscoacoustic forward modeling using staggered-grid velocity-stress finite-difference formulation.
"""
struct FdtdAcouVisco end
"""
2-D Linearized forward modeling using a finite-difference simulation of the acoustic wave-equation.
"""
struct FdtdAcouBorn end

# define supersource, source and receiver structs
struct Srcs
    n::Int
end
Srcs() = Srcs(0)
struct SSrcs
    n::Int
end
SSrcs() = SSrcs(0)
struct Recs
    n::Int
end
Recs() = Recs(0)


text:
export read_quake

@doc """
    Ev = read_quake(fmt, file [, keywords])

Read data in file format `fmt` from `file` into SeisEvent object `Ev`.

| Format      | String          | Notes                                       |

|:---         |:---             |:----                                        |
| PC-SUDS     | suds            |                                             |
| QuakeML     | qml, quakeml    | only reads first event from file            |
| UW          | uw              |                                             |

|KW       | Used By   | Type    | Default   | Meaning                         |
|:---     |:---       |:---     |:---       |:---                             |
| full    | suds, uw  | Bool    | false     | read full header into `:misc`?  |
| v       | all       | Int64   | 0         | verbosity                       |

### Notes on Functionality
* No "in-place" version of `read_quake` exists because earthquake data are
usually discrete, self-contained files.
* `read_quake` doesn't use file wildcards. See `?UW.readuwevt` for help with UW
file string syntax.

See also: read_data, get_data, read_meta
""" read_quake
function read_quake(fmt::String, fname::String;
  full    ::Bool    = false,              # full header
  v       ::Int64   = KW.v                # verbosity level
  )

  Ev = (
    if fmt == "suds"
      SUDS.readsudsevt(fname, full=full, v=v)
    elseif fmt == "uw"
      UW.readuwevt(fname, full=full, v=v)
    elseif fmt in ("qml", "quakeml")
      hdr, source = read_qml(fname)
      SeisEvent(hdr = hdr[1], source = source[1])
    end
    )
  return Ev
end


text:
## ExoplanetsSysSim/src/orbit.jl
## (c) 2015 Eric B. Ford

#export Orbit

immutable Orbit
  P::Float64             # days             # QUERY:  Should we store P or a here?
  ecc::Float64
  incl::Float64          # radians
  omega::Float64         # radians
  asc_node::Float64      # radians
  mean_anom::Float64     # radians          # QUERY:  Should we store t_0 or mean_anom here?
end

#Orbit() = Orbit(0.0,0.0,0.0,0.0,0.0,0.0)   # Comment out, so don't accidentally have invalid orbits

# This will only work if Orbit were mutable.  Is that better or worse?  Let's test and see....
function set!(o::Orbit, P::Float64, e::Float64, i::Float64, w::Float64, asc_node::Float64, M::Float64)
 o.P = P
 o.ecc = e
 o.incl = i
 o.omega = w
 o.asc_node = asc_node
 o.mean_anom = M
 return o
end

function test_orbit_constructors()
  #orb = Orbit()
  orb = Orbit(1.0, 0.03, 0.5*pi,0.0,0.0,pi)
  if !isimmutable(orb)
     set!(orb,1.0, 0.03, 0.5*pi,0.0,0.0,pi)
  end
end


text:
"""
$(SIGNATURES)

Optimizer for the entropy closure: `argmin(<η(α*m)> - α*u)`
"""
function optimize_closure(α, m, ω, u, η::Function; optimizer = Newton())
    loss(_α) = sum(η.(_α' * m)[:] .* ω) - dot(_α, u)
    res = Optim.optimize(loss, α, optimizer)

    return res
end

"""
$(SIGNATURES)

Reconstruct realizable moments based on Lagrange multiplier α
"""
function realizable_reconstruct(α, m, ω, η_dual_prime::Function)
    u = zero(α)
    for i in eachindex(u)
        u[i] = sum(η_dual_prime.(α' * m)[:] .* ω .* m[i, :])
    end

    return u
end


text:
using ProjAssistant
@quickactivate "Chemostat_Rath2017"

@time begin

    using Plots
    using Statistics
    using ProgressMeter
    using Base.Threads

    import Chemostat
    import Chemostat.MetNets
    import Chemostat.MetLP
    import Chemostat.MetEP
    const Ch = Chemostat

    import Chemostat_Rath2017
    const ChR = Chemostat_Rath2017
    const Rd = ChR.RathData
    const H1 = ChR.Human1
    const HG = H1.HumanGEM
    const AG = H1.HartGEMs

end

## ---------------------------------------------------------------------
# DESCRIPTION
# The intake of glucose must be the limiting external flux of the network

## ---------------------------------------------------------------------
let
    tissue = "GBM"
    modelid = "fva_base"
    model = AG.load_model(;modelid, tissue, uncompress = true)
    objidx = MetNets.rxnindex(model, HG.HUMAN_BIOMASS_IDER)

    exchis = filter(MetNets.exchanges(model)) do exchi
        MetNets.lb(model, exchi) < 0
    end
    @show length(exchis)

    # Open all but glc
    for exchi in exchis
        MetNets.lb!(model, exchi, -1000)
    end
    MetNets.lb!(model, HG.HUMAN_GLC_EX_IDER, -0.5)

    results = Dict()
    @showprogress for exchi in exchis

        exch = model.rxns[exchi]
        dat = get!(results, exch, Float64[])

        MetNets.summary(model, model.rxns[exchi])
        println()

        lb_ = MetNets.lb(model, exchi)

        for f in 0.0:0.1:1.0
            MetNets.lb!(model, exchi, lb_ * f)
            fbaout = MetLP.fba!(model, objidx)
            biom = MetLP.objval(fbaout)
            push!(dat, biom)
        end

        MetNets.lb!(model, exchi, lb_)
    end

    sdat(AG, results,
        "glc_limited_study", ".jls";
        verbose = true
    )

end;

## ---------------------------------------------------------------------
let
    met_readable = HG.load_met_readable_ids()
    met_map = HG.load_exch_met_map()
    results = ldat(AG, "glc_limited_study", ".jls")

    li = 5
    bioms = Float64[]
    mets = String[]
    for (exch, dat) in results
        met = get(met_map, exch, exch)
        met = get(met_readable, met, met)
        push!(bioms, dat[li])
        push!(mets, met)
    end

    rang_ = 1:10
    sis = sortperm(bioms)
    bar(mets[sis][rang_], bioms[sis][rang_]; xrotation = 45)

    # sfig(AG, ps,
    #     @fileid, "glc_limited_study", ".png";
    #     layout = (3,3)
    # )
end


text:
module TestLocks

using ArgCheck: @check
using Julio
using Test
using ..Utils: trywait

function test_trivial()
    lck = Julio.Lock()
    lock(lck)
    unlock(lck)
    @test true
end

function test_reentrancy()
    lck = Julio.Lock()
    Julio.withtaskgroup() do tg
        lock(lck)
        t = Julio.spawn!(tg) do
            lock(lck)
            unlock(lck)
        end
        @test !trywait(t)
        lock(lck)
        lock(lck)
        unlock(lck)
        unlock(lck)
        @test !trywait(t)
        unlock(lck)
        @test trywait(t, 3)
    end
    @test true
end

function test_exclusivity()
    lck = Julio.Lock()
    i1, o1 = Julio.channel()
    i2, o2 = Julio.channel()
    begin
        begin
            Julio.withtaskgroup() do tg
                Julio.spawn!(tg) do
                    lock(lck) do
                        put!(i1, :acquired)
                        put!(i1, :releasing)
                        put!(i2, :t1_releasing)
                    end
                end
                @test take!(o1) === :acquired
                @check !trylock(lck)
                Julio.spawn!(tg) do
                    put!(i2, :t2_locking)
                    lock(lck) do
                        put!(i2, :t2_acquired)
                    end
                end
                @test take!(o2) === :t2_locking
                @test take!(o1) === :releasing
                @test take!(o2) === :t1_releasing
                @test take!(o2) === :t2_acquired
            end
        end
    end
end

function test_repeat(; ntasks = Threads.nthreads() * 4, nrepeat = 1000)
    lck = Julio.Lock()
    Julio.withtaskgroup() do tg
        for _ in 1:ntasks
            Julio.spawn!(tg) do
                for i in 1:nrepeat
                    lock(lck) do
                    end
                    # i % 16 == 0 && Julio.yield()
                end
            end
        end
    end
    @test true
end

function test_cancel()
    scope = Julio.cancelscope()
    sept, rept = Julio.channel()
    lck = Julio.Lock()
    Julio.withtaskgroup() do tg
        Julio.spawn!(tg) do
            open(scope) do
                lock(lck) do
                    put!(sept, 111)
                    put!(sept, 222)
                end
            end
        end
        @test take!(rept) == 111
        Julio.spawn!(tg) do
            lock(lck) do
                put!(sept, 333)
            end
        end
        Julio.cancel!(scope)
        @test take!(rept) == 333
    end
end

end  # module


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

module EyeTrackingUtils

export make_clean_data

include(joinpath(@__DIR__, "prep.jl"))

end # end module


text:
module TestIntegration

using TimeZoneLookup
using TimeZoneLookup: V, query
using Setfield
using Combinatorics
using Test

function calcpoints(tree, maxn = 3)
    res = ntuple(i -> 0, maxn)
    for i in 0:50
        for j in 0:50
            face = query(tree, V(i, j)).face
            face == -1 && continue
            @set! res[face] += 1
        end
    end

    return res
end

@testset "Integrations" begin
@testset "Simple parallelogram" begin
    S = [EData(E(V(10, 10), V(30, 10)), 1, -1),
         EData(E(V(10, 10), V(20, 20)), -1, 1),
         EData(E(V(30, 10), V(40, 20)), 1, -1),
         EData(E(V(20, 20), V(40, 20)), -1, 1)
        ]

    for sp in permutations(S)
        tree = buildsearch(sp, V(0, 0), V(50, 50))
        res = calcpoints(tree)
        @test res[1] == 200
    end
end

@testset "Monotonous flag" begin
    S = [EData(E(V(10, 10), V(30, 10)), 1, -1),
         EData(E(V(10, 10), V(15, 15)), -1, 1),
         EData(E(V(11, 20), V(15, 15)), 1, -1),
         EData(E(V(30, 10), V(40, 20)), 1, -1),
         EData(E(V(11, 20), V(40, 20)), -1, 1)
        ]
    for sp in permutations(S)
        tree = buildsearch(sp, V(0, 0), V(50, 50))
        res = calcpoints(tree)
        @test res[1] == 224
    end
end

@testset "Parallelogram with a hole" begin
    S = [EData(E(V(10, 10), V(30, 10)), 1, -1),
         EData(E(V(10, 10), V(20, 40)), -1, 1),
         EData(E(V(30, 10), V(40, 40)), 1, -1),
         EData(E(V(20, 40), V(40, 40)), -1, 1),
         EData(E(V(14, 12), V(20, 18)), 1, -1),
         EData(E(V(14, 12), V(24, 12)), -1, 1),
         EData(E(V(20, 18), V(24, 12)), 1, -1),
        ]
    for sp in permutations(S)
        tree = buildsearch(sp, V(0, 0), V(50, 50))
        res = calcpoints(tree)
        @test res[1] == 571
    end
end

@testset "Two triangles" begin
    S = [EData(E(V(10, 10), V(30, 10)), 1, -1),
         EData(E(V(10, 10), V(20, 40)), -1, 1),
         EData(E(V(20, 40), V(40, 40)), -1, 2),
         EData(E(V(30, 10), V(40, 40)), 2, -1),
         EData(E(V(20, 40), V(30, 10)), 2, 1)
        ]

    for sp in permutations(S)
        tree = buildsearch(sp, V(0, 0), V(50, 50))
        res = calcpoints(tree)
        @test res[1] == 300
        @test res[2] == 300
    end
end

@testset "Degenerate point" begin
    S = [EData(E(V(10, 10), V(20, 40)), -1, 1),
         EData(E(V(20, 40), V(30, 10)), -1, 1),
         EData(E(V(20, 40), V(30, 45)), -1, 1),
         EData(E(V(20, 40), V(40, 40)), 1, -1),
        ]
    for sp in permutations(S)
        tree = buildsearch(sp, V(0, 0), V(50, 50))
        res = calcpoints(tree)
        @test res[1] == 640
    end
end

@testset "Incomplete Eye" begin
    S = [EData(E(V(50, 90), V(90, 50)), 1, -1),
         EData(E(V(10, 50), V(30, 20)), -1, 1),
         EData(E(V(20, 80), V(50, 90)), 1, -1),
         EData(E(V(30, 20), V(60, 10)), -1, 1),
        ]
    for sp in permutations(S)
        tree = buildsearch(sp, V(0, 0), V(100, 100))
        res = calcpoints(tree)
        @test res[1] == 1094
    end

end

@testset "Bad vertical split" begin
    S = [EData(E(V(10, 10), V(60, 10)), -1, 1),
         EData(E(V(40, 80), V(60, 80)), 1, -1),
         EData(E(V(50, 70), V(60, 80)), -1, 1),
         EData(E(V(20, 50), V(70, 50)), 1, -1),
        ]
    for sp in permutations(S)
        tree = buildsearch(sp, V(0, 0), V(100, 100))
        res = calcpoints(tree)
        @test res[1] == 451
    end
end

end

end #module


text:
using CCDReduction
using Test
using Statistics
using FITSIO

include("data.jl")

function test_header(ccd1::CCDData, ccd2::CCDData)
    header1 = ccd1.hdr
    header2 = ccd2.hdr
    @test keys(header1) == keys(header2)
    for (k1, k2) in zip(keys(header1), keys(header2))
        @test header1[k1] == header2[k2]
    end
end

@testset "basic methods" begin include("methods.jl") end
@testset "FITS interface" begin include("fits.jl") end
@testset "collection" begin include("collection.jl") end


text:
abstract type AbstractReadOnlySet{T,P} <: AbstractSet{T} end

Base.parent(x::AbstractReadOnlySet) = getfield(x, :parent)

# isfrozenvalue(::Type{<:AbstractReadOnlySet}) = true
isfrozenindex(::Type{<:AbstractReadOnlySet}) = true

# meltvalue!(x::AbstractReadOnlySet) = parent(x)
meltindex!(x::AbstractReadOnlySet) = parent(x)

# Workaround: `Base.copymutable(::AbstractSet)` returns an array
ownmutablevalue(x::AbstractSet) = Set(x)
ownmutableindex(x::AbstractSet) = Set(x)

Base.IteratorSize(::Type{<:AbstractReadOnlySet{T,P}}) where {T,P} =
    Base.IteratorSize(P)
Base.IteratorEltype(::Type{<:AbstractReadOnlySet{T,P}}) where {T,P} =
    Base.IteratorEltype(P)
Base.eltype(::Type{<:AbstractReadOnlySet{T,P}}) where {T,P} = eltype(P)
Base.iterate(x::AbstractReadOnlySet, args...) = iterate(parent(x), args...)
Base.length(x::AbstractReadOnlySet) = length(parent(x))

Base.copy(x::AbstractReadOnlySet) = copy(parent(x))
Base.copymutable(x::AbstractReadOnlySet) = Base.copymutable(parent(x))

Base.getproperty(x::AbstractReadOnlySet, name::Symbol) =
    constructorof(typeof(x))(getproperty(parent(x), name))
Base.getproperty(x::AbstractReadOnlySet, field) =
    constructorof(typeof(x))(getproperty(parent(x), field))

Base.broadcastable(x::AbstractReadOnlySet) = Base.broadcastable(parent(x))

struct ReadOnlySet{T,P<:AbstractSet{T}} <: AbstractReadOnlySet{T,P}
    parent::P
end

# isimmutablevalue(::Type{<:ReadOnlySet}) = true
isimmutableindex(::Type{<:ReadOnlySet}) = true
# isfrozenvalue(::Type{<:ReadOnlySet}) = false
isfrozenindex(::Type{<:ReadOnlySet}) = false
readonly_impl(x::AbstractSet) = ReadOnlySet(x)

struct ImmutableSet{T,P<:AbstractSet{T}} <: AbstractReadOnlySet{T,P}
    parent::P
end

freeze!(x::AbstractSet) = ImmutableSet(x)

struct AppendOnlySet{T,P<:AbstractSet{T}} <: AbstractReadOnlySet{T,P}
    parent::P
end

isfrozenindex(::Type{<:AppendOnlySet}) = false
Base.push!(x::AppendOnlySet, args...) = (push!(parent(x), args...); x)

freezevalue!(x::AbstractSet) = AppendOnlySet(x)

factoryof(::ReadOnlySet) = readonly
factoryof(::ImmutableSet) = freeze
factoryof(::AppendOnlySet) = freezevalue

function Base.showarg(io::IO, x::AbstractReadOnlySet, toplevel::Bool)
    print(io, nameof(factoryof(x)), '(')
    Base.showarg(io, parent(x), false)
    print(io, ')')
end


text:
# Copyright (c) 2015: Changhyun Kwon.
# This file is part of https://github.com/chkwon/TrafficAssignment.jl/blob/master/src/load_network.jl

# Traffic Assignment Data structure
mutable struct TA_Data
    network_name::String

    number_of_zones::Int
    number_of_nodes::Int
    first_thru_node::Int
    number_of_links::Int

    init_node::Array{Int,1}
    term_node::Array{Int,1}
    capacity::Array{Float64,1}
    link_length::Array{Float64,1}
    free_flow_time::Array{Float64,1}
    b::Array{Float64,1}
    power::Array{Float64,1}
    speed_limit::Array{Float64,1}
    toll::Array{Float64,1}
    link_type::Array{Int64,1}

    total_od_flow::Float64

    travel_demand::Array{Float64,2}
    od_pairs::Array{Tuple{Int64,Int64},1}

    toll_factor::Float64
    distance_factor::Float64

    best_objective::Float64
end

search_sc(s,c) = something(findfirst(isequal(c), s), 0)

function load_ta_network(network_name, network_data_file, trip_table_file; best_objective=-1.0, toll_factor=0.0, distance_factor=0.0)
    @assert ispath(network_data_file)
    @assert ispath(trip_table_file)

    ##################################################
    # Network Data
    ##################################################
    number_of_zones = 0
    number_of_links = 0
    number_of_nodes = 0
    first_thru_node = 0

    n = open(network_data_file, "r")

    while (line=readline(n)) != ""
        if occursin("<NUMBER OF ZONES>", line)
            number_of_zones = parse(Int, line[ search_sc(line, '>')+1 : end ] )
        elseif occursin("<NUMBER OF NODES>", line)
            number_of_nodes = parse(Int, line[ search_sc(line, '>')+1 : end ] )
        elseif occursin("<FIRST THRU NODE>", line)
            first_thru_node = parse(Int, line[ search_sc(line, '>')+1 : end ] )
        elseif occursin("<NUMBER OF LINKS>", line)
            number_of_links = parse(Int, line[ search_sc(line, '>')+1 : end ] )
        elseif occursin("<END OF METADATA>", line)
            break
        end
    end

    @assert number_of_links > 0

    init_node = Array{Int64}(undef, number_of_links)
    term_node = Array{Int64}(undef, number_of_links)
    capacity = zeros(number_of_links)
    link_length = zeros(number_of_links)
    free_flow_time = zeros(number_of_links)
    b = zeros(number_of_links)
    power = zeros(number_of_links)
    speed_limit = zeros(number_of_links)
    toll = zeros(number_of_links)
    link_type = Array{Int64}(undef, number_of_links)

    idx = 1
    while !eof(n)
      line = readline(n)
        if occursin("~", line) || line == ""
            continue
        end

        if occursin(";", line)
            line = strip(line, [' ', '\n', ';'])
            line = replace(line, ";" => "")

            numbers = split(line)
            init_node[idx] = parse(Int64, numbers[1])
            term_node[idx] = parse(Int64, numbers[2])
            capacity[idx] = parse(Float64, numbers[3])
            link_length[idx] = parse(Float64, numbers[4])
            free_flow_time[idx] = parse(Float64, numbers[5])
            b[idx] = parse(Float64, numbers[6])
            power[idx] = parse(Float64, numbers[7])
            speed_limit[idx] = parse(Float64, numbers[8])
            toll[idx] = parse(Float64, numbers[9])
            link_type[idx] = Int(round(parse(Float64, numbers[10])))

            idx = idx + 1
        end
    end

    ##################################################
    # Trip Table
    ##################################################

    number_of_zones_trip = 0
    total_od_flow = 0

    f = open(trip_table_file, "r")

    while (line=readline(f)) != ""
        if occursin("<NUMBER OF ZONES>", line)
            number_of_zones_trip = parse(Int, line[ search_sc(line, '>')+1 : end ] )
        elseif occursin("<TOTAL OD FLOW>", line)
            total_od_flow = parse(Float64, line[ search_sc(line, '>')+1 : end ] )
        elseif occursin("<END OF METADATA>", line)
            break
        end
    end

    @assert number_of_zones_trip == number_of_zones # Check if number_of_zone is same in both txt files
    @assert total_od_flow > 0

    travel_demand = zeros(number_of_zones, number_of_zones)
    od_pairs = Array{Tuple{Int64, Int64}}(undef, 0)

    origin = -1

    while !eof(f)
        line = readline(f)

        if line == ""
            origin = -1
            continue
        elseif occursin("Origin", line)
            origin = parse(Int, split(line)[2] )
        elseif occursin(";", line)
            pairs = split(line, ";")
            for i=1:size(pairs)[1]
                if occursin(":", pairs[i])
                    pair = split(pairs[i], ":")
                    destination = parse(Int64, strip(pair[1]) )
                    od_flow = parse(Float64, strip(pair[2]) )

                    # println("origin=$origin, destination=$destination, flow=$od_flow")

                    travel_demand[origin, destination] = od_flow
                    push!(od_pairs, (origin, destination))
                end
            end
        end
    end

    # Preparing data to return
    ta_data = TA_Data(
        network_name,
        number_of_zones,
        number_of_nodes,
        first_thru_node,
        number_of_links,
        init_node,
        term_node,
        capacity,
        link_length,
        free_flow_time,
        b,
        power,
        speed_limit,
        toll,
        link_type,
        total_od_flow,
        travel_demand,
        od_pairs,
        toll_factor,
        distance_factor,
        best_objective)

    return ta_data

end # end of load_network function


text:
using Plots

abstract type DescentMethod end

struct GradientDescent <: DescentMethod
    α
end

# init! function does nothing.
init!(M::GradientDescent, f, ∇f, x) = M;

function step!(M::GradientDescent, f, ∇f, x)
    α, g = M.α, ∇f(x);
    return x - α*g;
end

# test

A = rand(2, 2);
b = rand(2);
c = rand(2);
x = rand(2);

f = (x) -> 1/2 * x'*A*x + b'*x .+ c
∇f = (x) -> 1/2 * (A + A') * x + b

p = contour(-10:0.1:10, -10:0.1:10, (x, y)->sum(f([x, y])))
scatter!(p, [x[1]], [x[2]])

m = GradientDescent(0.5);
init!(m, f, ∇f, x);
x_res = [];
y_res = [];

for i in 1:100
    x = step!(m, f, ∇f, x)
    push!(x_res, [x[1]])
    push!(y_res, [x[2]])
end

x_res
y_res

scatter!(p, x_res, y_res)

x_res
y_res


text:
# Wrapper for level-1 BLIS fused-vector routines.
#

@blis_ccall_group(axpy2v,
                  Cvoid,
                  conjx,  BliConj,
                  conjy,  BliConj,
                  m,      BliDim,
                  αx,     Ptr{xType},
                  αy,     Ptr{xType},
                  x,      Ptr{xType}, incx, BliInc,
                  y,      Ptr{xType}, incy, BliInc,
                  z,      Ptr{xType}, incz, BliInc)
@blis_ccall_group(dotaxpyv,
                  Cvoid,
                  conjxt, BliConj,
                  conjx,  BliConj,
                  conjy,  BliConj,
                  m,      BliDim,
                  α,      Ptr{xType},
                  x,      Ptr{xType}, incx, BliInc,
                  y,      Ptr{xType}, incy, BliInc,
                  ρ,      Ptr{xType},
                  z,      Ptr{xType}, incz, BliInc)
@blis_ccall_group(axpyf,
                  Cvoid,
                  conja,  BliConj,
                  conjx,  BliConj,
                  m,      BliDim,
                  b,      BliDim,
                  α,      Ptr{xType},
                  a,      Ptr{xType}, inca, BliInc, lda, BliInc,
                  x,      Ptr{xType}, incx, BliInc,
                  y,      Ptr{xType}, incy, BliInc)
@blis_ccall_group(dotxf,
                  Cvoid,
                  conjat, BliConj,
                  conjx,  BliConj,
                  m,      BliDim,
                  b,      BliDim,
                  α,      Ptr{xType},
                  a,      Ptr{xType}, inca, BliInc, lda, BliInc,
                  x,      Ptr{xType}, incx, BliInc,
                  β,      Ptr{xType},
                  y,      Ptr{xType}, incy, BliInc)
@blis_ccall_group(dotxaxpyf,
                  Cvoid,
                  conjat, BliConj,
                  conja,  BliConj,
                  conjw,  BliConj,
                  conjx,  BliConj,
                  m,      BliDim,
                  b,      BliDim,
                  α,      Ptr{xType},
                  a,      Ptr{xType}, inca, BliInc, lda, BliInc,
                  w,      Ptr{xType}, incw, BliInc,
                  x,      Ptr{xType}, incx, BliInc,
                  β,      Ptr{xType},
                  y,      Ptr{xType}, incy, BliInc,
                  z,      Ptr{xType}, incz, BliInc)


text:
#include("/Users/ericagol/Courses/PlanetaryDynamics/ExoJulia/ExoJulia/Transit/occultquad.jl")
include("occultquad.jl")
include("../../../src/transit_poly_struct.jl")

# Plot accuracy of occultquad vs. transit_poly:

nr = 10000
r=range(0.001, stop=1.0, length=nr)
delta = 1e-12
#u1 = 0.0; u2 = 0.0
u1 = 1.0; u2 = 0.0
#u1 = 2.0; u2 = -1.0
#u1 = 0.5; u2 = 0.5

b1 = 1.0 .- r .- delta
f_MA = occultquad.(b1,u1,u2,r)
f_AL = zeros(Float64,nr); f_AL_big = zeros(Float64,nr)
for i=1:nr
  f_AL[i] = transit_poly(r[i],b1[i],[u1;u2])
  f_AL_big[i] = convert(Float64,transit_poly(big(r[i]),big(b1[i]),[big(u1);big(u2)]))
end

b2 = r .- delta
f_MA2 = occultquad.(b2,u1,u2,r);
f_AL2 = zeros(Float64,nr); f_AL_big2 = zeros(Float64,nr)
for i=1:nr
  f_AL2[i] = transit_poly(r[i],b2[i],[u1;u2])
  f_AL_big2[i] = convert(Float64,transit_poly(big(r[i]),big(b2[i]),[big(u1);big(u2)]))
end

b3 = 1.0 .- r .+ delta
f_MA3 = occultquad.(b3,u1,u2,r)
f_AL3 = zeros(Float64,nr); f_AL_big3 = zeros(Float64,nr)
for i=1:nr
  f_AL3[i] = transit_poly(r[i],b3[i],[u1;u2])
  f_AL_big3[i] = convert(Float64,transit_poly(big(r[i]),big(b3[i]),[big(u1);big(u2)]))
end

b4 = r .+ delta
f_MA4 = occultquad.(b4,u1,u2,r);
f_AL4 = zeros(Float64,nr); f_AL_big4 = zeros(Float64,nr)
for i=1:nr
  f_AL4[i] = transit_poly(r[i],b4[i],[u1;u2])
  f_AL_big4[i] = convert(Float64,transit_poly(big(r[i]),big(b4[i]),[big(u1);big(u2)]))
end
using PyPlot

fig,axes = subplots(2,1, figsize=(8,8))
ax = axes[1]

ax[:plot](r,f_AL,label="ALFM (2019)", color="C0", lw=2)
ax[:plot](r,f_MA,label="Mandel & Agol (2002)", color="C1", lw=1, ls="--")
ax[:plot](r,f_AL2, color="C0", lw=2)
ax[:plot](r,f_MA2, color="C1", lw=1, ls="--")

#ax[:plot](r,f_MA3,label="MA (2002), b=1-r+1e-12")
#ax[:plot](r,f_AL3,label="AL (2018), b=1-r+1e-12")
#ax[:plot](r,f_MA4,label="MA (2002), b=r+1e-12")
#ax[:plot](r,f_AL4,label="AL (2018), b=r+1e-12")
ax[:legend](loc="lower left",fontsize=10)
ax[:set_ylabel]("Flux", fontsize=16)
ax[:axis]([-0.01,1.01,0,1.01])

ax[:annotate](L"$b=1-r-\epsilon$", xy=(0.7, 0.45), xycoords="data",
              xytext=(0, 0), textcoords="offset points",
              ha="center", va="center", rotation=-35, fontsize=10)

ax[:annotate](L"$b=r-\epsilon$", xy=(0.7, 0.68), xycoords="data",
              xytext=(0, 0), textcoords="offset points",
              ha="center", va="center", rotation=-5, fontsize=10)

ax = axes[2]
ax[:semilogy](r,abs.(f_AL-f_AL_big),color="C0", lw=2)
ax[:semilogy](r,abs.(f_MA-f_AL_big),color="C1", lw=1, ls="--")

ax[:semilogy](r,abs.(f_AL2-f_AL_big2),color="C0", lw=2)
ax[:semilogy](r,abs.(f_MA2-f_AL_big2),color="C1", lw=1, ls="--")

#ax[:semilogy](r,abs.(f_MA3-f_AL_big3),label="MA (2002), b=1-r+1e-12")
#ax[:semilogy](r,abs.(f_AL3-f_AL_big3),label="AL (2018), b=1-r+1e-12")
#ax[:semilogy](r,abs.(f_MA4-f_AL_big4),label="MA (2002), b=r+1e-12")
#ax[:semilogy](r,abs.(f_AL4-f_AL_big4),label="AL (2018), b=r+1e-12")
ax[:set_xlabel](L"$r$", fontsize=18)
ax[:set_ylabel]("Error", fontsize=16)
#ax[:legend](loc="left",fontsize=8)
ax[:axis]([-0.01,1.01,1e-18,1e-4])

ax[:annotate](L"$b=1-r-\epsilon$", xy=(0.7,1e-14), xycoords="data",
              xytext=(0, 0), textcoords="offset points",
              ha="center", va="center", rotation=0, fontsize=10)

ax[:annotate](L"$b=r-\epsilon$", xy=(0.7, 1e-8), xycoords="data",
              xytext=(0, 0), textcoords="offset points",
              ha="center", va="center", rotation=0, fontsize=10)

ax[:annotate](L"$both$", xy=(0.7, 1e-17), xycoords="data",
              xytext=(0, 0), textcoords="offset points", color="w",
              ha="center", va="center", rotation=0, fontsize=10)

savefig("compare_MA2002.pdf", bbox_inches="tight")

max_stars_repo_name: dfm/limbdark

max_issues_repo_name: dfm/limbdark

max_forks_repo_name: dfm/limbdark


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

using Markdown
using InteractiveUtils

# ╔═╡ 33c40310-f05d-11eb-066c-87511fe4fa53
begin
	import Pkg
	Pkg.activate(joinpath(@__DIR__))
	using SimLib, Plots, PlutoUI, Statistics, LinearAlgebra, LaTeXStrings, XXZNumerics
	gr()
	LinearAlgebra.BLAS.set_num_threads(8)
	cd(joinpath(@__DIR__, ".."))
end

# ╔═╡ 194bd5e8-3953-44d2-bec0-72e8cba07e81
md"""
### Noisy chain

Consider two atoms with lattice spacing $s$, blockade radius $B$ and position range $\pm \sigma$. So we can say:

$$\begin{align}
	X &\sim \mathcal{U}[-\sigma, \sigma]\\
	Y &\sim \mathcal{U}[s-\sigma, s+\sigma]\\
\end{align}$$

Thus we can calculate the probability of these two atoms interfering to be:

$$\begin{align}
	P(|X-Y|\leq B) = \frac{1}{2σ} \left( \min(B,2\sigma + s) + \max(-B, \min(B,2\sigma-s) \right)
\end{align}$$
"""

# ╔═╡ 4726c0cf-0b39-4bd2-b373-8cd539a3946f
begin
	B=1
	P(s,σ) = (min(B,2σ+s) + max(-B, min(B,2σ-s)))/(4σ)
end

# ╔═╡ 56af5dae-c489-4c7f-9131-1c8223623a6a
let s_range = 0.5:(3B/100):7B,
	σ_range = 0:(3B/50):4B
	heatmap(s_range, σ_range, P.(s_range', σ_range); xlabel="s", ylabel="σ")
end

# ╔═╡ b3105074-35c9-4ddf-b94e-9f92abb506ac
plot(0:(3B/50):4B, 1 .- P.(0:(3B/50):4B, 0:(3B/50):4B))

# ╔═╡ 17b0fd1a-6213-457d-ad42-ae3176e313c1
begin
	function successrate(geom, N; trys=10000)
		count = trys
		for _ in 1:trys
			try
				sample_blockaded(geom, N)
			catch e;
				count -= 1
			end
		end
		count / trys
	end
	function chain_successrate(N, s, σ; trys=1000)
		successrate(NoisyChainPBC(N, s, σ), N)
	end
end

# ╔═╡ c4c39853-6609-4309-9408-07a0cc64fe7c
plot(0:(3B/50):4B, chain_successrate.(Ref(2), 0:(3B/50):4B, 0:(3B/50):4B))

# ╔═╡ 3bc58646-2091-41b3-983a-1354b9a8f3be
successrate(geometry_from_density(:box_pbc, 1.45, 7,1), 7)

# ╔═╡ 709214e0-de51-482d-96d7-84fec5a0f1d7
chain_successrate(10,1.2,0.5)

# ╔═╡ Cell order:
# ╠═33c40310-f05d-11eb-066c-87511fe4fa53
# ╠═194bd5e8-3953-44d2-bec0-72e8cba07e81
# ╠═4726c0cf-0b39-4bd2-b373-8cd539a3946f
# ╠═56af5dae-c489-4c7f-9131-1c8223623a6a
# ╠═b3105074-35c9-4ddf-b94e-9f92abb506ac
# ╠═c4c39853-6609-4309-9408-07a0cc64fe7c
# ╠═17b0fd1a-6213-457d-ad42-ae3176e313c1
# ╠═3bc58646-2091-41b3-983a-1354b9a8f3be
# ╠═709214e0-de51-482d-96d7-84fec5a0f1d7


text:
type FilterDesc
    ptr::Ptr{Void}

    function FilterDesc{T,N}(x::CuArray{T,N}; format=CUDNN_TENSOR_NCHW)
        p = Ptr{Void}[0]
        cudnnCreateFilterDescriptor(p)
        desc = new(p[1])
        finalizer(desc, cudnnDestroyFilterDescriptor)
        csize = Cint[size(x,i) for i=N:-1:1]
        cudnnSetFilterNdDescriptor(desc, datatype(T), format, N, csize)
        desc
    end
end

Base.unsafe_convert(::Type{Ptr{Void}}, desc::FilterDesc) = desc.ptr


text:
import Base: drop

type Drop <: Op; dropout; xdrop; y; Drop(d)=new(d); end

# TODO: implement Dropout using Mul2: however forw conditional on predict?
## maybe we can make rand conditional and mul2 treat nothing as identity?
# TODO: be careful about corrupting the xdrop matrix in Net

params(::Drop)=Any[]
ninputs(::Drop)=1
ysize(::Drop,x)=size(x)
overwrites(::Drop)=true
back_reads_x(::Drop)=false
back_reads_y(::Drop)=false

function forw(l::Drop, x; y=x, train=true, xdrop=nothing, seed=nothing, o...)
    issimilar(x,y) || error("Input mismatch")
    if train && (l.dropout > 0)
        similar!(l, :xdrop, y)
        if xdrop != nothing
            copy!(l.xdrop, xdrop)
        elseif seed != nothing
            # setseed(seed)
            # rand!(l.xdrop)
            # we do this instead of rand! on gpu to eliminate cpu/gpu difference:
            srand(seed)
            copy!(l.xdrop, rand(eltype(l.xdrop), size(l.xdrop)))
        else
            rand!(l.xdrop)
        end
        drop(x, y, l.xdrop, l.dropout, 1/(1-l.dropout))
    else
        y===x || copy!(y,x)
    end
    return y
end

function back(l::Drop, dy; dx=dy, returndx=true, o...)
    returndx || return
    if l.dropout > 0
        @assert issimilar(dy, l.xdrop) "$(summary(dy)) !~ $(summary(l.xdrop)) $(size(l.xdrop))"
        drop(dy, dx, l.xdrop, l.dropout, 1/(1-l.dropout))
    else
        dx===dy || copy!(dx,dy)
    end
    return dx
end

drop(x::Array, y::Array, xdrop::Array, dropout::Number, scale::Number)=(for i=1:length(y); y[i] = (xdrop[i] < dropout ? zero(x[i]) : scale * x[i]); end; y)
drop(x::KUdense, y::KUdense, xdrop::KUdense, dropout::Number, scale::Number)=(drop(x.arr, y.arr, xdrop.arr, dropout, scale); y)

GPU && (drop(x::CudaArray{Float32}, y::CudaArray{Float32}, xdrop::CudaArray{Float32}, dropout, scale)=ccall((:drop32,libknet),Void,(Cint,Ptr{Float32},Ptr{Float32},Ptr{Float32},Cdouble,Cdouble),length(x),x,y,xdrop,dropout,scale))
GPU && (drop(x::CudaArray{Float64}, y::CudaArray{Float64}, xdrop::CudaArray{Float64}, dropout, scale)=ccall((:drop64,libknet),Void,(Cint,Ptr{Float64},Ptr{Float64},Ptr{Float64},Cdouble,Cdouble),length(x),x,y,xdrop,dropout,scale))


text:
"""
    solve_cubic(a, b, c, d)

Finds real roots of a cubic equation
```math
a x^3 + b x^2 + c x + d = 0
```

Return Tuple{T, T, T} where first roots are real,
complex roots are defined to `NaN`

Reference: J. F. Blinn, "How to Solve a Cubic Equation, Part 5: Back to Numerics",
in IEEE Computer Graphics and Applications, vol. 27, no. 3, pp. 78-89, May-June 2007.
doi: 10.1109/MCG.2007.60
"""
function solve_cubic(a, b, c, d)
    # convert to Ax³ + 3Bx² + 3Cx + D = 0
    A, B, C, D = promote(float(a), float(b), float(c), float(d)) ./ (1, 3, 3, 1)
    δ₁ = A * C - B * B
    δ₂ = A * D - B * C
    δ₃ = B * D - C * C
    d13 = δ₁ * δ₃
    d22 = δ₂ * δ₂
    Δ = 4 * d13 - d22
    nanvalue = zero(A) / zero(A)

    if Δ < 0
        At = Cb = Db = zero(A)  # A-tilde, C-bar, D-bar
        if B^3 * D ≥ A * C^3
            At, Cb, Db = A, δ₁, -2 * B * δ₁ + A * δ₂
        else
            At, Cb, Db = D, δ₃, -D * δ₂ + 2 * C * δ₃
        end
        T₀ = -copysign(At, Db) * sqrt(-Δ)
        T₁ = -Db + T₀
        p = cbrt(T₁ / 2)
        q = T₁ == T₀ ? -p : -Cb / p
        x₁ = Cb ≤ 0 ? p + q : -Db / (p^2 + q^2 + Cb)
        x, w = B^3 * D ≥ A * C^3 ? (x₁ - B, A) : (-D, x₁ + C)
        return (x/w, nanvalue, nanvalue)
    else
        δ₁ == δ₂ == δ₃ == 0 && return (-B/A, -B/A, -B/A)
        sΔ = sqrt(Δ)
        θA, θD = abs.(atan.((A*sΔ, 2*B*δ₁ - A*δ₂, D*sΔ, D*δ₂ - 2*C*δ₃)) ./ 3)
        sCA, sCD = sqrt.(.-min.((δ₁, δ₃), 0))
        x₁A, x₁D = 2 .* (sCA, sCD) .* cos.((θA, θD))
        x₃A, x₃D = .-((sCA, sCD)) .* (cos.((θA, θD)) .+ sqrt(3) .* sin.((θA, θD)))
        xlt = (x₁A + x₃A > 2 * B) ? x₁A : x₃A
        xst = (x₁D + x₃D < 2 * C) ? x₁D : x₃D
        xl, wl = xlt - B, A
        xs, ws = -D, xst + C
        E = wl * ws
        F = -xl * ws - wl * xs
        G = xl * xs
        xm, wm = C * F - B * G, C * E - B * F
        return (xs/ws, xm/wm, xl/wl)
    end
end

"""
    solve_cubic!(roots_::AbstractVector, a, b, c, d)

Fill first three items of `roots_` with the real roots of equation or `NaN`s.
```math
a x^3 + b x^2 + c x + d = 0
```
See also: [`solve_cubic`](@ref).
"""
Base.@propagate_inbounds function solve_cubic!(roots_::AbstractVector, a, b, c, d)
    @view(roots_[1:3]) .= solve_cubic(a, b, c, d)
    return roots_
end

max_stars_repo_name: vvpisarev/CubicEoS


text:
using BioMart, GenomicFeatures

# affyids=c("202763_at","209310_s_at","207500_at")
# getBM(attributes=c('affy_hg_u133_plus_2', 'entrezgene_id'),
#       filters = 'affy_hg_u133_plus_2',
#       values = affyids,
#       mart = ensembl)

affyids = ["202763_at", "209310_s_at", "207500_at"]
q = BioMart.Query(
    BioMart.Dataset("hsapiens_gene_ensembl"),
    BioMart.Attributes("affy_hg_u133_plus_2", "entrezgene_id"),
    BioMart.Filter(affy_hg_u133_plus_2 = affyids),
)
q()

# entrez=c("673","837")
# goids = getBM(attributes = c('entrezgene_id', 'go_id'),
#     filters = 'entrezgene_id',
#     values = entrez,
#     mart = ensembl)

entrez = ["673", "837"]
q = BioMart.Query(
    BioMart.Dataset("hsapiens_gene_ensembl"),
    BioMart.Attributes("entrezgene_id", "go_id"),
    BioMart.Filter(entrezgene_id = entrez),
)
q()

# getBM(attributes = c('affy_hg_u133_plus_2','ensembl_gene_id'),
#       filters = c('chromosome_name','start','end'),
#       values = list(16,1100000,1250000),
#       mart = ensembl)

q = BioMart.Query(
    BioMart.Dataset("hsapiens_gene_ensembl"),
    Interval("16", 1100000, 1250000),
    BioMart.Attributes("affy_hg_u133_plus_2", "ensembl_gene_id"),
)
q()


text:
using AxisArrays
using HTTP
using JSON
using JuMP
using MathOptInterface
using MutableArithmetics


text:
using StatsFuns: logsumexp
using AdvancedHMC
using ForwardDiff
using ReverseDiff

"""
    discrete_enumeration(model::MinijyroModel, site_name::Any)

Computes the marginal distribution of `site_name` by enumerating all possible
execution traces. This requires that all sample sites have discrete values.

# Arguments
- `model::MinijyroModel`: model which to perform inference in
- `site_name::Any`: site name which to compute marginal distribution of; use `:_return`
    to compute marginal distribution of the return value

# Returns
- `Distributions.DiscreteNonParametric`: the marginal distribution
"""
function discrete_enumeration(model::MinijyroModel, site_name::Any)
    # Site name is the name of the variable we are interested in.
    empty_trace = Dict()
    empty_trace[:msgs] = Dict()
    q = [empty_trace]

    log_joint!(model)
    enum_model = queue(model, q)
    samples = Dict{Int,Float64}()
    while length(q) > 0
        trace = enum_model()
        val = site_name == :_return ? trace[:_return] : trace[:msgs][site_name][:value]
        if !haskey(samples, val)
            samples[val] = 0.0
        end
        samples[val] += exp(trace[:logjoint])
    end

    # Normalize scores so that they sum to 1.
    vals = collect(keys(samples))
    probs = collect(values(samples))
    probs = probs ./ sum(probs)
    return DiscreteNonParametric(vals, probs)
end

function get_param_info(model::MinijyroModel, model_args::Tuple)
    # For each unobserved site get information on whether it is a scalar or not
    # and if it is not a scalar then get its shape and the number of elements in
    # there.
    t = trace(model)(model_args...)
    param_info = Dict()
    num_params = 0 # Total number of parameters
    for msg in values(t[:msgs])
        if !msg[:is_observed]
            val = msg[:value]
            d = Dict{Symbol,Any}(
                :is_scalar => isa(val, Number)
            )
            if !d[:is_scalar]
                @assert isa(val, AbstractArray) "Sampled values must either be Number or AbstractArray but found $(typeof(val))."
                d[:shape] = size(val)
                d[:num_elems] = prod(size(val))
            end
            param_info[msg[:name]] = d
            num_params += d[:is_scalar] ? 1 : d[:num_elems]
        end
    end

    return (param_info, num_params)
end

function load_params!(params_dict::Dict, param_info::Dict, params::AbstractArray)
    # Load parameters from params into params_dict.
    i = 1
    for name in keys(param_info)
        if param_info[name][:is_scalar]
            params_dict[name] = params[i]
            i += 1
        else
            num_elems = param_info[name][:num_elems]
            params_dict[name] = reshape(params[i:i+num_elems-1], param_info[name][:shape])
            i += num_elems
        end
    end
end

function get_grad_fn(density::Function, autodiff::Symbol)
    if autodiff == :reverse
        function grad_density(x)
            results = DiffResults.GradientResult(x)
            ReverseDiff.gradient!(results, density, x)
            return (DiffResults.value(results), DiffResults.gradient(results))
        end
        return grad_density
    elseif autodiff == :forward
        return ForwardDiff
    else
        error("Autodiff option $autodiff not supported")
    end
end

"""
    hmc(args...; kwargs...)

Uses `AdvancedHMC` to run static Hamiltonian Monte Carlo with a fixed integrator
step size and fixed number of integrator steps.

# Arguments
- `model::MinijyroModel`: model to do inference in
- `model_args::Tuple`: arguments for the `model` to be called with
- `step_size::Float64`: integrator step size
- `n_leapfrog::Int`: number of steps taken by leapfrog integrator
- `n_samples::Int`: number of samples to take from posterior

# Keywords
- `autodiff::Symbol=:forward`: either `:forward` or `:reverse` to select autodiff mode
"""
function hmc(
    model::MinijyroModel,
    model_args::Tuple,
    step_size::Float64,
    n_leapfrog::Int,
    n_samples::Int;
    autodiff::Symbol=:forward
)
    # NOTE: This assumes a "static" model with fixed dimensionality.
    param_info, num_params = get_param_info(model, model_args)
    initial_params = randn(num_params)

    function density(params)
        params_dict = Dict()
        load_params!(params_dict, param_info, params)
        m = log_joint(condition(model, params_dict))
        return m(model_args...)[:logjoint]
    end

    grad_fn = get_grad_fn(density, autodiff)

    # Construct Hamiltonian system
    metric = UnitEuclideanMetric(length(initial_params))
    hamiltonian = Hamiltonian(metric, density, grad_fn)

    integrator = Leapfrog(step_size)
    proposal = StaticTrajectory(integrator, n_leapfrog)
    return AdvancedHMC.sample(hamiltonian, proposal, initial_params, n_samples)
end

"""
    nuts(args...; kwargs...)

Uses `AdvancedHMC` to run the NUTS algorithm to get posterior samples from `model`.

# Arguments
- `model::MinijyroModel`: model to do inference in
- `model_args::Tuple`: arguments for the `model` to be called with
- `num_samples::Int`: number of samples to take from posterior

# Keywords
- `num_warmup::Int=1000`: number of warmup steps to burn-in the chain; these samples
    will be discarded
- `autodiff::Symbol=:forward`: either `:forward` or `:reverse` to select autodiff mode
- `progress::Bool=true`: whether to display progress meter or not
"""
function nuts(
    model::MinijyroModel,
    model_args::Tuple,
    num_samples::Int;
    num_warmup::Int=1000,
    autodiff::Symbol=:forward,
    progress::Bool=true
)
    # NOTE: This assumes a "static" model with fixed dimensionality.
    param_info, num_params = get_param_info(model, model_args)
    initial_params = randn(num_params)

    function density(params)
        params_dict = Dict()
        load_params!(params_dict, param_info, params)
        m = log_joint(condition(model, params_dict))
        return m(model_args...)[:logjoint]
    end

    grad_fn = get_grad_fn(density, autodiff)

    metric = DiagEuclideanMetric(length(initial_params))
    hamiltonian = Hamiltonian(metric, density, grad_fn)

    initial_step_size = find_good_eps(hamiltonian, initial_params)
    integrator = Leapfrog(initial_step_size)

    proposal = NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)
    adaptor = StanHMCAdaptor(
        Preconditioner(metric), NesterovDualAveraging(0.8, integrator)
    )

    return AdvancedHMC.sample(
        hamiltonian,
        proposal,
        initial_params,
        num_samples,
        adaptor,
        num_warmup;
        progress=progress
    )
end


text:
# julia 1.0

# Import the Halite SDK, which will let you interact with the game.
include("hlt/Halite.jl")
using Halite

""" <<<Game Begin>>> """
# This game object contains the initial game state.
game = Game()
# At this point "game" variable is populated with initial map data.
# This is a good place to do computationally expensive start-up pre-processing.
# As soon as you call "ready" function below, the 2 second per turn timer will start.
ready(game, "MyJuliaBot")

# Now that your bot is initialized, save a message to yourself in the log file with some important information.
#   Here, you log here your id, which you can always fetch from the game object by using my_id.
info("Successfully created bot! My Player ID is $(game.my_id)")

""" <<<Game Loop>>> """

while true
    # This loop handles each turn of the game. The game object changes every turn, and you refresh that state by
    #   running update_frame().
    update_frame(game)
    # You extract player metadata and the updated map metadata here for convenience.
    me = game.me
    game_map = game.game_map

    # A command queue holds all the commands you will run this turn. You build this list up and submit it at the
    #   end of the turn.
    command_queue = []

    for ship in get_ships(me)
        # For each of your ships, move randomly if the ship is on a low halite location or the ship is full.
        #   Else, collect halite.
        if game_map[ship.position].halite_amount < MAX_HALITE / 10 || ship.is_full:
            append!(command_queue,move(ship,
                    random.choice([ Direction.North, Direction.South, Direction.East, Direction.West ])))
        else:
            append!(command_queue, stay_still(ship))
        end
    end

    # If the game is in the first 200 turns and you have enough halite, spawn a ship.
    # Don't spawn a ship if you currently have a ship at port, though - the ships will collide.
    if game.turn_number <= 200 && me.halite_amount >= SHIP_COST && !game_map[me.shipyard].is_occupied:
        append!(command_queue, me.shipyard.spawn())
    end

    # Send your moves back to the game environment, ending this turn.
    end_turn(game, command_queue)
end

max_stars_repo_name: Affos/Halite-III

max_issues_repo_name: Affos/Halite-III

max_forks_repo_name: Affos/Halite-III


text:
"""
    Ricker(; <keyword arguments>)

Create a Ricker wavelet.

# Keyword arguments
* `dt::Real=0.002`: sampling interval in secs.
* `f0::Real=20.0`: central frequency in Hz.

# Examples
```julia
julia> w = Ricker(); plot(w);
julia> w = Ricker(dt=0.004, f0=20); plot(w);
```

# Reference
Sheriff, Robert, 2002, Encyclopedic Dictionary of Applied Geophysics, fourth
ed.: Society of Exploration Geophysicists. Geophysical Reference Series No. 13.
"""

function Ricker(; dt::Real=0.002, f0::Real=20.0)

    nw = 2.0/(f0*dt)
    nc = floor(Int, nw/2)
    t = dt*collect(-nc:1:nc)
    b = (pi*f0*t).^2
    w = (1.-2.*b).*exp(-b)

end


text:
module LeafPhotosynthesisMod
using MathToolsMod
using Parameters
using Leaf

export fluxes
"Tolerance thhreshold for Ci iterations"
tol = 0.1
vpd_min = 0.1

" Just a placeholder for now"
@with_kw mutable struct fluxes{TT<:Number}
         APAR::TT = 500
         gbc::TT = 100
         gbv::TT = 100
         cair::TT = 400
         ceair::TT = 1400
         eair::TT = 1400
         je::TT = 1100
         #gs::TT = 0.0
         ac::TT = 0
         aj::TT = 0
         ai::TT = 0
         ap::TT = 0
         ag::TT = 0
         an::TT = 0
         cs::TT = 0
         #ci::TT = 0
         rd::TT = 0
         Je_pot::TT = 0
         Ja::TT = 0
         Je_red::TT = 0
         φ::TT = 0
end

struct atmos
         o2air              # Atmospheric O2 (mmol/mol)
         co2air             # Atmospheric CO2 (μmol/mol)
         eair               # Vapor pressure of air (Pa)
end

"""
    LeafPhotosynthesis(flux::fluxes, leaf::leaf_params,T::Number)

Compute net assimilation rate A, fluorescence F using biochemical model

# Arguments
- `flux::fluxes`: fluxes structure.
- `leaf::leaf_params`: leaf_params structure.
- `T::Number`: Leaf Temperature
"""
function LeafPhotosynthesis(flux::fluxes, leaf::leaf_params,T::Number)
    # Adjust rates to leaf Temperature (C3 only for now):
    setLeafT!(leaf, T)
    # Compute max PSII efficiency here (can later be used with a variable Kn!)
    if !leaf.dynamic_state
        leaf.Kn = 0.0
    end
    leaf.Kp = 4.0

    φ_PSII = leaf.Kp/(leaf.Kp+leaf.Kf+leaf.Kd+leaf.Kn)

    # Save leaf respiration
    flux.rd = leaf.rdleaf
    # Calculate potential electron transport rate (assuming no upper bound, proportional to absorbed light!):
    flux.Je_pot = 0.5 * leaf.maxPSII * flux.APAR;                          # potential electron transport rate (important for later)
    flux.Je_red = 0.5 * φ_PSII * flux.APAR;                                # Includes Kn here
    # Some bound constraint on VPD:
    flux.ceair = min(max(flux.eair, 0.03*leaf.esat), leaf.esat )

    # Electron transport rate for C3 plants
    # Actual colimited potential Je (curvature and Jmax)
    #flux.je = minimum(quadratic(leaf.θ_j, -(flux.Je_red + leaf.jmax), flux.Je_red * leaf.jmax))    # Bonan eq. 11.21
    flux.je = min(flux.Je_red,leaf.jmax)
    #flux.je = minimum(quadratic(leaf.θ_j, -(flux.Je_pot + leaf.jmax), flux.Je_pot * leaf.jmax))    # Bonan eq. 11.21

    # Ci calculation
    # Medlyn or Ball-Berry:
    if leaf.dynamic_state # Save actual gs
        gs_actual = leaf.gs
    end

    if (leaf.gstyp <= 1)
        Ci_0 = leaf.C3 ? 0.7*flux.cair : 0.4*flux.cair
        # Solve iterative loop:
        leaf.Ci = hybrid(flux,leaf, CiFunc!, Ci_0, 0.7*Ci_0, tol)
    elseif leaf.gstyp == 2 # Needed for Bonan Stomatal optimization model
        leaf.Ci = CiFuncGs!(leaf.gs, flux,leaf)
    end
    if leaf.dynamic_state
        leaf.gs_ss = leaf.gs
        leaf.gs = gs_actual
        leaf.Ci = CiFuncGs!(leaf.gs, flux,leaf)
    end

    # Rate of actual CO2 per electron, incl. photorespiration (not using effcon here for now)
    leaf.CO2_per_electron = (leaf.Ci-leaf.Γstar)/(4leaf.Ci+8leaf.Γstar) #* leaf.effcon;

    # Actual effective ETR:
    flux.Ja = max(0,flux.ag / leaf.CO2_per_electron);
    flux.Ja = min(flux.Ja,flux.Je_pot )

    # Effective photochemical yield:
    flux.φ = leaf.maxPSII*flux.Ja/flux.Je_pot;
    #println(flux.Ja, " ", flux.Je_pot)
    flux.φ = min(1/leaf.maxPSII,flux.φ)
    x   = max(0,  1-flux.φ/leaf.maxPSII);       # degree of light saturation: 'x' (van der Tol e.a. 2014)
    Fluorescencemodel!(flux.φ,x,leaf)

end # LeafPhotosynthesis (similar to biochem in SCOPE)

"""
    CiFunc!(Ci::Number, flux::fluxes, leaf::leaf_params)

Compute Assimilation using Ci as input

# Arguments
- `Ci::Number`: Ci.
- `flux::fluxes`: fluxes structure.
- `leaf::leaf_params`: leaf_params structure.
"""
function CiFunc!(Ci::Number, flux::fluxes, leaf::leaf_params)

    if leaf.C3
        # C3: Rubisco-limited photosynthesis; still need to check CO2 mixing ratios vs partial pressures.
        # still need to include ppm2bar (but can be done on leaf structure!)
        flux.ac = leaf.vcmax * max(Ci-leaf.Γstar, 0.0) / (Ci + leaf.kc*(1.0+leaf.o₂/leaf.ko)) # Bonan eq. 11.28
        # C3: RuBP-limited photosynthesis (this is the NADPH requirement stochiometry)
        flux.aj = flux.je * max(Ci-leaf.Γstar, 0.0) / (4.0*Ci + 8.0*leaf.Γstar)               # Bonan eq. 11.29

        # for C3, set ap to Inf
        flux.ap = Inf
    else #C4 Photosynthesis, still to be implemented
        flux.ac = flux.aj = flux.ap = 0.0
    end
    # Net photosynthesis as the minimum or co-limited rate
    if leaf.use_colim
        flux.ai = minimum(quadratic(leaf.C3 ? 0.98 : 0.80, -(flux.ac + flux.aj), flux.ac * flux.aj))
        if leaf.C3
            flux.ag = flux.ai
        else # C4 colimitation with ap
            flux.ag = minimum(quadratic(0.95, -(flux.ai + flux.ap), flux.ai*flux.ap))                 # Bonan Eq 11.33
        end
    else
        flux.ag = min(flux.ac,flux.aj,flux.ap)
    end
    # Prevent photosynthesis from ever being negative
    flux.ag = max(0,flux.ag)
    flux.ai = max(0,flux.ai)
    flux.aj = max(0,flux.aj)
    flux.ap = max(0,flux.ap)

    # Net photosynthesis
    flux.an = flux.ag - leaf.rdleaf

    # CO2 at leaf surface # might need to be changed
    flux.cs = flux.cair - flux.an / flux.gbc

    # Stomatal constraint function (not sure we "need" the quadratic colimitations here, why not just use BB or Medlyn?)
    if (leaf.gstyp == 1) # Ball-Berry
        if flux.an >0.0
            leaf.gs = maximum(quadratic(flux.cs, flux.cs*(flux.gbv - leaf.g0) - leaf.g1*flux.an, -flux.gbv * (flux.cs*leaf.g0 + leaf.g1*flux.an*flux.ceair/leaf.esat)))
            #leaf.g1 * flux.an * flux.ceair/leaf.esat/flux.cs  + leaf.g0;
            # println(leaf.gs)
        else
            leaf.gs = leaf.g0
        end
    elseif (leaf.gstyp == 0) # Medlyn
        if flux.an >0.0
            # Not sure how this all works, copied from Bonan's ML canopy model
            vpd_term = max((leaf.esat - flux.ceair), vpd_min) * 0.001
            term = 1.6 * flux.an / flux.cs
            leaf.gs = maximum(quadratic(1.0, -(2.0 * (leaf.g0 + term) + (leaf.g1 * term)^2 / (flux.gbv * vpd_term)), leaf.g0 * leaf.g0 + (2.0 * leaf.g0 + term * (1.0 - leaf.g1 * leaf.g1 / vpd_term)) * term))
        else
            leaf.gs = leaf.g0
        end
    end
    # Diffusion (supply-based) photosynthetic rate - Calculate Ci from the diffusion rate
    gleaf = 1.0 / (1.0/flux.gbc + 1.6/leaf.gs + 1.0/leaf.gm)
    cinew = flux.cair - flux.an / gleaf

    # CiFunc returns the difference between the current Ci and the new Ci
    leaf.Ci = cinew
    #return flux.an<0. ? 0.0 : cinew - Ci
    return cinew - Ci
end

"""
    CiFuncGs!(gs::Number, flux::fluxes, leaf::leaf_params)

Compute Assimilation using fixed stomatal conductance gs.

# Arguments
- `gs::Number`: Stomatal conductance.
- `flux::fluxes`: fluxes structure.
- `leaf::leaf_params`: leaf_params structure.
"""
function CiFuncGs!(gs::Number, flux::fluxes, leaf::leaf_params)
    # Compute overall conductance (Boundary layer, stomata and mesophyll)
    gleaf = 1.0/(1.0/flux.gbc + 1.6/gs + 1.0/leaf.gm)
    if gleaf<eps() gleaf=eps() end

    #flux.ac = leaf.vcmax * max(Ci-leaf.Γstar, 0.0) / (Ci + leaf.kc*(1.0+leaf.o₂/leaf.ko)) # Bonan eq. 11.28
    # C3: RuBP-limited photosynthesis (this is the NADPH requirement stochiometry)
    #flux.aj = flux.je * max(Ci-leaf.Γstar, 0.0) / (4.0*Ci + 8.0*leaf.Γstar)               # Bonan eq. 11.29

    if leaf.C3
        # C3 Rubisco Limited Photosynthesis co-limited by gs
        a0 = leaf.vcmax
        e0 = 1.0
        d0 = leaf.kc*(1.0+leaf.o₂/leaf.ko)
        flux.ac = minimum(quadratic(1.0/gleaf, -(e0*flux.cair + d0) - (a0 - e0*leaf.rdleaf) / gleaf, a0 * (flux.cair - leaf.Γstar) - leaf.rdleaf * (e0*flux.cair + d0)))+leaf.rdleaf

        # C3: RuBP-limited photosynthesis
        a0 = flux.je
        e0 = 4.0
        d0 = 8.0*leaf.Γstar
        flux.aj = minimum(quadratic(e0 / gleaf, -(e0*flux.cair + d0) - (a0 - e0*leaf.rdleaf) / gleaf, a0 * (flux.cair - leaf.Γstar) - leaf.rdleaf * (e0*flux.cair + d0)))+leaf.rdleaf

        # C3: Product-limited photosynthesis
        flux.ap = Inf
    # C4 to be implemented
    elseif !leaf.C3
        flux.ac = flux.aj = flux.ap = 0.0
    end
    if leaf.use_colim
        flux.ai = minimum(quadratic(leaf.C3 ? 0.98 : 0.80, -(flux.ac + flux.aj), flux.ac * flux.aj))   # Bonan Eq 11.33
        # Ap limitation only for C4 here:
        if leaf.C3
            flux.ag = flux.ai
        else # C4 colimitation with ap
            flux.ag = minimum(quadratic(0.95, -(flux.ai + flux.ap), flux.ai*flux.ap))                  # Bonan Eq 11.33
        end
    else
        flux.ag = min(flux.ac,flux.aj,flux.ap)
    end
    flux.ag = max(0,flux.ag)
    flux.ai = max(0,flux.ai)
    flux.aj = max(0,flux.aj)
    flux.ap = max(0,flux.ap)
    # Compute net Photosynthesis
    flux.an = flux.ag - leaf.rdleaf
    # Compute CO2 at leaf surface
    flux.cs = flux.cair - flux.an / flux.gbc

    # Compute Ci (included Mesophyll as well in principle)
    ci_val = flux.cair - flux.an / gleaf
    #leaf.CO2_per_electron = (ci_val-leaf.Γstar)./(ci_val+2.0*leaf.Γstar) .* leaf.effcon;
end # Function CiFuncGs!

"""
    Fluorescencemodel!(ps,x,leaf::leaf_params )

Compute Fluorescence yields, Kn and Kp.

# Arguments
- `ps::Number`: PSII yield.
- `x::Number`: Degree of light saturation: [0-1] .
- `leaf::leaf_params`: leaf_params structure.
"""
function Fluorescencemodel!(ps::Number,x::Number,leaf::leaf_params )
    const Kp_max = 4.0
    x_alpha = exp(log(x)*leaf.Knparams[2]); # this is the most expensive operation in this fn; doing it twice almost doubles the time spent here (MATLAB 2013b doesn't optimize the duplicate code)
    #println(x_alpha)
    leaf.Kn_ss = leaf.Knparams[1] * (1+leaf.Knparams[3])* x_alpha/(leaf.Knparams[3] + x_alpha);
    if !leaf.dynamic_state
        leaf.Kn = leaf.Kn_ss
    end
    Kf = leaf.Kf
    Kn = leaf.Kn
    Kd = leaf.Kd
    leaf.Kp   = min(max(0,-ps*(Kf+Kd+Kn)/(ps-1)),Kp_max);
    Kp = leaf.Kp

    leaf.Fo   = Kf/(Kf+Kp_max+Kd);
    leaf.Fo′  = Kf/(Kf+Kp_max+Kd+Kn);
    leaf.Fm   = Kf/(Kf   +Kd);
    leaf.Fm′  = Kf/(Kf   +Kd+Kn);
    leaf.ϕs   = leaf.Fm′*(1-ps);
    leaf.eta  = leaf.ϕs/leaf.Fo;
    leaf.qQ   = 1-(leaf.ϕs-leaf.Fo′)/(leaf.Fm′-leaf.Fo′);
    leaf.qE   = 1-(leaf.Fm-leaf.Fo′)/(leaf.Fm-leaf.Fo);
    leaf.NPQ  = Kn/(Kf+Kd);
end

end #Module

max_stars_repo_name: ruillercoaster/LSM-SPAM

max_issues_repo_name: ruillercoaster/LSM-SPAM

max_forks_repo_name: ruillercoaster/LSM-SPAM


text:
include("Include.jl")

# extra:
using Flux
using Flux: @epochs
using BSON: @save

# load the data -
path_data_file_TF = joinpath(_PATH_TO_DATA, "Training-Thrombin-TF.csv")
full_data_table_TF = load(path_data_file_TF)
path_data_file_TF_TM = joinpath(_PATH_TO_DATA, "Training-Thrombin-TF-TM.csv")
full_data_table_TF_TM = load(path_data_file_TF_TM)
full_training_data_frame = vcat(full_data_table_TF, full_data_table_TF_TM)

# filter -
has_TM_flag = 0
experimental_data_table = filter([:visitid, :TM] => (x, y) -> (x == 2 || x == 3), full_training_data_frame)

# get input and output data -
input_data = convert.(Float32, Matrix(experimental_data_table[!, 3:14]))
output_data = convert.(Float32, Matrix(experimental_data_table[!, 15:19]))
training_data = [(transpose(input_data), transpose(output_data))]

# build a model architecture -
deep_coag_model = Chain(Dense(12, 12, σ), Dense(12, 5));

# setup a loss function -
loss(x, y) = Flux.Losses.mae(deep_coag_model(x), y; agg = mean)

# pointer to params -
ps = Flux.params(deep_coag_model)

# # use old school gradient descent -
opt = Momentum(0.25, 0.95)

# # train -
@epochs 10000 Flux.train!(loss, ps, training_data, opt)

# save -
model_file_path = joinpath(_PATH_TO_MODELS, "deep_coag_model.bson")
@save model_file_path deep_coag_model

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
module Polygons

using Compat
using Compat.Statistics: mean
using Compat: reverse
import Base: length, show, isinf

export Polygon,vertex,interiorangle,isinpoly,naca4

"""
    Polygon(x::Vector{Float64}, y::Vector{Float64})

A polygon defined by its vertices, which must be provided in
counter-clockwise order.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0])
Polygon with 4 vertices at
             (-1.0,-1.0) (0.2,-1.0) (1.0,0.5) (-1.0,1.0)
             interior angles/π = [0.5, 0.656, 0.422, 0.422]
```
"""
struct Polygon
  vert :: Vector{ComplexF64}
  angle :: Vector{Float64}

  Polygon(vert,angle) = abs(vert[end]-vert[1])<eps() ? new(vert[1:end-1],angle) : new(vert,angle)
end

Polygon(x::T,y::T,angle) where T = Polygon(x+im*y,angle)

Polygon(x::T,y::T) where T = Polygon(x+im*y,interiorangle(x+im*y))

"""
    Polygon(w::Vector{Complex128})

Sets up a polygon with the coordinates of the vertices specified
with complex vector `w`. As usual, these must be supplied in
counter-clockwise order.

# Example

```jldoctest
julia> p = Polygon([-1.0-1.0im,0.2-1.0im,1.0+0.5im,-1.0+1.0im])
Polygon with 4 vertices at
             (-1.0,-1.0) (0.2,-1.0) (1.0,0.5) (-1.0,1.0)
             interior angles/π = [0.5, 0.656, 0.422, 0.422]
```
"""
Polygon(w::Vector{ComplexF64}) = Polygon(w,interiorangle(w))

"""
    vertex(p::Polygon) -> Vector{Complex128}

Returns the vector of vertices of the polygon `p`, in complex form.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> vertex(p)
4-element Array{Complex{Float64},1}:
 -1.0-1.0im
  0.2-1.0im
  1.0+0.5im
 -1.0+1.0im
```
"""
vertex(p::Polygon) = p.vert

"""
    isinf(p::Polygon) -> Bool

Returns `true` if any vertex in polygon `p` is at infinity.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> isinf(p)
false
```
"""
Base.isinf(p::Polygon) = any(isinf.(vertex(p)))

"""
    length(p::Polygon) -> Integer

Returns the number of vertices of the polygon `p`.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> length(p)
4
```
"""
Base.length(p::Polygon) = length(vertex(p))

"""
    interiorangle(p::Polygon) -> Vector{Float64}

Returns the vector of interior angles (divided by \$\\pi\$) of the polygon `p`.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> interiorangle(p)
4-element Array{Float64,1}:
 0.5
 0.655958
 0.422021
 0.422021
```
"""
interiorangle(p::Polygon) = length(p.angle) != 0 ? p.angle : interiorangle(p.vertex)

function interiorangle(w::Vector{ComplexF64})
  if length(w)==0
    return []
  end

  atinf = isinf.(w)
  mask = .~(atinf .| circshift(atinf,-1) .| circshift(atinf,1))

  dw = w - circshift(w,1)
  dwshift = circshift(dw,-1)
  beta = fill(NaN,length(w))
  beta[mask] = mod.(angle.( -dw[mask].*conj.(dwshift[mask]) )/π,2)

  mods = abs.(beta .+ 1) .< 1e-12
  beta[mods] = fill!(similar(beta[mods]), 1)

  return beta

end

function isinpoly(z::ComplexF64,w::Vector{ComplexF64},beta::Vector{Float64},tol)

  index = 0.0

  scale = mean(abs.(diff(circshift(w,-1))))
  if ~any(scale > eps())
        return
  end
  w = w/scale
  z = z/scale
  d = w .- z
  d[abs.(d) .< eps()] .= eps()
  ang = angle.(circshift(d,-1)./d)/π
  tangents = sign.(circshift(w,-1)-w)

  # Search for repeated points and fix these tangents
  for p = findall( tangents .== 0 )
    v = [w[p+1:end];w]
    tangents[p] = sign(v[findfirst(v.!=w[p])]-w[p])
  end

  # Points close to an edge
  onbdy = abs.(imag.(d./tangents)) .< 10*tol
  onvtx = abs.(d) .< tol
  onbdy = onbdy .& ( (abs.(ang) .> 0.9) .| onvtx .| circshift(onvtx,-1) )

  if ~any(onbdy)
    index = round(sum(ang)/2)
  else
    S = sum(ang[.~onbdy])
    b = beta[onvtx]
    augment = sum(onbdy) - sum(onvtx) - sum(b)

    index = round(augment*sign(S) + S)/2
  end
  return index==1

end

isinpoly(z,w,beta) = isinpoly(z,w,beta,eps())

"""
    isinpoly(z::Complex128,p::Polygon) -> Bool

Returns `true` or `false` depending on whether `z` is inside
or outside polygon `p`.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> isinpoly(0.0+0.0im,p)
true

julia> isinpoly(1.0+2.0im,p)
false
```
"""
isinpoly(z,p::Polygon) = isinpoly(z,p::Polygon,eps())

"""
    isinpoly(z::Complex128,p::Polygon,tol::Float64) -> Bool

Returns `true` if `z` is inside or within distance `tol` of polygon `p`.

# Example

```jldoctest
julia> p = Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0]);

julia> isinpoly(-1.01+0.0im,p)
false

julia> isinpoly(-1.01+0.0im,p,1e-2)
true
```
"""
isinpoly(z,p::Polygon,tol) = isinpoly(z,p.vert,p.angle,tol)

winding(z,x...) = float.(isinpoly(z,x...))

#=  some specific shape families =#

"""
    naca4(cam,pos,t[;np=20][,Zc=0.0+0.0im][,len=1.0]) -> Vector{Complex128}

Generates the vertices of a NACA 4-digit airfoil of chord length 1. The
relative camber is specified by `cam`, the position of
maximum camber (as fraction of chord) by `pos`, and the relative thickness
by `t`.

The optional parameter `np` specifies the number of points on the upper
or lower surface. The optional parameter `Zc` specifies the mean position of
the vertices (which is set to the origin by default). The optional parameter
`len` specifies the chord length.

# Example

```jldoctest
julia> w = naca4(0.0,0.0,0.12);

julia> p = Polygon(w);
```
"""
function naca4(cam::Number,pos::Number,t::Number;np=20,Zc=0.0+0.0im,len=1.0)

# Here, cam is the fractional camber, pos is the fractional chordwise position
# of max camber, and t is the fractional thickness.

npan = 2*np-2

# Trailing edge bunching
an = 1.5
anp = an+1
x = zeros(np)

θ = zero(x)
yc = zero(x)

for j = 1:np
    frac = Float64((j-1)/(np-1))
    x[j] = 1 - anp*frac*(1-frac)^an-(1-frac)^anp;
    if x[j] < pos
        yc[j] = cam/pos^2*(2*pos*x[j]-x[j]^2)
        if pos > 0
            θ[j] = atan(2*cam/pos*(1-x[j]/pos))
        end
    else
        yc[j] = cam/(1-pos)^2*((1-2*pos)+2*pos*x[j]-x[j]^2)
        if pos > 0
            θ[j] = atan(2*cam*pos/(1-pos)^2*(1-x[j]/pos))
        end
    end
end

xu = zero(x)
yu = xu
xl = xu
yl = yu

yt = t/0.20*(0.29690*sqrt.(x)-0.12600*x-0.35160*x.^2+0.28430*x.^3-0.10150*x.^4)

xu = x-yt.*sin.(θ)
yu = yc+yt.*cos.(θ)

xl = x+yt.*sin.(θ)
yl = yc-yt.*cos.(θ)

rt = 1.1019*t^2;
θ0 = 0
if pos > 0
    θ0 = atan(2*cam/pos)
end
# Center of leading edge radius
xrc = rt*cos(θ0)
yrc = rt*sin(θ0)
θle = collect(0:π/50:2π)
xlec = xrc .+ rt*cos.(θle)
ylec = yrc .+ rt*sin.(θle)

# Assemble data
coords = [xu yu xl yl x yc]
cole = [xlec ylec]

# Close the trailing edge
xpanold = [0.5*(xl[np]+xu[np]); reverse(xl[2:np-1], dims = 1); xu[1:np-1]]
ypanold = [0.5*(yl[np]+yu[np]); reverse(yl[2:np-1], dims = 1); yu[1:np-1]]

xpan = zeros(npan)
ypan = zeros(npan)
for ipan = 1:npan
    if ipan < npan
        xpan1 = xpanold[ipan]
        ypan1 = ypanold[ipan]
        xpan2 = xpanold[ipan+1]
        ypan2 = ypanold[ipan+1]
    else
        xpan1 = xpanold[npan]
        ypan1 = ypanold[npan]
        xpan2 = xpanold[1]
        ypan2 = ypanold[1]
    end
    xpan[ipan] = 0.5*(xpan1+xpan2)
    ypan[ipan] = 0.5*(ypan1+ypan2)
end
w = ComplexF64[1;reverse(xpan, dims = 1)+im*reverse(ypan,dims = 1)]*len
w .-= sum(w)/length(w)
return w .+ Zc

end

function show(io::IO, p::Polygon)
    println(io, "Polygon with $(length(p.vert)) vertices at")
    print(io,   "             ")
    for i = 1:length(p.vert)
        print(io,"($(real(p.vert[i])),$(imag(p.vert[i]))) ")
    end
    println(io)
    println(io, "             interior angles/π = $(round.(p.angle, digits=3))")
end

end


text:
#################### Raftery and Lewis Diagnostic ####################

function rafterydiag(x::Vector{T}; q::Real=0.025, r::Real=0.005,
                      s::Real=0.95, eps::Real=0.001,
                      range::AbstractRange=1:1:length(x)) where {T<:Real}
  nx = length(x)
  phi = sqrt(2.0) * erfinv(s)
  nmin = ceil(Int, q * (1.0 - q) * (phi / r)^2)
  if nmin > nx
    @warn("At least $nmin samples are needed for specified q, r, and s")
    kthin = burnin = total = NaN
  else
    dichot = Int[(x .<= quantile(x, q))...]
    kthin = 0
    bic = 1.0
    local test, ntest
    while bic >= 0.0
      kthin += 1
      test = dichot[1:kthin:nx]
      ntest = length(test)
      temp = test[1:(ntest - 2)] + 2 * test[2:(ntest - 1)] + 4 * test[3:ntest]
      trantest = reshape(counts(temp, 0:7), 2, 2, 2)
      g2 = 0.0
      for i1 in 1:2, i2 in 1:2, i3 in 1:2
        tt = trantest[i1, i2, i3]
        if tt > 0
          fitted = sum(trantest[:, i2, i3]) * sum(trantest[i1, i2, :]) /
                   sum(trantest[:, i2, :])
          g2 += 2.0 * tt * log(tt / fitted)
        end
      end
      bic = g2 - 2.0 * log(ntest - 2.0)
    end
    tranfinal = counts(test[1:(ntest - 1)] + 2 * test[2:ntest], 0:3)
    alpha = tranfinal[3] / (tranfinal[1] + tranfinal[3])
    beta = tranfinal[2] / (tranfinal[2] + tranfinal[4])
    kthin *= step(range)
    m = log(eps * (alpha + beta) / max(alpha, beta)) /
        log(abs(1.0 - alpha - beta))
    burnin = kthin * ceil(m) + range.start - 1
    n = ((2.0 - alpha - beta) * alpha * beta * phi^2) /
        (r^2 * (alpha + beta)^3)
    keep = kthin * ceil(n)
    total = burnin + keep
  end
  [kthin, burnin, total, nmin, total / nmin]
end

function rafterydiag(c::AbstractChains; q::Real=0.025, r::Real=0.005,
                     s::Real=0.95, eps::Real=0.001)
  _, p, m = size(c.value)
  vals = Array{Float64}(undef, p, 5, m)
  for j in 1:p, k in 1:m
    vals[j, :, k] = rafterydiag(c.value[:, j, k], q=q, r=r, s=s, eps=eps,
                                range=c.range)
  end
  hdr = header(c) * "\nRaftery and Lewis Diagnostic:\n" *
        "Quantile (q) = $q\nAccuracy (r) = $r\nProbability (s) = $s\n"
  ChainSummary(vals, c.names, ["Thinning", "Burn-in", "Total", "Nmin",
                               "Dependence Factor"], hdr)
end


text:
# Time stepping

function init_timestep(param::Param)
   nI, nJ, nK = param.nI, param.nJ, param.nK
   time_G = Array{Float64,3}(undef, nI, nJ, nK)
end

function calc_timestep!(param::Param, speedFlux::SpeedFlux, time_G)

   CFL = param.CFL
   CellSize_D = param.CellSize_D
   Cmax_XF = speedFlux.Cmax_XF
   Cmax_YF = speedFlux.Cmax_YF
   Cmax_ZF = speedFlux.Cmax_ZF

   nI, nJ, nK = param.nI, param.nJ, param.nK

   time_G .= 0.0
   if nI > 1
      @inbounds for k=1:nK, j=1:nJ, i=1:nI
         time_G[i,j,k] += CFL / max(Cmax_XF[i+1,j,k], Cmax_XF[i,j,k]) * CellSize_D[1]
      end
   end

   if nJ > 1
      @inbounds for k=1:nK, j=1:nJ, i=1:nI
         time_G[i,j,k] += CFL / max(Cmax_YF[i,j+1,k], Cmax_YF[i,j,k]) * CellSize_D[2]
      end
   end

   if nK > 1
      @inbounds for k=1:nK, j=1:nJ, i=1:nI
         time_G[i,j,k] += CFL / max(Cmax_ZF[i,j,k+1], Cmax_ZF[i,j,k]) * CellSize_D[3]
      end
   end

   dt = param.TimeAccurate ? minimum(time_G) : 0.0

   return dt
end


text:
# function eval_metric(::Val{:mse}, pred::AbstractMatrix{T}, Y::AbstractVector{T}, α=0.0) where T <: AbstractFloat
#     eval = mean((pred .- Y) .^ 2)
#     return eval
# end
function eval_metric(::Val{:mse}, pred::AbstractMatrix{T}, y::AbstractVector{T}, α) where T <: AbstractFloat
    eval = zero(T)
    @inbounds for i in 1:length(pred)
        eval += (pred[i,1] - y[i]) ^ 2
    end
    eval /= length(pred)
    return eval
end

function eval_metric(::Val{:rmse}, pred::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where T <: AbstractFloat
    eval = zero(T)
    @inbounds for i in 1:length(pred)
        eval += (pred[i,1] - y[i]) ^ 2
    end
    eval = sqrt(eval/length(pred))
    return eval
end

function eval_metric(::Val{:mae}, pred::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where T <: AbstractFloat
    eval = zero(T)
    @inbounds for i in 1:length(pred)
        eval += abs(pred[i,1] - y[i])
    end
    eval /= length(pred)
    return eval
end

function eval_metric(::Val{:logloss}, pred::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where T <: AbstractFloat
    eval = zero(T)
    @inbounds for i in 1:length(y)
        eval -= y[i] * log(max(1e-8, sigmoid(pred[i,1]))) + (1 - y[i]) * log(max(1e-8, 1 - sigmoid(pred[i,1])))
    end
    eval /= length(y)
    return eval
end

function eval_metric(::Val{:mlogloss}, pred::AbstractMatrix{T}, y::Vector{S}, α=0.0) where {T <: AbstractFloat, S <: Integer}
    eval = zero(T)
    L = length(y)
    K = size(pred,2)
    # pred = pred - maximum.(pred)
    @inbounds for i in 1:L
        pred[i] = pred[i] .- maximum(pred[i])
        soft_pred = exp.(pred[i]) / sum(exp.(pred[i]))
        eval -= log(soft_pred[y[i]])
    end
    eval /= length(y)
    return eval
end

function eval_metric(::Val{:poisson}, pred::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where T <: AbstractFloat
    eval = zero(T)
    @inbounds for i in 1:length(y)
        eval += exp(pred[i,1]) * (1 - y[i]) + log(factorial(y[i]))
    end
    eval /= length(y)
    return eval
end

# gaussian
# pred[i][1] = μ
# pred[i][2] = log(σ)
function eval_metric(::Val{:gaussian}, pred::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where {L, T <: AbstractFloat}
    eval = zero(T)
    @inbounds for i in 1:length(y)
        eval += pred[i,2] + (y[i] - pred[i,1])^2 / (2*max(1e-8, exp(2*pred[i,2])))
    end
    eval /= length(y)
    return eval
end

function eval_metric(::Val{:quantile}, pred::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where T <: AbstractFloat
    eval = zero(T)
    for i in 1:length(y)
        eval += α * max(y[i] - pred[i,1], zero(T)) + (1-α) * max(pred[i,1] - y[i], zero(T))
    end
    eval /= length(y)
    return eval
end


text:
# ---
# title: 62. Unique Paths
# id: problem62
# author: AquaIndigo
# date: 2020-12-09
# difficulty: Medium
# categories: Array, Dynamic Programming
# link: <https://leetcode.com/problems/unique-paths/description/>
# hidden: true
# ---
#
# A robot is located at the top-left corner of a `m x n` grid (marked 'Start' in
# the diagram below).
#
# The robot can only move either down or right at any point in time. The robot
# is trying to reach the bottom-right corner of the grid (marked 'Finish' in the
# diagram below).
#
# How many possible unique paths are there?
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2018/10/22/robot_maze.png)
#
#
#
#     Input: m = 3, n = 7
#     Output: 28
#
#
# **Example 2:**
#
#
#
#     Input: m = 3, n = 2
#     Output: 3
#     Explanation:
#     From the top-left corner, there are a total of 3 ways to reach the bottom-right corner:
#     1. Right -> Down -> Down
#     2. Down -> Down -> Right
#     3. Down -> Right -> Down
#
#
# **Example 3:**
#
#
#
#     Input: m = 7, n = 3
#     Output: 28
#
#
# **Example 4:**
#
#
#
#     Input: m = 3, n = 3
#     Output: 6
#
#
#
#
# **Constraints:**
#
#   * `1 <= m, n <= 100`
#   * It's guaranteed that the answer will be less than or equal to `2 * 109`.
#
#
## @lc code=start
using LeetCode

unique_paths(m::Int, n::Int) = binomial(m + n - 2, n - 1)
## @lc code=end