text:
export AbstractCuArray, AbstractCuVector, AbstractCuMatrx, AbstractCuVecOrMat
abstract type AbstractCuArray{T,N} end

const AbstractCuVector{T} = AbstractCuArray{T,1}
const AbstractCuMatrix{T} = AbstractCuArray{T,2}
const AbstractCuVecOrMat{T} = Union{AbstractCuVector{T},AbstractCuMatrix{T}}

Base.similar(x::AbstractCuArray{T}) where T = CuArray{T}(size(x))
Base.similar(x::AbstractCuArray{T}, dims::NTuple) where T = CuArray{T}(dims)
Base.similar(x::AbstractCuArray{T}, dims::Int...) where T = similar(x, dims)
Base.length(x::AbstractCuArray) = prod(size(x))
Base.isempty(x::AbstractCuArray) = length(x) == 0
Base.ndims(x::AbstractCuArray{T,N}) where {T,N} = N
Base.eltype(x::AbstractCuArray{T}) where T = T
Base.zeros(x::AbstractCuArray{T,N}) where {T,N} = zeros(CuArray{T}, size(x))
Base.ones(x::AbstractCuArray{T}) where T = ones(CuArray{T}, size(x))
Base.copy(src::AbstractCuArray) = copyto!(similar(src), src)

function Base.copyto!(dest::AbstractCuArray{T}, src::AbstractCuArray{T}) where T
    @assert length(dest) == length(src)
    if iscontigious(dest) && iscontigious(src)
        if isa(dest, CuSubArray)
            doffs = dest.offset + 1
            dest = dest.parent
        else
            doffs = 1
        end
        if isa(src, CuSubArray)
            soffs = src.offset + 1
            src = src.parent
        else
            soffs = 1
        end
        #if isa(dest, CuSubArray)
        #    dest = unsafe_wrap(CuArray, pointer(dest), size(dest))
        #end
        #if isa(src, CuSubArray)
        #    src = unsafe_wrap(CuArray, pointer(src), size(src))
        #end
        copyto!(dest, doffs, src, soffs, length(dest))
    else
        copyto!(CuDeviceArray(dest), CuDeviceArray(src))
    end
    dest
end

function addto!(dest::AbstractCuArray{T}, src::AbstractCuArray{T}) where T
    @assert length(dest) == length(src)
    if iscontigious(dest) && iscontigious(src)
        p_dest = pointer(dest)
        p_src = pointer(src)
        axpy!(length(dest), T(1), p_src, p_dest)
    else
        addto!(CuDeviceArray(dest), CuDeviceArray(src))
    end
    dest
end


text:
using LightGraphs

@testset "Constructors" begin
   @test @inferred(StorageGraph()) == StorageGraph()
   @test eltype(StorageGraph{UInt8}()) == UInt8
   @test eltype(StorageGraph()) == eltype(SimpleDiGraph())
   @test edgetype(StorageGraph()) == edgetype(SimpleDiGraph())
   @test nv(StorageGraph()) == 0
   @test ne(StorageGraph()) == 0
   @test is_directed(StorageGraph)
   @test is_directed(StorageGraph())
   @test eltype(@inferred(StorageGraph{UInt8}())) == UInt8
   @test @inferred(SimpleDiGraph(StorageGraph())) == SimpleDiGraph()

   g = SimpleDiGraph(3)
   add_edge!(g, 1, 2)
   add_edge!(g, 2, 3)
   data = Dict(1=>(a=1,),2=>(a=2,),3=>(a=3,))
   paths = Dict(Edge(1, 2)=>Set(1), Edge(2, 3)=>Set(1))
   sg = StorageGraph(g, data, paths)
   @test @inferred(StorageGraph(g, data, paths)) == sg
   @test nv(sg) == nv(g)
   @test ne(sg) == ne(g)
   @test sg.data == data
   @test sg.paths == paths
   @test sg.maxid[] == 2
   @test sg.index == Dict((a=1,)=>1,(a=2,)=>2,(a=3,)=>3)

   @test @inferred(StorageGraphs.fadj(sg, 2)) == LightGraphs.SimpleGraphs.fadj(g, 2)
   @test @inferred(StorageGraphs.badj(sg, 2)) == LightGraphs.SimpleGraphs.badj(g, 2)

   @test outneighbors(sg, 1) == outneighbors(sg, (a=1,)) == [2]

   ng = StorageGraph{UInt8}(sg)
   @test eltype(ng) == UInt8
   @test eltype(ng.data) == Pair{UInt8,NamedTuple}
   @test eltype(ng.maxid) == UInt8
end

@testset "Elementary operations" begin
   g = StorageGraph()
   @test add_vertex!(g, (a=1,))
   @test nv(g) == 1
   @test ne(g) == 0
   @test g.data == Dict(1=>(a=1,))
   @test g.index == Dict((a=1,)=>1)
   @test has_prop(g, 1, :a)
   @test sprint(show, g) == "{1, 0} $(eltype(g)) storage graph"

   @test add_vertex!(g, (b=2,))
   @test nv(g) == 2
   @test ne(g) == 0
   @test g.data == Dict(1=>(a=1,), 2=>(b=2,))
   @test g.index == Dict((a=1,)=>1, (b=2,)=>2)
   @test has_prop(g, 2, :b)
   @test !set_prop!(g, 3, (x="doesn't exist",))
   @test !set_prop!(g, 1, 2, 0)

   @test add_edge!(g, 1, 2, 1)
   @test nv(g) == 2
   @test ne(g) == 1
   @test g.paths == Dict(Edge(1,2)=>Set(1))
   @test has_prop(g, 1, 2, 1)

   @test g.maxid[] == get_prop(g) == 1
   set_prop!(g, 2)
   @test get_prop(g) == 2

   @test rem_edge!(g, 1, 2)
   @test ne(g) == 0
   @test !has_prop(g, 1, 2, 1)
   @test rem_vertex!(g, 2)
   @test nv(g) == 1
   @test !has_prop(g, 2, :b)
   @test rem_vertex!(g, 1)
   @test nv(g) == 0
   @test !has_prop(g, 1, :a)
   @test g == StorageGraph()
end

ext: jl


text:
#!/usr/bin/env julia

# julia script to run 1 simulation
# argument: integer that will be mapped to a combination of parameters
# assumption: directory "simresults/" created prior to running the script

# This is to be used with slurm, which defines a shell variable SLURM_ARRAY_TASK_ID.
# This variable will serve the integer argument to the script,
# to determine the combination of parameters that each slurm task will run.

# example use, with default of 2 replicates per parameter combination, combination 14:
# julia onesimulation.jl 14
# or simply:
# ./onesimulation.jl 14
# then look at result file in the result directory

# example use with 200 replicates per parameter combination:
# julia onesimulation.jl 14 200

resultdirectory = "simresults"

"""
    onesimulation(N, n_rare_cat, mu, seed)

Simulate data from a multinomial distribution with 3 categories
and `N` units total.
- `n_rare_cat` = 1 or 2 rare categories, with very small probabilities
- `mu`: expected (or mean) number of observations in each rare category.
  The probability that each observation is in a given rare category is
  `p=mu/N`.

With 2 rare categories, the probabilities are [p,p,1-2p].
With 1 rare category, the probabilities are [p,(1-p)/2,(1-p)/2].

The data (counts in each categories) are then analyzed for a
goodness-of-fit test:
1. first using the Pearson test statistics: sum of (oi - ei)^2 / ei
   where ei = expected count in category i, oi = observed count
2. second using Qlog statistics: 2 * sum of L(ni,ei) * (oi - ei)/ei
   where L(x,y) is the logarithmic mean of x and y:
   L(x,y) = (x-y) / (log(x)-log(y)).
   L(x,y) is always between x and y, as expected for a mean.
   See for example Lorenzen (1995) for a justification for using Qlog.

Both test statistics are compared to a chi-square distribution
with 2 degrees of freedom to obtain p-values.

output: (pearson_stat, Qlog_stat, pearson_pvalue, Qlog_pvalue)

examples:

```julia
onesimulation( 50, 1, 0.1)
onesimulation(500, 2, 2.0)
```
"""
function onesimulation(N, nrarecat, mu, seed)
  Random.seed!(seed)
  p = mu/N
  @assert p>=0.0 && p<=0.5 "the rare category/ies cannot have a negative or >.5 probability"
  if nrarecat==1
    probs = [p,(1.0-p)/2,(1.0-p)/2]
    @assert p<=1.0
  else
    probs = [p, p, (1.0-2p)]
  end
  df = 2 # because 3 categories
  observed = rand(Multinomial(N,probs), 1) # 3x1 array
  phat = observed./N # estimated probabilities of the 3 categories
  pearson = N * sum((phat - probs).^2 ./probs)
  Qlog = 0.0
  for i in 1:3
    if observed[i] == 0 continue; end # no need to do xxx / log(0) = 0.0
    ratio = phat[i]/probs[i]
    ratio ≉ 1.0 || continue
    dev = ratio - 1.0 # deviation on ratio scale
    # if ratio ≈ 1 then dev/log(ratio) ≈ 1, but NaN if ratio = 1 exactly
    Qlog += probs[i] * dev^2 / log(ratio)
  end
  Qlog *= 2N
  # cdf = cumulative distribution function: P(random value <= x)
  # ccdf = complementary cdf:               P(random value >  x)
  p_pearson = ccdf(Chisq(df), pearson)
  p_Qlog    = ccdf(Chisq(df), Qlog)
  return pearson,Qlog, p_pearson,p_Qlog
end

"""
    arrayID_to_parameters(arrayID, total_n_reps)

Convert a single integer into 4 parameters to be used for the simulation:
`sample_size, n_rare_categories, mu, replicate_number`.
The integer is meant to be used as the "array ID" when submitting all these
simulations with slurm and with a slurm array.

Assumption: parameters vary in
- sample_size: 30, 1000
- n_rare_categories: 1, 2
- mu: 0.1, 1.0, 2.0
- replicate_number: from 1 to total_n_reps (input: to start small when developing the script)
"""
function arrayID_to_parameters(arrayID::Int, Nrep::Int)
  # start with Nrep = 2 to develop the julia script + slurm script,
  # then ramp up to Nrep = 200 or more
  samplesizes = [30, 1000]
  nrarecats = [1, 2]
  mus = [0.1, 1., 2.]

  Nssize    = length(samplesizes)
  Nnrarecat = length(nrarecats)
  Nmu       = length(mus)

  ntotal = Nrep * Nssize * Nnrarecat * Nmu
  @assert arrayID <= ntotal "ID must be between 1 and $ntotal"
  I = CartesianIndices( (Nrep,Nssize,Nnrarecat,Nmu) )
  I = I[arrayID].I # coordinates corresponding to the linear index arrayID
  rep = I[1]
  samplesize = samplesizes[I[2]]
  nrarecat = nrarecats[I[3]]
  mu = mus[I[4]]
  return rep, samplesize, nrarecat, mu
end

#---------------------------------------#
#  finally: read script argument        #
#           and run the simulation      #
#---------------------------------------#

# parse the integer argument
@assert length(ARGS)>0 "need 1 parameters: arrayID"
arrayID = parse(Int, ARGS[1])
if length(ARGS)==1
  Nreps = 2
  @info "number of replicates will be 2 by default"
else
  Nreps = parse(Int, ARGS[2])
  @info "number of replicates: $Nreps"
end
rep, samplesize, nrarecat, mu = arrayID_to_parameters(arrayID, Nreps)
@info ""  arrayID replicate=rep samplesize number_rare_categories=nrarecat mu

# load the packages that the simulation function needs
using Distributions # needs to be installed prior
using Random        # no install: standard library
using Printf        # standard library

# run the simulation. use arrayID as seed: will be different for each simulation
pearson, qlog, p_pearson, p_qlog = onesimulation(samplesize, nrarecat, mu, arrayID)

# save the result in tiny csv-formatted file
# (later, all these files will be concatenated with "cat")
outputfile = joinpath(resultdirectory, "simulation_" * @sprintf("%04d", arrayID) * ".csv")
open(outputfile, "w") do g
  write(g, "$arrayID,$rep,$samplesize,$nrarecat,$mu,") # input
  write(g, "$pearson,$qlog,$p_pearson,$p_qlog\n")      # output
end
@info "ID $arrayID is done"

max_stars_repo_name: morrowcj/computingtools

max_issues_repo_name: morrowcj/computingtools

max_forks_repo_name: morrowcj/computingtools


text:
cd(joinpath(Pkg.dir("KernSmooth"), "deps"))
run(`make libkernsmooth.so`)


text:
using DiffEqFlux, OrdinaryDiffEq, Flux, Optim, Test

# General loss function to compare single shooting and multiple shooting predictions
function general_loss_function(result_neuralode)
	return sum(abs2, (ode_data[:,:] .- Array(prob_neuralode(u0, result_neuralode.minimizer)) ))
end

# Define initial conditions and timesteps
datasize = 30
u0 = Float32[2.0, 0.0]
tspan = (0.0f0, 5.0f0)
tsteps = range(tspan[1], tspan[2], length = datasize)

# Get the data
function trueODEfunc(du, u, p, t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end
prob_trueode = ODEProblem(trueODEfunc, u0, tspan)
ode_data = Array(solve(prob_trueode, Tsit5(), saveat = tsteps))

# Define the Neural Network
dudt2 = FastChain((x, p) -> x.^3,
                  FastDense(2, 16, tanh),
                  FastDense(16, 2))
prob_neuralode = NeuralODE(dudt2, (0.0,5.0), Tsit5(), saveat = tsteps)

function loss_neuralode(p)
    pred = Array(prob_neuralode(u0, p))
    loss = sum(abs2, (ode_data[:,1:size(pred,2)] .- pred))
    return loss, pred
end

result_neuralode = DiffEqFlux.sciml_train(loss_neuralode, prob_neuralode.p,
                                          ADAM(0.05),
										  maxiters = 300)

single_shoot_loss = general_loss_function(result_neuralode)
println("single_shoot_loss: ",single_shoot_loss)

# Define parameters for Multiple Shooting
grp_size_param = 1
loss_multiplier_param = 100

neural_ode_f(u,p,t) = dudt2(u,p)
prob_param = ODEProblem(neural_ode_f, u0, tspan, initial_params(dudt2))

function loss_function_param(ode_data, pred):: Float32
	return sum(abs2, (ode_data .- pred))^2
end

function loss_neuralode_param(p)
	return multiple_shoot(p, ode_data, tsteps, prob_param, loss_function_param, grp_size_param, loss_multiplier_param)
end

multiple_shoot_result_neuralode_1 = DiffEqFlux.sciml_train(loss_neuralode_param, prob_neuralode.p,
                                          ADAM(0.05),
                                          maxiters = 300)

multiple_shoot_loss_1 = general_loss_function(multiple_shoot_result_neuralode_1)
println("multiple_shoot_loss_1: ",multiple_shoot_loss_1)

# test for grp_size = 1
@test multiple_shoot_loss_1 < single_shoot_loss

# test for grp_size = 5
grp_size_param = 5
multiple_shoot_result_neuralode_2 = DiffEqFlux.sciml_train(loss_neuralode_param, prob_neuralode.p,
                                          ADAM(0.05),
                                          maxiters = 300)

multiple_shoot_loss_2 = general_loss_function(multiple_shoot_result_neuralode_2)
println("multiple_shoot_loss_2: ",multiple_shoot_loss_2)

@test multiple_shoot_loss_2 < single_shoot_loss


text:
using Gadfly, DataFrames
df = readtable("medical.csv")
print("size is ", size(df)) #> size is (50,7)
df[1:3, 1:size(df,2)]
# data sample:
# IX	Sex	Age	sBP	    dBP	  Drink	BMI
# 0	1	39	106.0	70.0	0	26.97
# 1	2	46	121.0	81.0	0	28.73
# 2	1	48	127.5	80.0	1	25.34

# transforming the data:
df[:Sex] = ifelse(df[:Sex].==1, "F", "M")
df[:Drink] = ifelse(df[:Drink].==1, "Y", "N")
df[1:3, 1:size(df,2)]
# transformed data sample:
# IX	Sex	Age	sBP	    dBP	  Drink	BMI
# 0	F	39	106.0	70.0	N	26.97
# 1	M	46	121.0	81.0	N	28.73
# 2	F	48	127.5	80.0	Y	25.34

describe(df)
# IX
# Min      0.0
# 1st Qu.  12.25
# Median   24.5
# Mean     24.5
# 3rd Qu.  36.75
# Max      49.0
# NAs      0
# NA%      0.0%

# Sex
# Length  50
# Type    ASCIIString
# NAs     0
# NA%     0.0%
# Unique  2

# Age
# Min      35.0
# 1st Qu.  42.0
# Median   46.0
# Mean     47.86
# 3rd Qu.  52.75
# Max      65.0
# NAs      0
# NA%      0.0%

# sBP
# Min      96.0
# 1st Qu.  119.5
# Median   131.5
# Mean     133.5
# 3rd Qu.  145.75
# Max      206.0
# NAs      0
# NA%      0.0%

# dBP
# Min      63.0
# 1st Qu.  76.375
# Median   84.25
# Mean     84.86
# 3rd Qu.  90.75
# Max      121.0
# NAs      0
# NA%      0.0%

# Drink
# Length  50
# Type    ASCIIString
# NAs     0
# NA%     0.0%
# Unique  2

# BMI
# Min      18.59
# 1st Qu.  23.25
# Median   26.15
# Mean     26.374200000000002
# 3rd Qu.  28.472499999999997
# Max      40.11
# NAs      0
# NA%      0.0%

# And the standard deviations
["Age" std(df["Age"]) "sBP" std(df["sBP"]) "dBP" std(df["dBP"]) "BMI" std(df["BMI"])]
# 1x8 Array{Any,2}:
#  "Age"  8.1941  "sBP"  22.6628  "dBP"  12.5921  "BMI"  4.41216

set_default_plot_size(20cm, 12cm)
plot(df, x="Age", y="sBP", color="Sex", Geom.smooth, Geom.bar(position=:dodge))
pl = plot(df, x="Age", y="sBP", color="Sex", Geom.smooth, Geom.bar(position=:dodge))
draw(PDF("medical.pdf", 6inch, 3inch), pl)

max_stars_repo_path: Module 1/Chapter10/medical.jl

max_issues_repo_path: Module 1/Chapter10/Old/medical.jl

max_forks_repo_path: Module 1/Chapter10/Old/medical.jl


text:
@testset "$TEST $G" begin

#### Graphs for testing
graphs = [
  # Graph with 8 vertices
  (8,
   [
     (1,2,10),(1,3,5),(1,4,15),(2,3,4),(2,5,9),
     (2,6,15),(3,4,4),(3,6,8),(4,7,16),(5,6,15),
     (5,8,10),(6,7,15),(6,8,10),(7,3,6),(7,8,10)
   ],
   1,8, # source/target
   3,   # answer for default capacity
   28,  # answer for custom capacity
   15,5 # answer for restricted capacity/restriction
  ),

  # Graph with 6 vertices
  (6,
   [
     (1,2,9),(1,3,9),(2,3,10),(2,4,8),(3,4,1),
     (3,5,3),(5,4,8),(4,6,10),(5,6,7)
   ],
   1,6, # source/target
   2,   # answer for default capacity
   12,  # answer for custom capacity
   8,5  # answer for restricted capacity/restriction
  )
]

for (nvertices,flow_edges,s,t,fdefault,fcustom,frestrict,caprestrict) in graphs
    flow_graph = DG(nvertices)
    capacity_matrix = zeros(Int,nvertices,nvertices)
    for e in flow_edges
        u,v,f = e
        add_edge!(flow_graph,u,v)
        capacity_matrix[u,v] = f
    end

    # Test DefaultCapacity
    d = Erdos.DefaultCapacity(flow_graph)
    @test typeof(d) <: AbstractMatrix{typeof(signed(V(0)))}
    @test d[s,t] == 0
    @test size(d) == (nvertices,nvertices)

    fdef1, Fdef1, labdef1 = maximum_flow(flow_graph,s,t)
    fdef2, Fdef2, labdef2 = maximum_flow(flow_graph,s,t, capacity_matrix)
    # Test all algorithms
    for ALGO in [EdmondsKarpAlgorithm, DinicAlgorithm, BoykovKolmogorovAlgorithm, PushRelabelAlgorithm]
        f, F, lab = maximum_flow(flow_graph,s,t,algorithm=ALGO())
        @test f == fdefault
        @test lab == labdef1

        f, F, lab = maximum_flow(flow_graph,s,t,capacity_matrix,algorithm=ALGO())
        @test f == fcustom
        @test lab == labdef2

        @test maximum_flow(flow_graph,s,t,capacity_matrix,algorithm=ALGO(),restriction=caprestrict)[1] == frestrict
    end
end

n = 30
s = n+1
t = n+2
g = random_regular_digraph(n, 3, DG)
c = spzeros(n+2,n+2)
for e in edges(g)
    c[src(e), dst(e)] = rand()
end
add_vertices!(g, 2)
for i=1:10
    j = rand(1:n)
    add_edge!(g, s, j)
    c[s, j] = rand()

    j = rand(1:n)
    add_edge!(g, j, t)
    c[j, t] = rand()
end

fdef, Fdef, labdef = maximum_flow(g,s,t)
fcustom, Fcustom, labcustom = maximum_flow(g,s,t, c)

enedef, cutdef, labd = minimum_cut(g,s,t)
enecustom, cutcustom, labc = minimum_cut(g,s,t, c)
@test enedef ≈ fdef
@test enecustom ≈ fcustom
@test labd == labdef
@test labc == labcustom
for e in cutcustom
    @test has_edge(g, e)
    @test labc[src(e)] == labc[s]
    @test labc[dst(e)] == labc[t]
end
for e in cutdef
    @test has_edge(g, e)
    @test labd[src(e)] == labd[s]
    @test labd[dst(e)] == labd[t]
end
for ALGO in [EdmondsKarpAlgorithm, DinicAlgorithm, BoykovKolmogorovAlgorithm, PushRelabelAlgorithm]
    f, F, lab = maximum_flow(g,s,t,algorithm=ALGO())
    @test f ≈ fdef
    @test lab == labdef
    f, F, lab = maximum_flow(g,s,t,c,algorithm=ALGO())
    @test f ≈ fcustom
    @test lab == labcustom
end

end # testset


text:
abstract type MatrixTriangulationPartition{D, T} <: AbstractTriangulationPartition{D, T} end
abstract type MatrixTriangulationPartitionFull{D, T} <: AbstractTriangulationPartitionFull{D, T} end

abstract type MutableMatrixTriangulationPartition{D, T} <: AbstractMutableTriangulationPartition{D, T} end
abstract type MutableMatrixTriangulationPartitionFull{D, T} <: AbstractMutableTriangulationPartitionFull{D, T} end

struct MatrixTriangulation{D, T} <: MatrixTriangulationPartition{D, T}
	data::AbstractArray{T, 2}
	simplexindices::DelaunayTriangulation
end

mutable struct MutableMatrixTriangulation{D, T} <: MutableMatrixTriangulationPartition{D, T}
	data::AbstractArray{T, 2}
	simplexindices::DelaunayTriangulation
end

# Full info
struct MatrixTriangulationFull{D, T} <: MatrixTriangulationPartitionFull{D, T}
	data::AbstractArray{T, 2}
	simplexindices::DelaunayTriangulation
    simplices::Vector{Simplex}
    radii::Vector{Float64}
    orientations::Vector{Float64}
    volumes::Vector{Float64}
    centroids::Vector{Vector{Float64}}
end

# Full info
struct MutableMatrixTriangulationFull{D, T} <: MutableMatrixTriangulationPartitionFull{D, T}
	data::AbstractArray{T, 2}
	simplexindices::DelaunayTriangulation
    simplices::Vector{Simplex}
    radii::Vector{Float64}
    orientations::Vector{Float64}
    volumes::Vector{Float64}
    centroids::Vector{Vector{Float64}}
end

MT = Union{MatrixTriangulationPartition, MatrixTriangulationPartitionFull,
    MutableMatrixTriangulationPartition, MutableMatrixTriangulationPartitionFull}
getpoints(mt::MT) = [mt.data[:, i] for i = 1:maximum(size(mt.data))]

getsimplices(mt::Union{MatrixTriangulationPartitionFull, MutableMatrixTriangulationPartitionFull}) =
    mt.simplices
getsimplices(mt::Union{MatrixTriangulationPartition, MutableMatrixTriangulationPartition}) =
    [Simplex(mt.data[:, mt.simplexindices[i]]) for i = 1:length(mt.simplexindices)]

####################
# Pretty printing
####################
function summarise(tp::Union{MatrixTriangulation, MutableMatrixTriangulation})
    _type = typeof(tp)
    n_simplices = nsimplices(tp)
    n_pts = maximum(size(tp.data))
    D = dimension(tp)
    summary = "$_type from $n_simplices simplices and $n_pts points"
end

Base.show(io::IO, tp::Union{MatrixTriangulation, MutableMatrixTriangulation}) =
    println(io, summarise(tp))

export
MatrixTriangulation,
MutableMatrixTriangulation,
MatrixTriangulationFull,
MutableMatrixTriangulationFull,
getpoints


text:
#=
distributions.jl

25/08/2015
Adham Beyki, odinay@gmail.com

TODO
 1. revise pdf and logpdf of Dirichlet, they are in fact expectations at the moment
 2. add check_args to inputs
=#

abstract Distribution

####################################
############ Gaussian1D ############
####################################
immutable Gaussian1D <: Distribution
    mu::Float64
    vv::Float64

    Gaussian1D(mu::Float64, vv::Float64) = new(mu, vv)
end
Gaussian1D(mu::Real, vv::Real) = Gaussian1D(convert(Float64, mu), convert(Float64, vv))

Base.deepcopy(me::Gaussian1D)  = Gaussian1D(me.mu, me.vv)

# displaying the Gaussian1D distribution
function Base.show(io::IO, me::Gaussian1D)
    println(io, "Gaussian1D distribution")
    println(io, "mean: mu=", me.mu, ", variance: vv=", me.vv)
end

# sampling
sample(me::Gaussian1D)    = me.mu + sqrt(me.vv)*randn()
sample(me::Gaussian1D, n) = me.mu + sqrt(me.vv)*randn(n)

# pdf
pdf(me::Gaussian1D, x::Float64)    = exp(-(x-me.mu)^2/(2*me.vv)) / (sqrt(2*pi*me.vv))
logpdf(me::Gaussian1D, x::Float64) = log(pdf(me, x))

##############################
######### Dirichlet ##########
##############################
type Dirichlet <: Distribution
    alpha::Vector{Float64}

    Dirichlet(alpha::Vector{Float64}) = new(alpha)
end

Dirichlet(alpha::Vector{Int})    = Dirichlet(convert(Vector{Float64}, alpha))
Dirichlet(alpha::Real, dim::Int) = Dirichlet(fill(alpha/dim, dim))
Base.deepcopy(me::Dirichlet)     = Dirichlet(deepcopy(me.alpha))

# displaying the Dirichlet distribution
function Base.show(io::IO, me::Dirichlet)
    println(io, "Dirichlet distribution")
    println(io, "cardinality = ", length(me.alpha))
    println(io, "alpha       = ", me.alpha)
end

Base.mean(me::Dirichlet) = me.alpha / sum(me.alpha)

pdf(me::Dirichlet, x::Int)    = mean(me)[x]
logpdf(me::Dirichlet, x::Int) = log(pdf(me, x))


text:
@inline function fit_cubic(; x0, xf, y0, yd0, yf, ydf)
    num_coeffs = Val(4)
    A = hcat(
        SVector(SUP.coefficient_gradient(x0, num_coeffs, Val(0))),
        SVector(SUP.coefficient_gradient(x0, num_coeffs, Val(1))),
        SVector(SUP.coefficient_gradient(xf, num_coeffs, Val(0))),
        SVector(SUP.coefficient_gradient(xf, num_coeffs, Val(1))),
    ) |> transpose
    b = SVector(tuple(y0, yd0, yf, ydf))
    Polynomial(Tuple(A \ b))
end

@inline function fit_quintic(; x0, xf, y0, yd0, ydd0, yf, ydf, yddf)
    num_coeffs = Val(6)
    A = hcat(
        SVector(SUP.coefficient_gradient(x0, num_coeffs, Val(0))),
        SVector(SUP.coefficient_gradient(x0, num_coeffs, Val(1))),
        SVector(SUP.coefficient_gradient(x0, num_coeffs, Val(2))),
        SVector(SUP.coefficient_gradient(xf, num_coeffs, Val(0))),
        SVector(SUP.coefficient_gradient(xf, num_coeffs, Val(1))),
        SVector(SUP.coefficient_gradient(xf, num_coeffs, Val(2))),
    ) |> transpose
    b = SVector(tuple(y0, yd0, ydd0, yf, ydf, yddf))
    Polynomial(Tuple(A \ b))
end


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

mutable struct IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    fieldsType::Any # spec type: Union{ Nothing, String } # spec name: fieldsType
    fieldsV1::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1FieldsV1 } # spec name: fieldsV1
    manager::Any # spec type: Union{ Nothing, String } # spec name: manager
    operation::Any # spec type: Union{ Nothing, String } # spec name: operation
    time::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: time

    function IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry(;apiVersion=nothing, fieldsType=nothing, fieldsV1=nothing, manager=nothing, operation=nothing, time=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("fieldsType"), fieldsType)
        setfield!(o, Symbol("fieldsType"), fieldsType)
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("fieldsV1"), fieldsV1)
        setfield!(o, Symbol("fieldsV1"), fieldsV1)
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("manager"), manager)
        setfield!(o, Symbol("manager"), manager)
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("operation"), operation)
        setfield!(o, Symbol("operation"), operation)
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("time"), time)
        setfield!(o, Symbol("time"), time)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry

const _property_map_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("fieldsType")=>Symbol("fieldsType"), Symbol("fieldsV1")=>Symbol("fieldsV1"), Symbol("manager")=>Symbol("manager"), Symbol("operation")=>Symbol("operation"), Symbol("time")=>Symbol("time"))
const _property_types_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("fieldsType")=>"String", Symbol("fieldsV1")=>"IoK8sApimachineryPkgApisMetaV1FieldsV1", Symbol("manager")=>"String", Symbol("operation")=>"String", Symbol("time")=>"IoK8sApimachineryPkgApisMetaV1Time")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry)
    true
end

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


text:
# integration with LinearAlgebra.jl

CuMatOrAdj{T} = Union{CuMatrix, LinearAlgebra.Adjoint{T, <:CuMatrix{T}}, LinearAlgebra.Transpose{T, <:CuMatrix{T}}}
CuOrAdj{T} = Union{CuVecOrMat, LinearAlgebra.Adjoint{T, <:CuVecOrMat{T}}, LinearAlgebra.Transpose{T, <:CuVecOrMat{T}}}

# matrix division

function Base.:\(_A::CuMatOrAdj, _B::CuOrAdj)
    A, B = copy(_A), copy(_B)
    A, ipiv = CUSOLVER.getrf!(A)
    return CUSOLVER.getrs!('N', A, ipiv, B)
end

# patch JuliaLang/julia#40899 to create a CuArray
# (see https://github.com/JuliaLang/julia/pull/41331#issuecomment-868374522)
if VERSION >= v"1.7-"
_zeros(::Type{T}, b::AbstractVector, n::Integer) where {T} = CUDA.zeros(T, max(length(b), n))
_zeros(::Type{T}, B::AbstractMatrix, n::Integer) where {T} = CUDA.zeros(T, max(size(B, 1), n), size(B, 2))
function Base.:\(F::Union{LinearAlgebra.LAPACKFactorizations{<:Any,<:CuArray},
                          Adjoint{<:Any,<:LinearAlgebra.LAPACKFactorizations{<:Any,<:CuArray}}},
                 B::AbstractVecOrMat)
    m, n = size(F)
    if m != size(B, 1)
        throw(DimensionMismatch("arguments must have the same number of rows"))
    end

    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
    FF = Factorization{TFB}(F)

    # For wide problem we (often) compute a minimum norm solution. The solution
    # is larger than the right hand side so we use size(F, 2).
    BB = _zeros(TFB, B, n)

    if n > size(B, 1)
        # Underdetermined
        copyto!(view(BB, 1:m, :), B)
    else
        copyto!(BB, B)
    end

    ldiv!(FF, BB)

    # For tall problems, we compute a least squares solution so only part
    # of the rhs should be returned from \ while ldiv! uses (and returns)
    # the complete rhs
    return LinearAlgebra._cut_B(BB, 1:n)
end
end

# qr

using LinearAlgebra: AbstractQ

# AbstractQ's `size` is the size of the full matrix,
# while `Matrix(Q)` only gives the compact Q.
# See JuliaLang/julia#26591 and JuliaGPU/CUDA.jl#969.
CuMatrix{T}(Q::AbstractQ{S}) where {T,S} = convert(CuArray, Matrix{T}(Q))
CuMatrix(Q::AbstractQ{T}) where {T} = CuMatrix{T}(Q)
CuArray{T}(Q::AbstractQ) where {T} = CuMatrix{T}(Q)
CuArray(Q::AbstractQ) = CuMatrix(Q)


text:
function multiplyNumbers(x, y)
    return x * y
end

z = multiplyNumbers(1 + 2im, 1 - 3im)
println(z)
println(z^2)


text:
include("load.jl")
TapirSchedulersTests.runtests()


text:
"""
This function selects a single implementation from `fnames` and associates it
with `fname`. It benchmarks every function in `fnames` against the single
type `Vector{T}` of length `N` and writes the resulting association into the
IOstream `file`.

    fname: name to be associated with Vectorize.(fname)
    fnames: a list of strings "Vectorize.VML.sin" to be benchmarked
    T: the type of the Vector{T} to pass to the function
    file: the IOStream to write the final string onto
    N: the length of the vector to benchmark over.
"""
function benchmarkSingleArgFunction(fname, fnames,
                                    T, file::IOStream, N::Integer)

    # dictionary to store elpased values
    val = Dict()
    println("TESTING: $(fname)(Array{$T})")

    # Iterate over all functions
    for fstr in fnames
        f = eval(Meta.parse(fstr)) # convert string to function
        X = randn(eltype(T), N) # function arguments
        f(X) # force compilation
        time = 0
        for i in 1:10 # benchmark ten times
            time += @elapsed f(X)
        end
        val[fstr] = time/10.0
    end

    # default to first function
    t = val[fnames[1]]
    fbest = fnames[1]

    # find fastest function
    for f in fnames
        if val[f] < t
            t = val[f]
            fbest = f
        end
    end

    # write chosen file into IOStream
    write(file, "\n$(fname)(X::$T) = $(fbest)(X)\n")
    # docs = @doc fbest
    # write(file, "\n\"\"\"\n", docs, "\n\"\"\"\n")
    # write(file, "$(fname)(X::Array{$T}) = $(fbest)(X)\n")
    println("BENCHMARK: $(fname)($T) mapped to $(fbest)()\n")
end

"""
This function selects a single implementation from `fnames` and associates it
with `fname`. It benchmarks every function in `fnames` against two `Vector{T}`'s
of length `N` and writes the resulting association into the
IOstream `file`.

    fname: name to be associated with Vectorize.(fname)
    fnames: a list of strings "Vectorize.VML.sin" to be benchmarked
    T: tuple of types of Vector{T} to pass to the function
    file: the IOStream to write the final string onto
    N: the length of the vector to benchmark over.
"""
function benchmarkTwoArgFunction(fname, fnames,
                                 T, file::IOStream, N::Integer)

    # dictionary of elapsed times for each implementation
    val = Dict()
    println("TESTING: $(fname)($(T[1]), $(T[2]))")

    # Iterate over all functions
    for fstr in fnames
        f = eval(Meta.parse(fstr)) # convert string to function
        X = randn(eltype(T[1]), N) # function arguments
        Y = randn(eltype(T[2]), N) # function arguments
        f(X, Y) # force compilation
        time = 0
        for i in 1:10 # benchmark ten times
            time += @elapsed f(X, Y)
        end
        val[fstr] = time/10.0
    end

    # default to first function
    t = val[fnames[1]]
    fbest = fnames[1]

    # find fastest function
    for f in fnames
        if val[f] < t
            t = val[f]
            fbest = f
        end
    end

    # write result into IOStream
    write(file, "\n$(fname)(X::$(T[1]), Y::$(T[2])) = $(fbest)(X, Y)\n")
    println("BENCHMARK: $(fname)($(T[1]), $(T[2])) mapped to $(fbest)()\n")
end

"""
This function selects a single implementation from `fnames` and associates it
with `fname`. It benchmarks every function in `fnames` against three `Array{T}`'s
of length `N` and writes the resulting association into the
IOstream `file`.

    fname: name to be associated with Vectorize.(fname)
    fnames: a list of strings "Vectorize.VML.sin" to be benchmarked
    T: tuple of types of Array{T} to pass to the function
    file: the IOStream to write the final string onto
    N: the length of the vector to benchmark over.
"""
function benchmarkThreeArgFunction(fname, fnames,
                                   T, file::IOStream, N::Integer)

    # dictionary to store elapsed values
    val = Dict()
    println("TESTING: $(fname)($(T[1]), $(T[2]), $(T[3]))")

    # Iterate over all functions
    for fstr in fnames
        f = eval(Meta.parse(fstr)) # convert string to function
        X = randn(eltype(T[1]), N) # function arguments
        Y = randn(eltype(T[2]), N) # function arguments
        Z = randn(eltype(T[3]), N) # function arguments
        f(X, Y, Z) # force compilation
        time = 0
        for i in 1:10 # benchmark 10 times
            time += @elapsed f(X, Y, Z)
        end
        val[fstr] = time/10.0
    end

    # default to first function
    t = val[fnames[1]]
    fbest = fnames[1]

    # find fastest function
    for f in fnames
        if val[f] < t
            t = val[f]
            fbest = f
        end
    end

        # write result into IOStream
    write(file, "\n$(fname)(X::$(T[1]), Y::$(T[2]), Z::$(T[3])) = $(fbest)(X, Y, Z)\n")
    println("BENCHMARK: $(fname)($(T[1]), $(T[2]), $(T[3])) mapped to $(fbest)()\n")
end

"""
`trycmd(cmd::Cmd, msg::String="", err::String="")::String`
This function attemps to run the shell command specified by `cmd` using `run`.
If `cmd` returns successfully, it will print `msg` using info().
If `run(cmd)` throws an exception (the shell command returns non-zero),
then the function throws an error and prints the `err` argument.

Note: This function does not return the value of the command to the caller. However,
prints to STD_OUT by the `run` command will print automatically from the shell process.
"""
function trycmd(cmd::Cmd; msg::String="", err::String="")
    try
        run(cmd)
        if msg == "" return else info(msg) end
    catch
        if err != ""
            error(err)
        end
    end
end

"""
`trycmd_read(cmd::Cmd, msg::String="", err::String="")::String`
This function attemps to run the shell command specified by `cmd` using `readall`.
If `cmd` returns successfully, it will return its result to the caller and print `msg`.
If `readall(cmd)`/`readstring(cmd)` throws an exception (the shell command returns non-zero),
then the function throws an error and prints the `err` argument.
"""
function trycmd_read(cmd::Cmd; msg::String="", err::String="")
    try
        result = readstring(cmd) # Julia v0.5 and above
        if msg == "" return result else info(msg) end
        return result
    catch
        if err != ""
            error(err)
        end
    end
end

"""
`prompt_yn(prompt::String)::Bool`

Prompts the user for a y or n response (y/n). Continually prompts
until user enters a valid response. Prompt string should be provided without
colon and whitespace.
"""
function prompt_yn(prompt::String)
    print(prompt, " [y/n]:  ")
    input = chomp(readline())
    while input != "y" && input != "n"
        print("You did not enter y or n. Please enter a valid response: [y/n]  ")
        input = chomp(readline())
    end

    return (input == "y") ? true : false
end


text:
#
# Compiler state.
#

#
# A combinator constructor.
#

abstract AbstractBinding

#
# Collection of bindings and other compiler state.
#

typealias BindingTuple Tuple{Vararg{AbstractBinding}}
typealias NullableBindingTuple Nullable{BindingTuple}
typealias NullableSymbol Nullable{Symbol}
typealias SymbolTable ImmutableDict{Tuple{Symbol,Int}, AbstractBinding}

immutable Scope
    db::AbstractDatabase
    domain::Type
    rev::Bool
    tag::NullableSymbol
    globals::SymbolTable
    locals::SymbolTable
    items::NullableBindingTuple
end

Scope(
    db::AbstractDatabase, domain::Type;
    rev::Bool=false,
    tag::Union{Symbol,NullableSymbol}=NullableSymbol(),
    globals::SymbolTable=SymbolTable(),
    locals::SymbolTable=SymbolTable(),
    items::Union{BindingTuple,NullableBindingTuple}=NullableBindingTuple()) =
    Scope(db, domain, rev, tag, globals, locals, items)

domain(scope::Scope) = scope.domain

addglobal(scope::Scope, key::Tuple{Symbol,Int}, binding::AbstractBinding) =
    Scope(scope.db, scope.domain, scope.rev, scope.tag, assoc(scope.globals, key, binding), scope.locals, scope.items)

addglobal(scope::Scope, key::Symbol, binding::AbstractBinding) =
    addglobal(scope, (key, 0), binding)

addlocal(scope::Scope, key::Tuple{Symbol,Int}, binding::AbstractBinding) =
    Scope(scope.db, scope.domain, scope.rev, scope.tag, scope.globals, assoc(scope.locals, key, binding), scope.items)

addlocal(scope::Scope, key::Symbol, binding::AbstractBinding) =
    addlocal(scope, (key, 0), binding)

settag(scope::Scope, tag::Symbol) =
    Scope(scope.db, scope.domain, scope.rev, tag, scope.globals, scope.locals, scope.items)

setrev(scope::Scope, rev::Bool) =
    Scope(scope.db, scope.domain, rev, scope.tag, scope.globals, scope.locals, scope.items)

setitems(scope::Scope, items::BindingTuple) =
    Scope(scope.db, scope.domain, scope.rev, scope.tag, scope.globals, scope.locals, items)

lookup(scope::Scope, key::Tuple{Symbol,Int}) =
    begin
        name, arity = key
        b = get(scope.locals, (name, arity))
        if !isnull(b)
            return b
        end
        b = get(scope.locals, (name, -1))
        if !isnull(b)
            return b
        end
        b = get(scope.globals, (name, arity))
        if !isnull(b)
            return b
        end
        return get(scope.globals, (name, -1))
    end
lookup(scope::Scope, key::Symbol) = lookup(scope, (key, 0))

nest(scope::Scope, domain::Type) = nest(scope.db, scope, domain)

#
# Encapsulates the compiler state and the combinator.
#

immutable Query
    scope::Scope
    pipe::AbstractPipe
    syntax::Nullable{AbstractSyntax}
    origin::Nullable{Query}
end

Query(
    scope::Scope,
    pipe::AbstractPipe;
    syntax::Union{AbstractSyntax,Nullable{AbstractSyntax}}=Nullable{AbstractSyntax}(),
    origin::Union{Query,Nullable{Query}}=Nullable{Query}()) =
    Query(scope, pipe, syntax, origin)

Query(
    q::Query;
    scope=nothing,
    pipe=nothing,
    syntax=nothing,
    origin=nothing) =
    Query(
        scope == nothing ? q.scope : scope,
        pipe == nothing ? q.pipe : pipe,
        syntax == nothing ? q.syntax : syntax,
        origin == nothing ? q.origin : origin)

# Initial compiler state.
Query(db::AbstractDatabase; params...) =
    let scope = scope(db)
        for (tag, T) in params
            output =
                T <: Vector ? Output(eltype(T), lunique=false, ltotal=false) :
                T <: Nullable ? Output(eltype(T), ltotal=false) : Output(T)
            scope = addglobal(scope, tag, ParamBinding(tag, output))
        end
        Query(scope, HerePipe(domain(scope)))
    end

# Extracts the compiler state.
scope(q::Query) = q.scope

lookup(q::Query, key) =
    lookup(scope(q), key)

addlocal(scope::Scope, key, q::Query) =
    addlocal(scope, key, SimpleBinding(q))

addglobal(scope::Scope, key, q::Query) =
    addglobal(scope, key, SimpleBinding(q))

# Extracts the combinator.
pipe(q::Query) = q.pipe
convert(::Type{AbstractPipe}, q::Query) = q.pipe

# The input type and structure.
input(q::Query) = input(q.pipe)
# The output type and structure.
output(q::Query) = output(q.pipe)

# Displays the query.
function show(io::IO, q::Query)
    print(io, isnull(q.syntax) ? "(?)" : get(q.syntax), " :: ", mapping(q))
end

# Compiles the query.
prepare(base, expr; params...) = prepare(Query(base; params...), syntax(expr))
prepare(base::Query, expr::AbstractSyntax) =
    if !isnull(base.origin)
        return prepare(get(base.origin), expr)
    else
        origin = compile(base, expr)
        return Query(optimize(format(base.scope, origin)), origin=origin)
    end

format(base::Scope, q::Query) = select(q)

# Executes the query.
execute(q::Query, args...; params...) =
    execute(pipe(optimize(select(q))), args...; params...)
call(q::Query, args...; params...) = execute(q, args...; params...)

# Builds initial execution pipeline.
#compile(base::Query, expr::AbstractSyntax) =
#    error("compile() is not implemented for $(typeof(expr))")

# Optimizes the execution pipeline.
optimize(q::Query) = q
#optimize(q::Query) = Query(q, pipe=optimize(q.pipe))

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
module SymbolicWedderburn

using LinearAlgebra
using Primes

using PermutationGroups
using Cyclotomics

include("gf.jl")
include("eigenspacedecomposition.jl")
include("cmmatrix.jl")

include("characters.jl")
include("characters_arith.jl")
include("powermap.jl")
include("dixon.jl")

include("actions.jl")
include("projections.jl")

end # module


text:
# This file was generated, do not modify it. # hide
#hideall
names = (:Taimur, :Catherine, :Maria, :Arvind, :Jose, :Minjie)
numbers = (1525, 5134, 4214, 9019, 8918, 5757)
println("@@simple-table")
println("Name | Number")
println(":--- | :---")
println.("$name | $number" for (name, number) in zip(names, numbers))
println("@@")
raw"""
~~~
<style>
.simple-table tr {
  padding:0;
  line-height:1em;
}
</style>
~~~
""" |> println

max_stars_repo_name: mcognetta/franklindocs

max_issues_repo_name: mcognetta/franklindocs

max_forks_repo_name: mcognetta/franklindocs


text:
#=
  Implementation de la linearisation de la norme euclidienne (dans IR^2)
  de Camino et al. (2019) applique au "continuous k-center problem"
  X. Gandibleux - Novembre 2019
=#

using JuMP, GLPK
using Printf
using PyPlot

function kCenterPb( XY, K=1, NDIR=28, XMIN=0, XMAX=250, YMIN=0, YMAX=250)

    # XY   : coordonnees (x,y) des lieux a visiter (exemple : XY = [78 97; 61 130; 74 139])
    # K    : nombre de centres a positionner (par defaut 1)
    # NDIR : nombre de directions a calculer dans l'approximation de L2 (par defaut 28)
    # XMIN, XMAX, YMIN, YMAX : # bornes sur X et Y de la boite consideree (par defaut [0,0] x [250,250])

    # ==========================================================================
    # Parametres de l'approximation lineaire de L2  (Camino et al. 2019)

    ΘMAX = π/NDIR
    Ux = zeros(NDIR)
    Uy = zeros(NDIR)
    for i=1:NDIR
        Ux[i] = cos( 2.0 * float(i-1) * π / float(NDIR) )
        Uy[i] = sin( 2.0 * float(i-1) * π / float(NDIR) )
    end

    # ==========================================================================
    # Calcul des majorants sur les lieux a considerer

    function dist(x1,y1,x2,y2)
        d = sqrt( (x1-x2)^2 + (y1-y2)^2)
        return d
    end

    N = size(XY,1)  # N : nombre de lieux a visiter
    M = zeros(N)    # M : vecteur des majorants
    for c=1:N
        M[c] = max( dist(XY[c,1],XY[c,2],XMIN,YMIN) ,
                    dist(XY[c,1],XY[c,2],XMIN,YMAX) ,
                    dist(XY[c,1],XY[c,2],XMAX,YMIN) ,
                    dist(XY[c,1],XY[c,2],XMAX,YMAX)   )
    end

    # ==========================================================================
    # Modele d'optimisation (Camino et al. 2019)
    kcenter = Model(with_optimizer(GLPK.Optimizer))

    @variable(kcenter, λ ≥ 0)
    @variable(kcenter, x[1:K] ≥ 0)
    @variable(kcenter, y[1:K] ≥ 0)
    @variable(kcenter, α[1:N,1:K], Bin)

    @objective(kcenter, Min, λ )

    # contraintes donnant les bornes de la boite
    @constraint(kcenter, limxmin[k=1:K], x[k] ≥ XMIN)
    @constraint(kcenter, limxmax[k=1:K], x[k] ≤ XMAX)
    @constraint(kcenter, limymin[k=1:K], y[k] ≥ YMIN)
    @constraint(kcenter, limymax[k=1:K], y[k] ≤ YMAX)

    # contraintes d'affectation lieu-centre ------------------------------------
    @constraint(kcenter, affectation[c=1:N], sum(α[c,k] for k in 1:K) == 1)

    # contraintes reformulant l'expression non lineaire
    for c=1:N # tous les lieux
        for k=1:K # tous les centres
            for i=1:NDIR # toutes les directions
                @constraint(kcenter, cos(ΘMAX) * λ ≥ ( ((x[k]-XY[c,1])*Ux[i])+((y[k]-XY[c,2])*Uy[i]) ) - (M[c]*(1-α[c,k])) )
            end
        end
    end

    # optimisation et recuperation des resultats
    optimize!(kcenter)
    zOpt = objective_value(kcenter)
    xOpt = value.(x)
    yOpt = value.(y)

    # ==========================================================================
    # Edition des resultats au terminal
    @printf("Distance λOpt=%8.2f \n",zOpt)
    for k=1:K
        print("Centre ", k, " : ")
        @printf(" x=%8.2f | y=%8.2f \n",xOpt[k],yOpt[k])
    end

    # ==========================================================================
    # Edition des resultats graphiquement
    figure("K-Center",figsize=(6,6)) # Create a new figure
    title("\$K=\$" * string(K) * " | \$N_{directions}=\$" * string(NDIR))
    xlabel("X")
    ylabel("Y")
    xlim(XMIN, XMAX)
    ylim(YMIN, YMAX)
    for i in 1:N
        scatter( XY[i,1] , XY[i,2] , s=15, color="blue")
    end
    for k in 1:K
        scatter( xOpt[k] , yOpt[k] , s=15, color="red")
    end
    theta = range(0, stop=2π, length=100)
    r=zOpt
    for k in 1:K
        x1 = r*cos.(theta) .+ xOpt[k]
        x2 = r*sin.(theta) .+ yOpt[k]
        scatter( x1,x2 , s=0.1, color="red")
    end
end

# ==============================================================================
# ==============================================================================

println(" ")
println("USAGE : kCenterPb( XY, K=1, NDIR=28, XMIN=0, XMAX=250, YMIN=0, YMAX=250)")

XY=[218 55; 78 97; 61 130; 74 139; 186 159]
@show XY
kCenterPb(XY, 2, 28, 0, 250, 0, 250)

max_stars_repo_name: xgandibleux/kcenter

max_issues_repo_name: xgandibleux/kcenter

max_forks_repo_name: xgandibleux/kcenter


text:
using Distributed
addprocs(3)

using Test
using Dagger

include("util.jl")
include("fakeproc.jl")

include("thunk.jl")

#= FIXME: Unreliable, and some thunks still get retained
# N.B. We need a few of these probably because of incremental WeakRef GC
@everywhere GC.gc()
@everywhere GC.gc()
@everywhere GC.gc()
@everywhere GC.gc()
sleep(1)
@test isempty(Dagger.Sch.EAGER_ID_MAP)
state = Dagger.Sch.EAGER_STATE[]
@test isempty(state.waiting)
@test_broken length(keys(state.waiting_data)) == 1
# Ensure that all cache entries have expired
@test_broken isempty(state.cache)
=#

include("scheduler.jl")
include("processors.jl")
include("ui.jl")
include("checkpoint.jl")
include("domain.jl")
include("array.jl")
try # TODO: Fault tolerance is sometimes unreliable
include("fault-tolerance.jl")
catch
end
println(stderr, "tests done. cleaning up...")
Dagger.cleanup()
#include("cache.jl")
println(stderr, "all done.")


text:
const ∫ = integrate

doc"""
Integrate a Taylor1Model.
`x0` is optional constant to add.
"""
function integrate(f::Taylor1Model, x0=0)

    p2 = integrate(f.p)

    Δ = integral_bound(f)

    g = Taylor1Model(f, [0; p2[1:end]], Δ)

    return g

end

function integral_bound(f::Taylor1Model)
    n = f.n
    high_order_term = f.p[n]

    coeff = bound(high_order_term)
    power = (f.I - f.x0)^n

    return ((coeff * power) + f.Δ) * diam(f.I)
end

"""
fs is an array of functions
This works in 2 variables currently.
"""
function Taylor_step(fs, n, u0, v0, t_interval)

    u = u0
    v = v0
    u_new = u0   # so exist outside loop
    v_new = v0

    # build up Taylor series by Picard:
    for i in 1:n+1   # how many iterations are required?
        u_new = u0 + ∫(fs[1](u, v))
        v_new = v0 + ∫(fs[2](u, v))

        u, v = u_new, v_new
    end

    # prepare Taylor Model:
    uu = Taylor1Model(n, Interval(t_interval.lo), t_interval, u, 0..0)
    vv = Taylor1Model(n, Interval(t_interval.lo), t_interval, v, 0..0)

    uΔ = integral_bound(fs[1](uu, vv))
    vΔ = integral_bound(fs[2](uu, vv))

    # make sure the intervals contain 0:
    m = mag(uΔ)
    uΔ = -m..m

    m = mag(vΔ)
    vΔ = -m..m

    @show uΔ, vΔ

    uu = Taylor1Model(n, Interval(t_interval.lo), t_interval, u, uΔ)
    vv = Taylor1Model(n, Interval(t_interval.lo), t_interval, v, vΔ)

    integ_u_bound = uΔ
    integ_v_bound = vΔ

    while ! ((integ_u_bound ⊆ uu.Δ) && (integ_v_bound ⊆ vv.Δ))
        uΔ *= 2
        vΔ *= 2

        @show uΔ, vΔ

        uu = Taylor1Model(n, Interval(t_interval.lo), t_interval, u, uΔ)
        vv = Taylor1Model(n, Interval(t_interval.lo), t_interval, v, vΔ)

        integ_u_bound = integral_bound(fs[1](uu, vv))
        integ_v_bound = integral_bound(fs[2](uu, vv))

    end

    # only need to bound the integral, not actually carry out
    # the whole integral

    # contract:

    uu_new = ∫( fs[1](uu, vv), u0[0] )
    vv_new = ∫( fs[2](uu, vv), v0[0] )

    for i in 1:10
       uu, vv = uu_new, vv_new
       uu_new = ∫( fs[1](uu, vv), u0[0] )
       vv_new = ∫( fs[2](uu, vv), v0[0] )
    end

    @show uu_new.Δ, vv_new.Δ

    U = uu_new(t_interval.hi)
    V = vv_new(t_interval.hi)

    return U, V, uu_new, vv_new
end


text:
# Routines to classify interval matrices

"""
    is_strongly_regular(A::AbstractMatrix{T}) where {T<:Interval}

Tests whether the square interval matrix ``A`` is strongly regular, i.e. if ``A_c^{-1}A`` is
an H-matrix, where ``A_c`` is the midpoint matrix of ``A```.
For more details see section 4.6 of [[HOR19]](@ref).

### Examples

```jldoctest
julia> A = [2..4 -2..1; -1..2 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-2, 1]
 [-1, 2]   [2, 4]

julia> is_strongly_regular(A)
true

julia> A = [0..2 1..1;-1.. -1 0..2]
2×2 Matrix{Interval{Float64}}:
   [0, 2]  [1, 1]
 [-1, -1]  [0, 2]

julia> is_strongly_regular(A)
false
```
"""
function is_strongly_regular(A::AbstractMatrix{T}) where {T<:Interval}
    m, n = size(A)
    m == n || return false

    Ac = mid.(A)
    rank(Ac) == n || return false
    return is_H_matrix(inv(Ac)*A)
end

"""
    is_H_matrix(A::AbstractMatrix{T}) where {T<:Interval}

Tests whether the square interval matrix A is an H-matrix, by testing that ``⟨A⟩^{-1}e>0``,
where ``e=[1, 1, …, 1]ᵀ``. Note that in practice it tests that a
_floating point approximation_ of ``⟨A⟩^{-1}e`` satisfies the condition.
For more details see section 4.4 of [[HOR19]](@ref).

### Examples

```jldoctest
julia> A = [2..4 -1..1; -1..1 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-1, 1]
 [-1, 1]   [2, 4]

julia> is_H_matrix(A)
true

julia> A = [2..4 -2..1; -1..2 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-2, 1]
 [-1, 2]   [2, 4]

julia> is_H_matrix(A)
false
```
"""
function is_H_matrix(A::AbstractMatrix{T}) where {T<:Interval}
    m, n = size(A)
    m == n || return false

    compA = comparison_matrix(A)
    rank(compA) == n || return false
    return all(compA\ones(n) .> 0)
end

"""
    is_strictly_diagonally_dominant(A::AbstractMatrix{T}) where {T<:Interval}

Checks whether the square interval matrix ``A`` of order ``n`` is stictly diagonally
dominant, that is if ``mig(Aᵢᵢ) > ∑_{k ≠ i} mag(Aᵢₖ)`` for ``i=1,…,n``.
For more details see section 4.5 of [[HOR19]](@ref).

### Examples

```jldoctest
julia> A = [2..4 -1..1; -1..1 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-1, 1]
 [-1, 1]   [2, 4]

julia> is_strictly_diagonally_dominant(A)
true

julia> A = [2..4 -2..1; -1..2 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-2, 1]
 [-1, 2]   [2, 4]

julia> is_strictly_diagonally_dominant(A)
false
```
"""
function is_strictly_diagonally_dominant(A::AbstractMatrix{T}) where {T<:Interval}
    m, n = size(A)
    m == n || return false

    @inbounds for i=1:m
        sum_mag = sum(Interval(mag(A[i, k])) for k=1:n if k ≠ i)
        mig(A[i, i]) ≤ inf(sum_mag) && return false

    end
      return true
end

"""
    is_Z_matrix(A::AbstractMatrix{T}) where {T<:Interval}

Checks whether the square interval matrix ``A`` is a Z-matrix, that is whether ``Aᵢⱼ≤0``
for all ``i≠j``. For more details see section 4.2 of [[HOR19]](@ref).

### Examples

```jldoctest
julia> A = [2..4 -2.. -1; -2.. -1 2..4]
2×2 Matrix{Interval{Float64}}:
   [2, 4]  [-2, -1]
 [-2, -1]    [2, 4]

julia> is_Z_matrix(A)
true

julia> A = [2..4 -2..1; -1..2 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-2, 1]
 [-1, 2]   [2, 4]

julia> is_Z_matrix(A)
false
```
"""
function is_Z_matrix(A::AbstractMatrix{T}) where {T<:Interval}

    m, n = size(A)
    m == n || return false

    @inbounds for j in 1:n
        for i in 1:m
            if i != j && sup(A[i, j]) > 0
                return false
            end
        end
    end
    return true
end

"""
    is_M_matrix(A::AbstractMatrix{T}) where {T<:Interval}

Checks whether the square interval matrix ``A`` is an M-matrix, that is a Z-matrix with
non-negative inverse. For more details see section 4.2 of [[HOR19]](@ref).

### Examples

```jldoctest
julia> A = [2..2 -1..0; -1..0 2..2]
2×2 Matrix{Interval{Float64}}:
  [2, 2]  [-1, 0]
 [-1, 0]   [2, 2]

julia> is_M_matrix(A)
true

julia> A = [2..4 -2..1; -1..2 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-2, 1]
 [-1, 2]   [2, 4]

julia> is_M_matrix(A)
false
```
"""
function is_M_matrix(A::AbstractMatrix{T}) where {T<:Interval}

    is_Z_matrix(A) || return false

    Ainf = inf.(A)
    e = ones(size(A, 1))
    u = Ainf\e
    all(u .> 0) || return false
    return all(A*u .> 0)
end


text:
################################
####### Define parameters ######
################################
"""
Mutable structure containing the variables required to solve the analytical approach. All the functions are solve using the internal values of the structure. You should declare a mutable structure to the perform the analytical estimations.

# Parameters
 - `gam_neg::Int64`: Selection coefficient for deleterious alleles
 - `gL::Int64`: Selection coefficient for weakly benefical alleles
 - `gH::Int64`: Selection coefficient for strongly benefical alleles
 - `al_low::Float64`: Proportion of α due to weak selection
 - `al_tot::Float64`: α
 - `θ_noncoding::Float64`: Mutation rate defining BGS strength
 - `θ_coding::Float64`: Mutation rate on coding region
 - `al::Float64`: DFE shape parameter
 - `be::Float64`: DFE scale parameter
 - `B::Float64`: BGS strength
 - `B_bins::Array{Float64,1}`: BGS values to simulate
 - `ppos_l::Float64`: Fixation probabily of weakly beneficial alleles
 - `ppos_h::Float64`: Fixation probabily of strongly beneficial alleles
 - `N::Int64`: Population size
 - `n::Int64`: Sample size
 - `flanking_length::Int64`: Flanking region length
 - `rho::Float64`: Recombination rate
 - `TE::Float64`

"""
@with_kw mutable struct parameters

	gam_neg::Int64             = -457
	gL::Int64                  = 10
	gH::Int64                  = 500
	al_low::Float64            = 0.2
	al_tot::Float64            = 0.4
	θ_noncoding::Float64       = 1e-3
	θ_coding::Float64          = 1e-3
	shape::Float64             = 0.184
	scale::Float64             = 0.000402
	B::Float64                 = 0.999
	B_bins::Array{Float64,1}   = push!(collect(0.1:0.025:0.975),0.999)
	ppos_l::Float64            = 0
	ppos_h::Float64            = 0.001
	N::Int64                   = 1000
	n::Int64                   = 500
	flanking_length::Int64      = 2*10^5
	ρ::Float64                 = 0.001
	TE::Float64                = 5.0
	diploid::Bool              = false

	NN::Int64 = 2*N
	nn::Int64 = 2*n

	dac::Array{Int64,1} = [2,4,5,10,20,50,200,500,700]
end

"""
Mutable structure containing the downsampled SFS.

# Returns
 - `bn::Dict`: SparseMatrixCSC containing the binomial convolution

"""
@with_kw mutable struct binomial_dict
	bn::Dict = Dict{Float64,SparseMatrixCSC{Float64,Int64}}()
end

################################
###### Solving parameters ######
################################
"""
	Br(Lmax,theta)

Expected reduction in nucleotide diversity. Explored at [Charlesworth B., 1994](https://doi.org/10.1017/S0016672300032365):

```math
\\frac{\\pi}{\\pi_{0}} = e^{\\frac{4\\muL}{2rL+t}}
```

# Arguments
 - `param::parameters`
 - `theta::Float64`
# Returns
 - `Float64`: expected reduction in diversity given a non-coding length, mutation rate and defined recombination.
"""
function Br(param::parameters,theta::Float64)

	ρ::Float64     = param.ρ
	t::Float64     = -1.0*param.gam_neg/(param.NN+0.0)
	μ::Float64     = theta/(2.0*param.NN)
	r::Float64     = ρ/(2.0*param.NN)

	out::Float64  = ℯ^(-4*μ*param.flanking_length/(2*param.flanking_length*r+t))
	return out
end

# Set mutation rate given the expected reduction in nucleotide diversity (B value ) in a locus.
"""
	set_θ!(param)

Find the optimum mutation given the expected reduction in nucleotide diversity (B value) in a locus.

# Returns
 - `adap.θ_noncoding::Float64`: changes adap.θ_noncoding value.
"""
function set_θ!(param::parameters)

	i(θ,p=param)      = Br(p,θ) - p.B
	theta_f           = find_zero(i,0.0)
	param.θ_noncoding = theta_f
end

function alpha_exp_sim_low(param::parameters,ppos_l::Float64,ppos_h::Float64)
	return fix_pos_sim(param,param.gL,0.5*ppos_l)/(fix_pos_sim(param,param.gL,0.5*ppos_l)+fix_pos_sim(param,param.gH,0.5*ppos_h) + fix_neg_b(param,0.5*ppos_l+0.5*ppos_h))
end

function alpha_exp_sim_tot(param::parameters,ppos_l::Float64,ppos_h::Float64)
	return (fix_pos_sim(param,param.gL,0.5*ppos_l)+fix_pos_sim(param,param.gH,0.5*ppos_h))/(fix_pos_sim(param,param.gL,0.5*ppos_l)+fix_pos_sim(param,param.gH,0.5*ppos_h)+fix_neg_b(param,0.5*ppos_l+0.5*ppos_h))
end

function solv_eqns(param,config)

	ppos_l,ppos_h = config
	return (alpha_exp_sim_tot(param,ppos_l,ppos_h)-param.al_tot,alpha_exp_sim_low(param,ppos_l,ppos_h)-param.al_low)
end

"""
	set_ppos!(param)

Find the probabilty of positive selected alleles given the model. It solves a equation system taking into account fixations probabilities of weakly and strong beneficial alleles.

# Returns
 - `Tuple{Float64,Float64}`: weakly and strong beneficial alleles probabilites.
"""
function set_ppos!(param::parameters)

	function f!(F,x,param=param)
		F[1] = alpha_exp_sim_tot(param,x[1],x[2])-param.al_tot
		F[2] = alpha_exp_sim_low(param,x[1],x[2])-param.al_low
	end

	ppos_l,ppos_h = nlsolve(f!,[0.0; 0.0]).zero

	if ppos_l < 0.0
		 ppos_l = 0.0
	end
	if ppos_h < 0.0
		 ppos_h = 0.0
	end

	param.ppos_l,param.ppos_h = ppos_l, ppos_h
end

"""
	binom_op(param)

Binomial convolution to sample the allele frequencies probabilites depending on background selection values and sample size.

# Arguments
 - `param::parameters`
 - `binom::binomial_dict`

# Returns
 - `Array{Float64,2}`: convoluted SFS for each B value defined in the model (param.B_bins). The estimations are saved at *binom*.
"""
function binom_op!(param::parameters,binom::Dict)

	for b in param.B_bins

		NN2          = convert(Int64,ceil(param.NN*b))
		samples      = collect(1:(param.nn-1))
		p_size        = collect(0:NN2)
		samples_freqs = permutedims(p_size/NN2)
		neutral_sfs   = @. 1/p_size
		replace!(neutral_sfs, Inf => 0.0)

		f(x,y=param.nn) = Binomial(y,x)
		z    = f.(samples_freqs)

		out  = pdf.(z,samples)
		out  = round.(out,digits=10)
		outS = SparseArrays.dropzeros(SparseArrays.sparse(out))
		binom[b] = outS
	end
end

function binom_op!(NN::Int64,nn::Int64,B_bins::Vector{Float64})

	bn = Dict(B_bins[i] => spzeros(nn+1,NN) for i in 1:length(B_bins))

	for b_val in B_bins

		NN2          = convert(Int64,ceil(NN*b_val))
		samples      = collect(1:(nn-1))
		p_size        = collect(0:NN2)
		samples_freqs = permutedims(p_size/NN2)
		neutral_sfs   = @. 1/p_size
		replace!(neutral_sfs, Inf => 0.0)

		f(x,y=nn) = Binomial(y,x)
		z    = f.(samples_freqs)

		out  = pdf.(z,samples)
		out  = round.(out,digits=10)
		out_S = SparseArrays.dropzeros(SparseArrays.sparse(out))

		bn[b_val] = out_S

	end
	return(bn)
end

"""

	Φ(param,gamma,ppos)

Reduction in fixation probabilty due to background selection and linkage. The formulas used have been subjected to several theoretical works ([Charlesworth B., 1994](https://doi.org/10.1017/S0016672300032365), [Hudson et al., 1995](https://www.genetics.org/content/141/4/1605), [Nordborg et al. 1995](https://doi.org/10.1017/S0016672300033619), [Barton NH., 1995](https://www.genetics.org/content/140/2/821)).

The fixation probabilty of selected alleles are reduce by a factor ``\\phi``:

```math
\\phi(t,s) = e^{[\\frac{-2\\mu}{t(1+\\frac{rL}{t}+\\frac{2s}{t})}]}
```

Multiplying across all deleterious linkes sites, we find:

```math
\\Phi = \\prod_{1}^{L} = \\phi(t,s) = e^{\\frac{-2t\\mu(\\psi[1,\\frac{r+2s+L}{r}] - \\psi[1,\\frac{r(L+1)+2s+t}{r}])}{r^2}}
```

```math
\\phi(t,s) = e^{\\frac{-2t\\mu(\\psi[1,\\frac{r+2s+L}{r}] - \\psi[1,\\frac{r(L+1)+2s+t}{r}])}{r^2}}
```

# Arguments
 - `gamma::Int64`: selection coefficient.

# Returns
 - `Float64`: expected rate of positive fixations under background selection.

"""
function Φ(param::parameters,gammaValue::Int64)

	S::Float64 = abs(param.gam_neg/(1.0*param.NN))
	r::Float64 = param.ρ/(2.0*param.NN)
	μ::Float64 = param.θ_noncoding/(2.0*param.NN)
	s::Float64 = gammaValue/(param.NN*1.0)

	Ψ0::Float64 = polygamma(1,(s+S)/r)
	Ψ1::Float64 = polygamma(1,(r+param.flanking_length*r+s+S)/r)
	CC::Float64 = 1.0

	out::Float64 = (ℯ^(-2.0*S*μ*(Ψ0-Ψ1)/(r^2)))
	return out
end

"""
	analytical_alpha(param, binom)

Analytical α(x) estimation. Solve α(x) generally. We used the expected rates of divergence and polymorphism to approach the asympotic value accouting for background selection, weakly and strong positive selection. α(x) can be estimated taking into account the role of positive selected alleles or not. In this way we explore the role of linkage to deleterious alleles in the coding region.

```math
\\mathbb{E}[\\alpha_{x}] =  1 - \\left(\\frac{\\mathbb{E}[D_{s}]}{\\mathbb{E}[D_{N}]}\\frac{\\mathbb{E}[P_{N}]}{\\mathbb{E}[P_{S}]}\\right)
```

# Arguments
 - `param::parameters`
 - `binom::binomial_dict`
# Returns
 - `Array{Float64,1}` α(x).
"""
function analytical_alpha(param::parameters,binom::Dict{Float64, SparseMatrixCSC{Float64, Int64}})

	################################################################
		# Solve the model similarly to original python mktest  #
	################################################################
	B = param.B

	set_θ!(param)
	θ_noncoding = param.θ_noncoding
	# Solve the probabilities of fixations without BGS
	## Set non-bgs
	param.B = 0.999
	## Solve the mutation rate
	set_θ!(param)
	## Solve the probabilities
	set_ppos!(param)
	# Return to the original values
	param.θ_noncoding = θ_noncoding
	param.B = B

	################################################################
	 # Accounting for positive alleles segregating due to linkage #
	################################################################

	# Fixation
	fN     = param.B*fix_neut(param)
	fNeg   = param.B*fix_neg_b(param,0.5*param.ppos_h+0.5*param.ppos_l)
	fPosL  = fix_pos_sim(param,param.gL,0.5*param.ppos_l)
	fPosH  = fix_pos_sim(param,param.gH,0.5*param.ppos_h)

	ds = fN
	dn = fNeg + fPosL + fPosH

	## Polymorphism
	neut = DiscSFSNeutDown(param,binom[param.B])

	selH::Array{Float64,1} = if isinf(exp(param.gH * 2))
		DiscSFSSelPosDownArb(param,param.gH,param.ppos_h,binom[param.B])
	else
		DiscSFSSelPosDown(param,param.gH,param.ppos_h,binom[param.B])
	end

	selL::Array{Float64,1} = DiscSFSSelPosDown(param,param.gL,param.ppos_l,binom[param.B])
	selN::Array{Float64,1} = DiscSFSSelNegDown(param,param.ppos_h+param.ppos_l,binom[param.B])

	splitColumns(matrix::Array{Float64,2}) = (view(matrix, :, i) for i in 1:size(matrix, 2));
	tmp = cumulative_sfs(hcat(neut,selH,selL,selN),false)

	neut, selH, selL, selN = splitColumns(tmp)
	sel = (selH+selL)+selN

	# ps = @. neut / (sel+neut)
	# pn = @. sel / (sel+neut)

	## Outputs
	α = @. 1 - ((ds/dn) * (sel/neut))

	#############################################################
	#		Accounting for neutral and deleterious alleles		#	#############################################################
	## Fixation
	fN_nopos     = fN*(param.θ_coding/2.)*param.TE*param.NN
	fNeg_nopos   = fNeg*(param.θ_coding/2.)*param.TE*param.NN
	fPosL_nopos  = fPosL*(param.θ_coding/2.)*param.TE*param.NN
	fPosH_nopos  = fPosH*(param.θ_coding/2.)*param.TE*param.NN

	ds_nopos = fN_nopos
	dn_nopos = fNeg_nopos + fPosL_nopos + fPosH_nopos

	## Polymorphism
	sel_nopos = selN
	ps_nopos = @. neut / (sel_nopos + neut)
	pn_nopos = @. sel_nopos / (sel_nopos + neut)

	α_nopos = 1 .- (ds_nopos/dn_nopos) .* (sel_nopos./neut)

	##########
	# Output #
	##########
	return (α,α_nopos)
end


text:
mutable struct Ship <: GameObject
    r::Real
    heading::Real
    rotation::Real
    pos::Vector{<:Real}
    vel::Vector{<:Real}
    isboosting::Bool
    @init function Ship()
        this.pos = [@width() / 2, @height() / 2]
        this.r = 20
        this.heading = 0
        this.rotation = 0
        this.vel = [0, 0]
        this.isboosting = false
    end
end

boosting!(this::Ship, val::Bool) = setproperty!(this, :isboosting, val)

function update!(this::Ship)
    this.isboosting && boost!(this)
    this.pos += this.vel
    this.vel *= 0.95
end

function boost!(this::Ship)
    ϕ = this.heading
    force = [cos(ϕ), sin(ϕ)] * 0.3
    this.vel += force
end

function render(this::Ship)
    r = this.r
    x, y = this.pos
    @layer begin
        translate(x, y)
        rotate(this.heading + π/2)

        poly([Point(-r, r), Point(r, r), Point(0, -r)], :path, close=true)
        sethue("#111"); fillpreserve()
        sethue("white"); strokepath()
    end
end

function setrotation!(this::Ship, dθ::Real)
    this.rotation = dθ
end

function turn!(this::Ship)
    this.heading += this.rotation
end

function hits(this::Ship, asteroid::Asteroid)
    d = distance(Point(this.pos...), Point(asteroid.pos...))
    return d < this.r + asteroid.r
end

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
using Revise
using ADCME
using ADCMEKit
using NNFEM
using PyPlot
using ProgressMeter
using AdFem
using SparseArrays
#----------------------------------------------------------------------------
# step 1: create a computational domain
# fixed on the bottom and push from right

m = 40
n = 20
h = 1/n
NT = 50
Δt = 0.01/NT

bdedge = bcedge("right", m, n, h)
bdnode = bcnode("lower", m, n, h)

λ = constant(8.641975308641973e9)
μ = constant(3.7037037037037034e9)
invη = 1/1e5*constant(ones(4*m*n))

M = 2000*compute_fem_mass_matrix1(m, n, h)
Zero = spzeros((m+1)*(n+1), (m+1)*(n+1))
M = SparseTensor([M Zero;Zero M])

## alpha-scheme
β = 1/4; γ = 1/2

# invη is a 4*m*n array
function make_matrix(invη)
  fn_G = invη->begin
    G = tensor([1/Δt+2/3*μ*invη -μ/3*invη 0.0
      -μ/3*invη 1/Δt+2/3*μ*invη 0.0
      0.0 0.0 1/Δt+μ*invη])
    invG = inv(G)
  end
  invG = map(fn_G, invη)
  S = tensor([2μ/Δt+λ/Δt λ/Δt 0.0
      λ/Δt 2μ/Δt+λ/Δt 0.0
      0.0 0.0 μ/Δt])

  H = invG*S
  K = compute_fem_stiffness_matrix(H, m, n, h)
  L = M + β*Δt^2*K
  # L = pl[1]*L
  L, Lbd = fem_impose_Dirichlet_boundary_condition_experimental(L, bdnode, m, n, h)

  return K, L, S, invG
end

function simulate()

  a = TensorArray(NT+1); a = write(a, 1, zeros(2(m+1)*(n+1))|>constant)
  v = TensorArray(NT+1); v = write(v, 1, zeros(2(m+1)*(n+1))|>constant)
  d = TensorArray(NT+1); d = write(d, 1, zeros(2(m+1)*(n+1))|>constant)
  U = TensorArray(NT+1); U = write(U, 1, zeros(2(m+1)*(n+1))|>constant)
  Sigma = TensorArray(NT+1); Sigma = write(Sigma, 1, zeros(4*m*n, 3)|>constant)
  Varepsilon = TensorArray(NT+1); Varepsilon = write(Varepsilon, 1,zeros(4*m*n, 3)|>constant)

  Forces = zeros(NT, 2(m+1)*(n+1))
  for i = 1:NT
    T = eval_f_on_boundary_edge((x,y)->1e6, bdedge, m, n, h)

    T = [-T zeros(size(T)...)]
    rhs = compute_fem_traction_term(T, bdedge, m, n, h)

    Forces[i, :] = rhs
  end
  Forces = constant(Forces)

  K, L, S, invG = make_matrix(invη)

  L = factorize(L)

  function condition(i, tas...)
    i <= NT
  end

  function body(i, tas...)
    a_, v_, d_, U_, Sigma_, Varepsilon_ = tas
    a = read(a_, i)
    v = read(v_, i)
    d = read(d_, i)
    U = read(U_, i)
    Sigma = read(Sigma_, i)
    Varepsilon = read(Varepsilon_, i)

    res = batch_matmul(invG/Δt, Sigma)
    F = compute_strain_energy_term(res, m, n, h) - K * U
    rhs = Forces[i] -  F

    td = d + Δt*v + Δt^2/2*(1-2β)*a
    tv = v + (1-γ)*Δt*a
    rhs = rhs - K*td
    rhs = scatter_update(rhs, constant([bdnode; bdnode.+(m+1)*(n+1)]), constant(zeros(2*length(bdnode))))

    ## alpha-scheme
    a = L\rhs # bottleneck
    d = td + β*Δt^2*a
    v = tv + γ*Δt*a
    U_new = d

    Varepsilon_new = eval_strain_on_gauss_pts(U_new, m, n, h)

    res2 = batch_matmul(invG * S, Varepsilon_new-Varepsilon)
    Sigma_new = res +  res2

    i+1, write(a_, i+1, a), write(v_, i+1, v), write(d_, i+1, d), write(U_, i+1, U_new),
          write(Sigma_, i+1, Sigma_new), write(Varepsilon_, i+1, Varepsilon_new)
  end

  i = constant(1, dtype=Int32)
  _, _, _, _, u, sigma, varepsilon = while_loop(condition, body,
                    [i, a, v, d, U, Sigma, Varepsilon])

  U = stack(u)
  Sigma = stack(sigma)
  Varepsilon = stack(varepsilon)

  U = set_shape(U, (NT+1, size(U,2)))
  Sigma = set_shape(Sigma, (NT+1, 4*m*n, 3))

  return U, Sigma
end

U, sigma = simulate()
sess = Session(); init(sess)
U, sigma = run(sess, [U, sigma])

##

elements = []
# prop = Dict("name"=> "PlaneStrain", "rho"=> 2000, "E"=> 1e10, "nu"=> 0.35)
prop = Dict("name"=> "ViscoelasticityMaxwell", "rho"=> 2000, "E"=> 1e10, "nu"=> 0.35, "eta"=>1e7)

coords = zeros((m+1)*(n+1), 2)
for j = 1:n
    for i = 1:m
        idx = (m+1)*(j-1)+i
        elnodes = [idx; idx+1; idx+1+m+1; idx+m+1]
        ngp = 2
        nodes = [
        (i-1)*h (j-1)*h
        i*h (j-1)*h
        i*h j*h
        (i-1)*h j*h
        ]
        coords[elnodes, :] = nodes
        push!(elements, SmallStrainContinuum(nodes, elnodes, prop, ngp))
    end
end

# fixed on the bottom, push on the right
EBC = zeros(Int64, (m+1)*(n+1), 2)
FBC = zeros(Int64, (m+1)*(n+1), 2)
g = zeros((m+1)*(n+1), 2)
f = zeros((m+1)*(n+1), 2)

for i = 1:m+1
    EBC[i+n*(m+1),:] .= -1
end

Edge_Traction_Data = Array{Int64}[]
for i = 1:n
  elem = elements[(i-1)*m + m]
  for k = 1:4
    if elem.coords[k,1]>2-0.001 && elem.coords[k+1>4 ? 1 : k+1,1]>2-0.001
      push!(Edge_Traction_Data, [(i-1)*m + m, k, 1])
    end
  end
end
Edge_Traction_Data = hcat(Edge_Traction_Data...)'|>Array

ndims = 2
domain = Domain(coords, elements, ndims, EBC, g, FBC, f, Edge_Traction_Data)
######

domain.history["state"] = []
domain.history["stress"] = []
for i = 1:NT+1
  push!(domain.history["state"], U[i,:])
  push!(domain.history["stress"], sigma[i,:,:])
end
# visualize_von_mises_stress_on_scoped_body(d_, domain, scale_factor=10.0)
p = visualize_total_deformation_on_scoped_body(U, domain, scale_factor=10.0)
saveanim(p, "adfem.gif")
close("all")

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
"""
    BackendServicesConfig{PK}(; kwargs...)

## Required Keyword Arguments:
- `iss::String`
- `private_key::PK`
- `sub::String`
- `token_endpoint::String`

## Optional Keyword Arguments:
- `scope::Union{String, Nothing}`. Default value: `nothing`.
"""
Base.@kwdef struct BackendServicesConfig{PK <: JSONWebTokens.Encoding}
    iss::String
    private_key::PK
    scope::Union{String, Nothing} = nothing
    sub::String
    token_endpoint::String
end

Base.@kwdef struct BackendServicesResult
    access_token::String
    access_token_is_jwt::Bool = false
    access_token_jwt_decoded::Dict{String, Any} = Dict{String, Any}()
    access_token_response::Dict{Symbol, Any}
end

"""
    HealthBase.get_fhir_access_token(result::BackendServicesResult)
"""
function HealthBase.get_fhir_access_token(result::BackendServicesResult)
    return result.access_token
end


text:
@deprecate testany(index::T, mask::T) where {T <: Integer} anyone(index, mask)
@deprecate testall(index::T, mask::T) where {T <: Integer} allone(index, mask)
@deprecate testval(index::T, mask::T, onemask::T) where {T <: Integer} ismatch(index, mask, onemask)
@deprecate takebit(index::T, ibit::Int) where T <: Integer  readbit(index, ibit)
@deprecate takebit(index::T, bits::Int...) where T <: Integer  readbit(index, bits...)
@deprecate breflect(nbits::Int, b::Integer)  breflect(b; nbits=nbits)
@deprecate breflect(nbits::Int, b::T, masks::Vector{T}) where T breflect(b, masks; nbits=nbits)
@deprecate bit(b; len) BitStr{len}(b)

@deprecate to_location(x::Integer) x+1
@deprecate to_location(x::BitStr) buffer(x)+1


text:
@testset "utfnorm" begin
    using DelimitedFiles

    hexcode2unicode(x) =  transcode(String, map(x->parse(UInt32, x;base=16), split(x)))

    nr = BytePairEncoding.UtfNormalizer(:NFKC)
    xnr = BytePairEncoding.UnNormalizer()

    tbl = readdlm(joinpath(dirname(@__FILE__), "data/nfkc.tsv"), '\t', String; comments=true)
    for i ∈ 1:size(tbl)[1]
        lhs = hexcode2unicode(tbl[i, 1])
        rhs = hexcode2unicode(tbl[i, 2])
        @test normalize(nr, lhs) == rhs
        @test normalize(xnr, lhs) == lhs
        @test normalize(xnr, rhs) == rhs
    end
end