text:
"""
Parallel Factor Analysis (PARAFAC) decomposition
"""
struct PARAFAC2{T<:Number, N} <: TensorDecomposition{T, N}
    factors::NTuple{N, Matrix{T}}
    B::Matrix{T}
    A::Matrix{T}
    props::Dict{Symbol, Any}

    function PARAFAC2(
                   X::Vector{Matrix{T}},
                   F::Matrix{T},
                   B::Matrix{T},
                   A::Matrix{T}) where T
        N = length(X)
        @assert size(B, 2) == N
        factors = ntuple(N) do i
            u = svd((F .* view(B, :, i)) * (X[i] * A)')
            return u.V * (u.U'F)
        end
        return new{T, N}(factors, B, A, Dict{Symbol,Any}())
    end

end

"""
PARAFAC2 model
"""
function parafac2(X::AbstractVector{<:StridedMatrix{T}},
                  r::Integer;
                  tol::Float64=1e-5,
                  maxiter::Integer=100,
                  compute_error::Bool=false,
                  verbose::Bool=true) where T
    m = length(X)
    n = size(X[1], 2)
    all(Xi -> size(Xi, 2) == n, X) ||
        throw(DimensionMismatch("All input matrices should have the same number of rows."))

    F = Matrix{T}(I, r, r)
    B = fill(T(1.0), m, r)
    A = eigen(Symmetric(sum(Xi -> Xi'Xi, X)), (n-r+1):n).vectors
    G1 = Matrix{T}(I, r, r); G2 = copy(G1); G3 = fill(T(m), r, r)
    H = Matrix{T}[(size(X[i], 1) > size(X[i], 2) ? qr(X[i]).R : X[i]) for i in 1:m]
    P = [Matrix{T}(undef, r, size(H[i], 1)) for i in 1:m]
    PH = Array{T, 3}(undef, r, n, m)

    niters = 0
    converged = false
    resid = resid0 = norm(X)
    pb = Progress(maxiter, "PARAFAC2 iterations ")
    while !converged && niters < maxiter
        update!(pb, niters)
        for i in eachindex(H)
            u = svd((F .* view(B, i, :)) * (H[i] * A)')
            mul!(P[i], u.U, u.Vt)
            mul!(view(PH, :, :, i), P[i], H[i])
        end
        mul!(F, _row_unfold(PH, 1), khatrirao(B, A) / (G3 .* G2))
        #F = _row_unfold(PH, 1) * khatrirao(B, A) / (G3 .* G2)
        mul!(G1, F', F)
        mul!(A, _row_unfold(PH, 2), khatrirao(B, F) / (G3 .* G1))
        #A = _row_unfold(PH, 2) * khatrirao(B, F) / (G3 .* G1)
        mul!(G2, A', A)
        mul!(B, _row_unfold(PH, 3), khatrirao(A, F) / (G2 .* G1))
        #B = _row_unfold(PH, 3) * khatrirao(A, F) / (G2 .* G1)
        mul!(G3, B', B)

        resid_old = resid
        resid = sum(k -> sum(abs2, H[k] .- P[k]' * ((F .* view(B, k, :)) * A')), 1:m)
        converged = abs(resid - resid_old) < tol * resid_old

        niters += 1
    end
    finish!(pb)
    verbose && _iter_status(converged, niters, maxiter)
    res = PARAFAC2(X, F, Matrix(B'), A)
    compute_error && _set_rel_residue(res, sqrt(resid) / resid0)
    return res
end


text:
"""
Point symmetry database module.
Provides matrix representations (in units of lattice vectors) and other information about the point symmetries
(Hermann-Mauguin notation, names of the elements, etc.), together with character tables and irreducible representations.
"""
module PointSymmetryDatabase

import Serialization
import Pkg

using YAML

import ..LatticeTools: PointSymmetry, read_point_symmetry, simplify_name

NUM_POINT_SYMMETRIES_2D = 11
POINT_SYMMETRY_DATABASE_2D = Vector{PointSymmetry}(undef, NUM_POINT_SYMMETRIES_2D)
POINT_SYMMETRY_LOOKUP_2D = Dict{Vector{String}, Int}()

NUM_POINT_SYMMETRIES_3D = 32
POINT_SYMMETRY_DATABASE_3D = Vector{PointSymmetry}(undef, NUM_POINT_SYMMETRIES_3D)
POINT_SYMMETRY_LOOKUP_3D = Dict{Vector{String}, Int}()

# COV_EXCL_START
function load_yaml()
    for group_index in 1:NUM_POINT_SYMMETRIES_2D
        psym = load_group_2d(group_index)
        global POINT_SYMMETRY_DATABASE_2D[group_index] = psym
        global POINT_SYMMETRY_LOOKUP_2D[sort(simplify_name.(psym.element_names))] = group_index
    end
    for group_index in 1:NUM_POINT_SYMMETRIES_3D
        psym = load_group_3d(group_index)
        global POINT_SYMMETRY_DATABASE_3D[group_index] = psym
        global POINT_SYMMETRY_LOOKUP_3D[sort(simplify_name.(psym.element_names))] = group_index
    end
end

function load()
    mutable_artifacts_toml = joinpath(@__DIR__, "..", "..", "MutableArtifacts.toml")
    try
        cache_hash = Pkg.Artifacts.artifact_hash("PointSymmetryDatabase", mutable_artifacts_toml)
        if !isnothing(cache_hash) && Pkg.Artifacts.artifact_exists(cache_hash)
            cache_filepath = joinpath(Pkg.Artifacts.artifact_path(cache_hash), "PointSymmetryDatabase.cache")
            cache = Serialization.deserialize(cache_filepath)
            global POINT_SYMMETRY_DATABASE_2D = cache.database2d
            global POINT_SYMMETRY_DATABASE_3D = cache.database3d
            global POINT_SYMMETRY_LOOKUP_2D = cache.lookup2d
            global POINT_SYMMETRY_LOOKUP_3D = cache.lookup3d
            return
        end
    catch e
        @warn "Failed to load PointSymmetryDatabase."
        @warn "$e"
        @warn "Loading YAML"
    end

    load_yaml()
    cache_hash = Pkg.Artifacts.create_artifact() do cache_directory
        global POINT_SYMMETRY_DATABASE_2D
        global POINT_SYMMETRY_DATABASE_3D
        global POINT_SYMMETRY_LOOKUP_2D
        global POINT_SYMMETRY_LOOKUP_3D
        cache_filepath = joinpath(cache_directory, "PointSymmetryDatabase.cache")
        Serialization.serialize(cache_filepath,
                                (database2d=POINT_SYMMETRY_DATABASE_2D,
                                    lookup2d=POINT_SYMMETRY_LOOKUP_2D,
                                    database3d=POINT_SYMMETRY_DATABASE_3D,
                                    lookup3d=POINT_SYMMETRY_LOOKUP_3D))
    end
    Pkg.Artifacts.bind_artifact!(mutable_artifacts_toml, "PointSymmetryDatabase", cache_hash; force=true)
end

function __init__()
    load()
end

function load_group_2d(group_index::Integer)
    data_directory = abspath(joinpath(@__DIR__, "..", "..", "data", "PointGroup2D"))
    file_path = joinpath(data_directory, "PointGroup2D-$group_index.yaml")
    data_yaml = YAML.load_file(file_path)
    point_symmetry = read_point_symmetry(data_yaml)
    return point_symmetry
end

function load_group_3d(group_index::Integer)
    data_directory = abspath(joinpath(@__DIR__, "..", "..", "data", "PointGroup3D"))
    file_path = joinpath(data_directory, "PointGroup3D-$group_index.yaml")
    data_yaml = YAML.load_file(file_path)
    point_symmetry = read_point_symmetry(data_yaml)
    return point_symmetry
end
# COV_EXCL_STOP

# 3D by default
get(args...) = get3d(args...)
find(args...) = find3d(args...)
findindex(args...) = findindex3d(args...)

"""
    get2d(group_index::Integer)

Get the `group_index`th two-dimensional point group symmetry.
Ref: [https://www.cryst.ehu.es]
"""
function get2d(group_index::Integer)
    if !(1 <= group_index <= NUM_POINT_SYMMETRIES_2D)
        throw(ArgumentError("Point group 2D #$group_index not found"))
    end
    return POINT_SYMMETRY_DATABASE_2D[group_index]
end

"""
    get3d(group_index::Integer)

Get the `group_index`th three-dimensional point group symmetry.
Ref: [https://www.cryst.ehu.es]
"""
function get3d(group_index::Integer)
    if !(1 <= group_index <= NUM_POINT_SYMMETRIES_3D)
        throw(ArgumentError("Point group 3D #$group_index not found"))
    end
    return POINT_SYMMETRY_DATABASE_3D[group_index]
end

"""
    find2d(group_name::AbstractString)

Find a two-dimensional point group symmetry by Hermann-Mauguin name.
Ref: https://www.cryst.ehu.es/
"""
function find2d(group_name::AbstractString)
    for i in 1:NUM_POINT_SYMMETRIES_2D
        psym = get2d(i)
        if group_name == psym.hermann_mauguin
            return psym
        end
    end
    return nothing
end

"""
    find3d(group_name::AbstractString)

Find a three-dimensional point group symmetry by Hermann-Mauguin name.
Ref: https://www.cryst.ehu.es/
"""
function find3d(group_name::AbstractString)
    for i in 1:NUM_POINT_SYMMETRIES_3D
        psym = get3d(i)
        if group_name == psym.hermann_mauguin
            return psym
        end
    end
    return nothing
end

"""
    find2d(element_names::AbstractVector{<:AbstractString})

Find a two-dimensional point group symmetry by element names.
The names will be "simplified", without any subscripts.
"""
function findindex2d(element_names::AbstractVector{<:AbstractString})
    simple_names = sort(simplify_name.(element_names))
    return POINT_SYMMETRY_LOOKUP_2D[simple_names]
end

function find2d(element_names::AbstractVector{<:AbstractString})
    return get2d(findindex2d(element_names))
end

"""
    find3d(element_names::AbstractVector{<:AbstractString})

Find a three-dimensional point group symmetry by element names.
The names will be "simplified", without any subscripts.

```julia
PointSymmetry
````
"""
function findindex3d(element_names::AbstractVector{<:AbstractString})
    simple_names = sort(simplify_name.(element_names))
    return POINT_SYMMETRY_LOOKUP_3D[simple_names]
end

function find3d(element_names::AbstractVector{<:AbstractString})
    return get3d(findindex3d(element_names))
end

end # module PointSymmetryDatabase

ext: jl


text:
#!/usr/bin/env julia

using BenchmarkTools

function caret(x, y)
  z = x ^ y

  return z
end

o = @benchmark caret(10, -3.)
println(minimum(o).time, " nanoseconds")

max_stars_repo_name: scienceopen/numba-examples

max_issues_repo_name: scienceopen/numba-examples

max_forks_repo_name: scienceopen/numba-examples


text:
module TimeManager

export currentTime
export programTime
export programTimeStr
export programStartTime
export OnTime

now() = Dates.time()
programStartTime = now()

currentTime(startTime::Real) = (now() - startTime)
programTime() = currentTime(programStartTime)
programTimeStr() = @sprintf("%.3f", programTime())

function OnTime(milisec::Float64, prevTime::Base.RefValue{Float64})
	time=now()
	r=(time - Base.getindex(prevTime)) >= milisec
	if r Base.setindex!(prevTime, time) end
	r
end

#Dates.millisecond(Dates.unix2datetime(Dates.time()))

end

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
push!(LOAD_PATH, string(pwd(),"/src"))
using ModifiedStokesSolver
using Base.Test

@testset "FMM" begin
    ## define source and target points
    ns = 100
    nt = 102

    srand(1)
    src = 1 - 2*rand(2,ns)
    fvec = 1 - 2*rand(2,ns)
    nvec = 1 - 2*rand(2,ns)
    targ = 1 - 2*rand(2,nt)

    # Compute using the kernels from ModifiedStokesSolver
    function mss_stokeslet_direct(src, targ, fvec, alpha)
        ns = size(src, 2)
        nt = size(targ, 2)
        u = zeros(2, nt)
        for i=1:ns
            for j=1:nt
                f = fvec[:, i]
                r = src[:,i]-targ[:,j]
                u[:,j] += stokeslet(r, f, alpha)
            end
        end
        return u
    end
    function mss_stresslet_direct(src, targ, fvec, nvec, alpha)
        ns = size(src, 2)
        nt = size(targ, 2)
        u = zeros(2, nt)
        for i=1:ns
            for j=1:nt
                r = src[:,i]-targ[:,j]
                u[:,j] += stresslet(r, fvec[:,i], nvec[:,i], alpha)
            end
        end
        return u
    end
    function mss_stresslet_direct_self(src, fvec, nvec, alpha)
        ns = size(src, 2)
        u = zeros(2, ns)
        for i=1:ns
            for j=1:ns
                i==j && continue
                r = src[:,i]-src[:,j]
                u[:,j] += stresslet(r, fvec[:,i], nvec[:,i], alpha)
            end
        end
        return u
    end

    ## Test with large alpha
    alpha = 20.0

    ## Stokeslet
    u1 = mss_stokeslet_direct(src, targ, fvec, alpha)
    u2 = fmm_stokeslet_direct(src, targ, fvec, alpha)
    u3 = fmm_stokeslet_targ(src, targ, fvec, alpha)
    println(" * Stokeslet, alpha=", alpha)
    Gerrdirect = norm(u1[:]-u2[:]) / norm(u1[:], Inf)
    Gerrfmm = norm(u1[:]-u3[:]) / norm(u1[:], Inf)
    @test Gerrdirect < 5e-14
    @test Gerrfmm < 5e-14

    ## Stresslet
    ## Test also with small alpha
    ## (Stokeslet has no special eval for small args)
    for alpha=[alpha, 1e-15]
        println(" * Stresslet, alpha=", alpha)
        # Stresslet targ
        d1 = mss_stresslet_direct(src, targ, fvec, nvec, alpha)
        d2 = fmm_stresslet_direct(src, targ, fvec, nvec, alpha)
        d3 = fmm_stresslet_targ(src, targ, fvec, nvec, alpha)
        # Stresslet self
        ds1 = mss_stresslet_direct_self(src, fvec, nvec, alpha)
        ds2 = fmm_stresslet_direct(src, src, fvec, nvec, alpha; self=true)
        ds3 = fmm_stresslet_self(src, fvec, nvec, alpha)
        # Test
        Terrdirect = norm(d1[:]-d2[:]) / norm(d1[:], Inf)
        Terrdirect_self = norm(ds1[:]-ds2[:]) / norm(d1[:], Inf)
        Terrfmm = norm(d1[:]-d3[:]) / norm(d1[:], Inf)
        Terrfmm_self = norm(ds1[:]-ds3[:]) / norm(ds1[:], Inf)
        @test Terrdirect < 5e-14
        @test Terrdirect_self < 5e-14
        @test Terrfmm < 5e-14
        @test Terrfmm_self < 5e-14
    end
end


text:
struct Tuners{F <: AbstractFloat}
  nu::Vector{MCMC.TuningParam{F}}
end

function Tuners(K, value=1.0)
  nu_tuner = [MCMC.TuningParam(value) for _ in 1:K]
  return Tuners(nu_tuner)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
"""
    from_json(::Type{T},s::AbstractString) where T
"""
function from_json(::Type{T},s::AbstractString) where T
  json_dict = JSON.parse(s)
  dict = _decode_json_dict(json_dict)
  from_dict(T,dict)
end

"""
    from_json_file(::Type{T},s::AbstractString) where T
"""
function from_json_file(::Type{T},s::AbstractString) where T
  json_dict = JSON.parsefile(s)
  dict = _decode_json_dict(json_dict)
  from_dict(T,dict)
end

"""
    to_json(object)
"""
function to_json(object)
  dict = to_dict(object)
  JSON.json(dict)
end

"""
    to_json_file(object,filename)
"""
function to_json_file(object,filename)
  dict = to_dict(object)
  open(filename,"w") do f
    JSON.print(f,dict)
  end
end

JSON.lower(object::GridapType) = to_dict(object)

function _decode_json_dict(json_dict::Dict{String,Any})
  dict = Dict{Symbol,Any}()
  for k in keys(json_dict)
    v = json_dict[k]
    if isa(v,Dict{String,Any})
      w = _decode_json_dict(v)
    elseif isa(v,AbstractVector)
      _w = []
      for vi in v
        if isa(vi,Dict{String,Any})
          wi = _decode_json_dict(vi)
        else
          wi = vi
        end
        push!(_w,wi)
      end
      w = collect(_w)
    else
      w = v
    end
    dict[Symbol(k)] = w
  end
  dict
end


text:
const Read = Dict("1"=>1,"2"=>2,"3"=>3,"4"=>4,"5"=>5,"6"=>6,"7"=>7,"8"=>8,"9"=>9,"a"=>10,
            "b"=>11,"c"=>12,"d"=>13,"e"=>14,"f"=>15,"g"=>16,"h"=>17,"i"=>18,"j"=>19,"k"=>20,
            "l"=>21,"m"=>22,"n"=>23,"o"=>24,"p"=>25,"q"=>26,"r"=>27,"s"=>28,"t"=>29,"u"=>30)
#="""
The structure `graph_coloring` define the colorations of a graph.
# Arguments
- `graph::SimpleGraph`: the graph.
- `num_cols::UInt8`: the number of colors.
- `current_color::Array{UInt8,1}`: the array of colors.
"""=#
mutable struct graph_coloring
    graph::SimpleGraph
    num_cols::UInt8
    current_color::Array{UInt8,1}
end

function graph_coloring(G::SimpleGraph, num_cols::UInt8)::graph_coloring
    return graph_coloring(G, num_cols, UInt8[])
end

mutable struct graph_coloring_from_file
    file_name::String
    #color_file::IOStream
    color_file_gz::CodecZlib.TranscodingStreams.TranscodingStream{GzipDecompressor, IOStream}
    current_color::Array{UInt8,1}
    current_aut::Int64
end

# function graph_coloring_from_file( file_name::String )
#     color_file = open( file_name )
#     return graph_coloring_from_file( file_name, color_file, [], 1 )
# end

function graph_coloring_from_file(file_name::String)::graph_coloring_from_file

    color_file_gz = GzipDecompressorStream(open(file_name))

    return graph_coloring_from_file(file_name, color_file_gz, UInt8[], 1)
end

function Base.iterate( GC::graph_coloring_from_file, ::Int64 = 0)::Union{Nothing, Tuple{Vector{UInt8}, Int64}}

    st = readline(GC.color_file_gz)
    if st == "STOP"
        close(GC.color_file_gz)
        return nothing
    end

    st = split(st, ",")
    #GC.current_color = [ parse( UInt8, Read[x] ) for x in split( st[1], "" )]
    GC.current_color = [UInt8(Read[x]) for x in split(st[1], "")]
    GC.current_aut = parse(Int64, st[2])

    return GC.current_color, 0
end

function exists_file_with_colorings(pruf_str::String, n::Int64)::Tuple{Bool,Union{String,Nothing}}

    dir = dirname(current_dir)*"/Data/Dimension$n/" #path of the folder containing the colorations

    if isdir(dir) && pruf_str*"0.gz" in readdir(dir)
        return true, dir*pruf_str*"0.gz"
    else
        return false, nothing
    end
end

function is_coloring(G::SimpleGraph, cols::Array{UInt8,1})::Bool

    num_v = nv( G )
    for i in 1:num_v
        for j in (i+1):num_v
            if cols[i] == cols[j] && has_edge(G, i, j)
                return false
            end
        end
    end

    return true
end

function smallest_coloring(G::SimpleGraph, num_cols::UInt8)::Array{UInt8,1}

    cols = ones(UInt8,nv(G))
    while true
        if is_coloring(G, cols)
            return cols
        end

        next_tuple!(cols, num_cols)
    end
end

function Base.iterate(GC::graph_coloring, ::Int64 = 0)::Union{Nothing, Tuple{Vector{UInt8}, Int64}}

    if isempty(GC.current_color)
        GC.current_color = smallest_coloring(GC.graph, GC.num_cols)
        return GC.current_color, 0
    end

    while true
        v = next_tuple!(GC.current_color, GC.num_cols)
        if v == false
            return nothing
        end

        if is_coloring(GC.graph, GC.current_color)
            return GC.current_color, 0
        end

    end
end

function next_tuple!(tuple::Array{UInt8,1}, max_entry::UInt8)::Union{Nothing, Bool}

    k = length(tuple)
    for k in k:-1:0

        if k == 0
            return false
        elseif tuple[k] < max_entry
            tuple[k] += 1
            break
        else
            tuple[k] = 1
        end
    end

    return nothing
end

#=
"""
Compute the graph with a specific Prufer sequence.
"""=#
function PruferToGraph(prufer::Vector{UInt8})::SimpleGraph

    num_v = length(prufer)+2 #number of vertices
    g = Graph(num_v) #initialize a graph with specific number of vertices

    degree::Vector{Int64} = zeros(Int64,num_v)

    for i in 1:num_v-2
        degree[prufer[i]] += 1
    end

    for i in 1:(num_v-2)
        for j in 1:num_v
            if degree[j] == 0 #If j is not present in prufer set
                add_edge!(g, j, prufer[i]) #add edge
                degree[j] = -1 #Remove from Prufer
                degree[prufer[i]] -= 1 #low the degree of the vertex i
                break
            end
        end
    end

    # For the last element, we look at those vertices still with degree
    i = indexin(0, degree)[1]
    j = num_v + 1 - indexin(0, reverse(degree))[1]
    add_edge!(g, i, j)

    return g
end
#=
"""
Extract a Prufer sequence and a number from a string. The number is the number of automorphisms of the graph with that Prufer sequence.
"""=#
function get_graph(str::String)::Tuple{SimpleGraph{Int64}, Int64}

    s = split(str, ',')
    g = PruferToGraph(UInt8[Read[string(s[1][i])] for i in 1:length(s[1])]) #get the graph
    a = parse(Int64,s[2])

    return (g, a)
end
#=
"""
Return all the arrays of length `l` such that the sum of all elements of the array is `d`.
"""=#
function get_weights(l::Int64, d::Int64)::Vector{Vector{Int64}}

    return vcat(unique.(permutations.(partitions(d,l)))...)
end


text:
using DelimitedFiles
using Test

#=
Refactor after seeing Tony Hallam's solution and wanting to tidy up with
functions and tests. Might try using doctests tomorrow.

I'm sure this can be solved in a really elegant recursive way but I need
to think about it some more. (Recursion is a state of mind, it seems.)
=#

"""
    part_1(data, target)

Search for the number making up to 2020.
"""
function part_1(data::Array, target::Integer=2020)::Integer
    for d in data
        if (target - d) in data
            return convert(Integer, d * (target - d))
        end
    end
end

test = readdlm("./data/day01_test.txt")
@test part_1(test) == 514579

data = readdlm("./data/day01.txt")
println("Part 1: $(part_1(data))")

"""
    part_2(data, target)

I'm sure this could be recursive but I can't quite get it...
"""
function part_2(data::Array, target::Integer=2020)::Integer
    for d in data
        for e in [i for i in data if ((d + i) < target)]
            if (target - d - e) in data
                return convert(Integer, d * e * (target - d - e))
            end
        end
    end
end

@test part_2(test) == 241861950

println("Part 2: $(part_2(data))")


text:
# Copyright 2019 Albin Severinson
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module FountainCodes

using Random, StatsBase, Statistics, Distributions
using LinearAlgebra, SparseArrays
using Primes, DataStructures

export CoefficientType, Binary, NonBinary, Code

# type system
abstract type AbstractErasureCode end

# For testing
export TestConstraints, subtract!, rand_nonzero
struct TestConstraints{Tc,Tv} <: AbstractVector{Tv}
    A::Matrix{Tc}
    Vs::Vector{Tv}
    function TestConstraints(A, Vs)
        _, n = size(A)
        length(Vs) == n || throw(DimensionMismatch("Vs has dimension $length(Vs)), but A has dimensions $(size(A))"))
        new{eltype(A),eltype(Vs)}(Matrix(A), deepcopy(Vs))
    end
end

Base.length(tc::TestConstraints) = length(tc.Vs)
Base.size(tc::TestConstraints) = (length(tc),)
Base.getindex(tc::TestConstraints, args...) = getindex(tc.Vs, args...)

function subtract!(tc::TestConstraints; coef, rpi_src, rpi_dst)
    tc.A[:, rpi_dst] .-= coef .* view(tc.A, :, rpi_src)
    tc.Vs[rpi_dst] -= coef * tc.Vs[rpi_src]
    return
end

function rand_nonzero(rng::AbstractRNG, T::Type)
    rv = rand(rng, T)
    while iszero(rv)
        rv = rand(rng, T)
    end
    rv
end

rand_nonzero(rng::AbstractRNG, T::Type{<:AbstractFloat}) = randn(rng, T)

function rand_nonzero(rng::AbstractRNG, T::Type, dims...)
    rv = zeros(T, dims...)
    for i in 1:length(rv)
        rv[i] = rand_nonzero(rng, T)
    end
    rv
end

export dimension, constraint_matrix, generator_matrix

include("GF256.jl")
include("QMatrix.jl") # currently unused
include("Decode.jl")
include("Gray.jl")
include("R10Tables.jl")
include("R10.jl")
include("LT.jl")
include("RQTables.jl")
include("RQ.jl")

end # module

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
#= Load these if running these tests in terminal
import Pkg; Pkg.activate("../")  # CytofDensityEstimation
include("runtests.jl")
=#
# ApproximateTwoSampleKSTest(a, b)

@testset "Test update!" begin
  K = 3
  # use_big_data = true
  use_big_data = false

  if use_big_data
    data = CDE.Model.Data([fill(-Inf, 30000); randn(70000)],
                          [fill(-Inf, 20000); randn(60000)])
  else
    data = CDE.Model.Data([-Inf, 1.0, 2.0], [-1.0, -Inf, 3.0, 4.0])
  end

  prior = CDE.Model.Prior(K)
  state = CDE.Model.State(data, prior)
  tuners = CDE.Model.Tuners(K)

  @time for i in ProgressBar(1:100)
    CDE.Model.update_state!(state, data, prior, tuners)
  end

  # Test flags
  flags = [CDE.UpdateBetaWithSkewT(), CDE.UpdateLambdaWithSkewT()]
  mu_orig = deepcopy(state.mu)
  nu_orig = deepcopy(state.nu)
  @time for i in ProgressBar(1:100)
    CDE.Model.update_state!(state, data, prior, tuners, fix=[:mu], flags=flags)
  end
  @test mu_orig == state.mu
  @test nu_orig != state.nu
end


text:
"""
    AnisotropicDiffusivity{NX, NY, NZ, KX, KY, KZ}

Parameters for anisotropic diffusivity models.
"""
struct AnisotropicDiffusivity{NX, NY, NZ, KX, KY, KZ} <: AbstractTensorDiffusivity
    νx :: NX
    νy :: NY
    νz :: NZ
    κx :: KX
    κy :: KY
    κz :: KZ
end

"""
    AnisotropicDiffusivity(; νx=ν₀, νy=ν₀, νz=ν₀, κx=κ₀, κy=κ₀, κz=κ₀,
                             νh=nothing, κh=nothing)

Returns parameters for a closure with a diagonal diffusivity tensor with heterogeneous
'anisotropic' components labeled by `x`, `y`, `z`.
Each component may be a number or function.
The tracer diffusivities `κx`, `κy`, and `κz` may be `NamedTuple`s with fields corresponding
to each tracer, or a single number or function to be a applied to all tracers.

If `νh` or `κh` are provided, then `νx = νy = νh`, and `κx = κy = κh`, respectively.

By default, a viscosity of `ν₀ = 1.05×10⁻⁶` m² s⁻¹ is used for all viscosity components
and a diffusivity of `κ₀ = 1.46×10⁻⁷` m² s⁻¹ is used for all diffusivity components for every tracer.
These values are the approximate viscosity and thermal diffusivity for seawater at 20°C
and 35 psu, according to Sharqawy et al., "Thermophysical properties of seawater: A review
of existing correlations and data" (2010).
"""
function AnisotropicDiffusivity(FT=Float64; νx=ν₀, νy=ν₀, νz=ν₀, κx=κ₀, κy=κ₀, κz=κ₀, νh=nothing, κh=nothing)
    if !isnothing(νh)
        νx = νh
        νy = νh
    end

    if !isnothing(κh)
        κx = κh
        κy = κh
    end

    if all(isa.((νx, νy, νz, κx, κy, κz), Number))
        κx = convert_diffusivity(FT, κx)
        κy = convert_diffusivity(FT, κy)
        κz = convert_diffusivity(FT, κz)
        return AnisotropicDiffusivity(FT(νx), FT(νy), FT(νz), κx, κy, κz)
    else
        return AnisotropicDiffusivity(νx, νy, νz, κx, κy, κz)
    end
end

function with_tracers(tracers, closure::AnisotropicDiffusivity)
    κx = tracer_diffusivities(tracers, closure.κx)
    κy = tracer_diffusivities(tracers, closure.κy)
    κz = tracer_diffusivities(tracers, closure.κz)
    return AnisotropicDiffusivity(closure.νx, closure.νy, closure.νz, κx, κy, κz)
end

calculate_diffusivities!(K, arch, grid, closure::AnisotropicDiffusivity, args...) = nothing

@inline ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, clock, closure::AnisotropicDiffusivity, U, args...) =
    ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid, clock, closure.νx, closure.νy, closure.νz, U.u)

@inline ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, clock, closure::AnisotropicDiffusivity, U, args...) =
    ∂ⱼνᵢⱼ∂ᵢv(i, j, k, grid, clock, closure.νx, closure.νy, closure.νz, U.v)

@inline ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, clock, closure::AnisotropicDiffusivity, U, args...) =
    ∂ⱼνᵢⱼ∂ᵢw(i, j, k, grid, clock, closure.νx, closure.νy, closure.νz, U.w)

@inline function ∇_κ_∇c(i, j, k, grid, clock, closure::AnisotropicDiffusivity,
                        c, ::Val{tracer_index}, args...) where tracer_index

    @inbounds κx = closure.κx[tracer_index]
    @inbounds κy = closure.κy[tracer_index]
    @inbounds κz = closure.κz[tracer_index]

    return ∂ⱼκᵢⱼ∂ᵢc(i, j, k, grid, clock, κx, κy, κz, c)
end

Base.show(io::IO, closure::AnisotropicDiffusivity) =
    print(io, "AnisotropicDiffusivity: " *
              "(νx=$(closure.νx), νy=$(closure.νy), νz=$(closure.νz)), " *
              "(κx=$(closure.κx), κy=$(closure.κy), κz=$(closure.κz))")


text:
using IOExtension, Test

path, io = mktemp()
println(io, "ABC")
println(io, "XYZ")
close(io)

path_out = tempname()

# open and work with various files
openmultiple((path, path_out), ("r", "w")) do io_in, io_out
    while !eof(io_in)
        s = readline(io_in)
        println(io_out, "copy: " * s)
    end
end

open(path_out) do io
    @test readline(io) == "copy: ABC"
    @test readline(io) == "copy: XYZ"
    @test eof(io)
end
rm(path_out)


text:
using NewsvendorModel, Distributions
using Test, Documenter

@testset "Testing NVModel initialization" begin include("nvmodel_test.jl") end

@testset "Testing metrics" begin include("newsvendorproblem_test.jl") end

doctest(NewsvendorModel)


text:
using Base.Test
using Compat
import IJulia: Msg

idents = ["idents"]
header = Dict("msg_id"=>"c673eed8-7c36-47f4-82af-df8ec546a87d",
              "msg_type"=>"comm_msg",
              "username"=>"username",
              "version"=>"5.0",
              "session"=>"980FF7B5F8A24ECF9876A91F38F4EC09")
content = Dict("comm_id"=>"d733be5d-fd21-4194-95da-db364346f2c9",
               "data"=>Dict(:value=>Dict("text/plain"=>"10")))
msg = Msg(idents, header, content)

@test 1 == length(msg.idents)
@test "c673eed8-7c36-47f4-82af-df8ec546a87d" == msg.header["msg_id"]


text:
using Wilsonloop
using Test

function test()

    println("plaq")
    plaq = make_plaq()
    display(plaq)
    for μ=1:4
        println("μ = $μ")
        staples = make_plaq_staple(μ)
        display(staples)
    end
    show(loops_staple_prime[(4,4)])

    #return

    loop = [(1,+1)]

    println(loop)
    w = Wilsonline(loop)
    loop2 = [(1,-1)]
    w2 = Wilsonline(loop2)

    V1 = derive_U(w,1)
    V2 = derive_U(w2,1)
    show(V1)
    println("d")
    show(V2)
    println("d")

    loop = [(1,+1),(2,+1),(1,-1),(2,-1)]
    println(loop)
    w = Wilsonline(loop)
    println("P: ")
    show(w)
    println("P^+: ")
    show(w')
    println("staple")
    for μ=1:4
        println("μ = $μ")
        V1 = make_staple(w,μ)
        V2 = make_staple(w',μ)
        show(V1)
        show(V2)
    end

    println("derive w")
    for μ=1:4
        dU = derive_U(w,μ)
        for i=1:length(dU)
            show(dU[i])
        end
    end

    println("-------------------------------------------------------")
    println("C and dC/dU")
    for μ=1:4
        C = make_Cμ(w,μ)
        #=
        V1 = make_staple(w,μ)
        V2 = make_staple(w',μ)
        C = eltype(V1)[]
        for i=1:length(V1)
            push!(C,V1[i]')
        end
        for i=1:length(V2)
            push!(C,V2[i]')
        end
        =#
        println("-------------------------------------------")
        println("μ = $μ")
        for i=1:length(C)
            println("---------------------------------------")
            println("C[$i]: ")
            show(C[i])
            for ν=1:4
                println("-----------------------------")
                println("ν = $ν")
                dCdU = derive_U(C[i],ν)
                println("dC_{$μ}/dU_{$ν}: ")
                for j=1:length(dCdU)
                    show(dCdU[j])
                end
            end
        end
    end

    println("-------------------------------------------------------")
    println("C and dC/dUdag")
    for μ=1:4
        C = make_Cμ(w,μ)
        #=
        V1 = make_staple(w,μ)
        V2 = make_staple(w',μ)
        C = eltype(V1)[]
        for i=1:length(V1)
            push!(C,V1[i]')
        end
        for i=1:length(V2)
            push!(C,V2[i]')
        end
        =#
        println("-------------------------------------------")
        println("μ = $μ")
        for i=1:length(C)
            println("---------------------------------------")
            println("C[$i]: ")
            show(C[i])
            for ν=1:4
                println("-----------------------------")
                println("ν = $ν")
                dCdU = derive_Udag(C[i],ν)
                println("dC_{$μ}/dUdag_{$ν}: ")
                for j=1:length(dCdU)
                    show(dCdU[j])
                end
            end
        end
    end

    #=
    w = Wilsonline(loop)
    println("P: ")
    display(w)
    println("P^+: ")
    display(w')
    println("staple")
    make_staples(w)
    =#
end

@testset "Wilsonloop.jl" begin
    test()
    @test true
    # Write your tests here.
end


text:
__precompile__()

module FEDVR

using SparseArrays
using RecipesBase

include("grid.jl")
include("lagrange.jl")
include("basis.jl")
include("derivatives.jl")
include("projections.jl")

end # module


text:
# Removing t = 0, such that Σ is invertible
t = Vector(0.1:0.1:10)
p = 2

# Creating generators U,V that result in a positive-definite matrix Σ
Ut, Vt = spline_kernel(t', p)
K = SymSemiseparableMatrix(Ut,Vt)
Σ    = Matrix(K)
chol = cholesky(Σ)

# Creating a symmetric extended generator representable semiseperable matrix
# Calculating its Cholesky factorization
L = cholesky(K)
# Creating a test vector
xt = ones(size(K,1),10)

# Testing size
@test size(L,1) == length(t)
@test size(L,2) == length(t)

# Testing inverses (Using Cholesky factorizations)
@test chol.L\xt  ≈ L\xt
@test chol.U\xt  ≈ L'\xt
@test chol.L*xt  ≈ L*xt
@test chol.L'*xt ≈ L'*xt

# Testing logdet and det
@test logdet(L) ≈ logdet(chol.L)
@test det(L) ≈ det(chol.L)

# Testing show
@test L.L ≈ tril(Ut'*L.Wt)
@test L.U ≈ triu(L.Wt'*Ut)
@test Matrix(L) ≈ tril(Ut'*L.Wt)
@test L[3,1] ≈ chol.L[3,1]
@test L[2,2] ≈ chol.L[2,2]
@test L[1,3] ≈ chol.L[1,3]


text:
"""
    volatility(returns; multiplier=1.0)

Calculates the volatility based on the standard deviation of the returns. The optional `multiplier` parameter allows for scaling the resulting volatility metric, i.e. for annualization.

# Formula

    Vol = std(returns) * sqrt(multiplier)

# Arguments
- `returns`:    Vector of asset returns (usually log-returns).
- `multiplier`: Optional scalar multiplier, i.e. use `12` to annualize monthly returns, and use `252` to annualize daily returns.
"""
function volatility(returns; multiplier=1.0)
    std(returns) * sqrt(multiplier)
end


text:
immutable AngleSpace end

rand(rng::AbstractRNG, ::AngleSpace) = rand(rng)*2*pi
actions(::PowseekerMDP) = AngleSpace()

immutable GPSOrAngleSpace
    p_gps::Float64
end
actions(::PowseekerPOMDP) = GPSOrAngleSpace(0.1)

function rand(rng::AbstractRNG, s::GPSOrAngleSpace)
    if rand(rng) < s.p_gps
        return GPSOrAngle(true, 0.0)
    else
        return GPSOrAngle(false, rand(rng, AngleSpace()))
    end
end


text:
# ### scatter plot 3d

#md # ![FILE_NAME.png](images/FILE_NAME.png)

using Random
Random.seed!(2018)

n = 50
x = rand(n)
y = rand(n)
z = rand(n)
ms = rand(50) * 30
scatter(x, y, z, markersize=ms)


text:
function change_side(side)
    if side == Nord
        return West
    end
    if side == West
        return Sud
    end
    if side == Sud
        return Ost
    end
    if side == Ost
        return Nord
    end
end

function rev(side)
    if side == Nord
        return Sud
    end
    if side == Sud
        return Nord
    end
    if side == West
        return Ost
    end
    if side == Ost
        return West
    end
end

function move_count(r,side,count)
    for i in 1:count
        if isborder(r, side) == false
            move!(r,side)
        else
            f, c = obhod(r, side)
            if f == false
                break
            elseif i + c > count - 1
                break
            end
        end
        if ismarker(r) == true
            break
        end
    end
end

function main_function19(r)
    count = 1
    side = Nord
    while ismarker(r)==false
        move_count(r,side,count)
        side = change_side(side)
        move_count(r,side,count)
        count += 1
        side = change_side(side)
    end
end
function change_side(side)
    if side == Nord
        return West
    end
    if side == West
        return Sud
    end
    if side == Sud
        return Ost
    end
    if side == Ost
        return Nord
    end
end

function rev(side)
    if side == Nord
        return Sud
    end
    if side == Sud
        return Nord
    end
    if side == West
        return Ost
    end
    if side == Ost
        return West
    end
end

function obhod(r,side)
    path = []
    count = 0
    count1 = 0
    f = true
    while isborder(r,side)==true
        if isborder(r,change_side(side))==true
            f = false
            break
        end
        move!(r,change_side(side))
        push!(path,rev(change_side(side)))
        count += 1
    end
    if f == true
        move!(r,side)
        push!(path,rev(side))
        while isborder(r,rev(change_side(side)))
            move!(r,side)
            push!(path,rev(side))
            count1 += 1
        end
        for i in 1:count
            move!(r,rev(change_side(side)))
            push!(path,change_side(side))
        end
    end
    return f, count1
end


text:
# convert folders in results into single csvs
# if the last line is not a result line => ignored
using DataFrames, CSV

function create_csv()
    root_dir = "graph_color/results"
    max_time = 1800

    for dir in readdir(root_dir)
        !isdir(joinpath(root_dir, dir)) && continue
        df = DataFrame(
            "instance" => String[],
            "status" => String[],
            "result" => Float64[],
            "time" => Float64[],
        )
        for (root, dirs, files) in walkdir(joinpath(root_dir, dir))
            pnames = joinpath.(root, files) # files is a Vector{String}, can be empty
            for pname in pnames
                parts = split(pname, "/")
                if parts[end] == "stdout"
                    instance = parts[end - 1]
                    lines = readlines(pname)
                    isempty(lines) && continue
                    result = lines[end]
                    result_parts = split(result, ", ")
                    length(result_parts) != 3 && continue
                    status, result_str, time_str = result_parts
                    result = parse(Float64, result_str)
                    time = parse(Float64, time_str)
                    push!(df, (instance, status, result, time))
                end
            end
            CSV.write(joinpath(root_dir, "$dir.csv"), df)
        end
    end
end

create_csv()


text:
struct PSISettings
    horizon::Base.RefValue{Int}
    use_forecast_data::Bool
    use_parameters::Bool
    warm_start::Base.RefValue{Bool}
    balance_slack_variables::Bool
    services_slack_variables::Bool
    initial_time::Base.RefValue{Dates.DateTime}
    PTDF::Union{Nothing, PSY.PTDF}
    optimizer::Union{Nothing, JuMP.MOI.OptimizerWithAttributes}
    constraint_duals::Vector{Symbol}
    system_to_file::Bool
    export_pwl_vars::Bool
    allow_fails::Bool
    ext::Dict{String, Any}
end

function PSISettings(
    sys;
    initial_time::Dates.DateTime = UNSET_INI_TIME,
    use_parameters::Bool = false,
    use_forecast_data::Bool = true,
    warm_start::Bool = true,
    balance_slack_variables::Bool = false,
    services_slack_variables::Bool = false,
    horizon::Int = UNSET_HORIZON,
    PTDF::Union{Nothing, PSY.PTDF} = nothing,
    optimizer::Union{Nothing, JuMP.MOI.OptimizerWithAttributes} = nothing,
    constraint_duals::Vector{Symbol} = Vector{Symbol}(),
    system_to_file = true,
    export_pwl_vars = false,
    allow_fails = false,
    ext = Dict{String, Any}(),
)
    return PSISettings(
        Ref(horizon),
        use_forecast_data,
        use_parameters,
        Ref(warm_start),
        balance_slack_variables,
        services_slack_variables,
        Ref(initial_time),
        PTDF,
        optimizer,
        constraint_duals,
        system_to_file,
        export_pwl_vars,
        allow_fails,
        ext,
    )
end

function copy_for_serialization(settings::PSISettings)
    vals = []
    for name in fieldnames(PSISettings)
        if name == :optimizer
            # Cannot guarantee that the optimizer can be serialized.
            val = nothing
        else
            val = getfield(settings, name)
        end

        push!(vals, val)
    end

    return deepcopy(PSISettings(vals...))
end

function restore_from_copy(
    settings::PSISettings;
    optimizer::Union{Nothing, JuMP.MOI.OptimizerWithAttributes},
)
    vals = []
    for name in fieldnames(PSISettings)
        if name == :optimizer
            val = optimizer
        else
            val = getfield(settings, name)
        end

        push!(vals, val)
    end

    return PSISettings(vals...)
end

function set_horizon!(settings::PSISettings, horizon::Int)
    settings.horizon[] = horizon
    return
end
get_horizon(settings::PSISettings)::Int = settings.horizon[]
get_use_forecast_data(settings::PSISettings) = settings.use_forecast_data
get_use_parameters(settings::PSISettings) = settings.use_parameters
function set_initial_time!(settings::PSISettings, initial_time::Dates.DateTime)
    settings.initial_time[] = initial_time
    return
end
get_initial_time(settings::PSISettings)::Dates.DateTime = settings.initial_time[]
get_PTDF(settings::PSISettings) = settings.PTDF
get_optimizer(settings::PSISettings) = settings.optimizer
get_ext(settings::PSISettings) = settings.ext
function set_warm_start!(settings::PSISettings, warm_start::Bool)
    settings.warm_start[] = warm_start
    return
end
get_warm_start(settings::PSISettings) = settings.warm_start[]
get_constraint_duals(settings::PSISettings) = settings.constraint_duals
get_balance_slack_variables(settings::PSISettings) = settings.balance_slack_variables
get_services_slack_variables(settings::PSISettings) = settings.services_slack_variables
get_system_to_file(settings::PSISettings) = settings.system_to_file
get_export_pwl_vars(settings::PSISettings) = settings.export_pwl_vars


text:
using LinearAlgebra

@testset "PointForce " begin

  @testset "Gaussian function" begin
    @test Gaussian!(1.0, 0.1) ≈ exp(-0.5*1.0^2*0.1^(-2))
    @test Gaussian(2.0)(0.5) ≈ Gaussian!(0.5, 2.0)
  end

  @testset "Spatial Gaussian distribution for NodeData" begin

    Re = 500; # Reynolds number
    xlim = (-1.0,5.0)
    ylim = (-2.0,2.0)
    ReΔ = 16
    Δx = ReΔ/Re
    Δt = min(0.25*Δx,0.5*Δx^2*Re);
    sys = Systems.NavierStokes(Re, Δx, xlim, ylim, Δt)

    w = Nodes(Dual,size(sys));
    wbuffer =  Nodes(Dual,size(sys));
    x0 = (1.0, 2.0)
    σx = 2.0

    SpatialGauss(wbuffer, w, x0, σx, sys)
    xx, yy = coordinates(w,dx=sys.grid.Δx,I0=Systems.origin(sys))
    @test wbuffer[10, 20] ≈ exp(-0.5*σx^(-2)*((xx[10]-x0[1])^2+(yy[20]-x0[2])^2))
  end

  @testset "Spatial Gaussian distribution for EdgeData" begin

    Re = 500; # Reynolds number
    xlim = (-1.0,5.0)
    ylim = (-2.0,2.0)
    ReΔ = 16
    Δx = ReΔ/Re
    Δt = min(0.25*Δx,0.5*Δx^2*Re);
    sys = Systems.NavierStokes(Re, Δx, xlim, ylim, Δt)

    q = Edges(Primal,size(sys));
    qbuffer =  Edges(Primal,size(sys));
    x0 = (1.0, 2.0)
    σx = 2.0

    SpatialGauss(qbuffer, q, x0, σx, sys)
    xx, xy, yx, yy = coordinates(q,dx=sys.grid.Δx,I0=Systems.origin(sys))
    @test qbuffer.u[10, 20] ≈ exp(-0.5*σx^(-2)*((xx[10]-x0[1])^2+(xy[20]-x0[2])^2))

    @test qbuffer.v[60, 41] ≈ exp(-0.5*σx^(-2)*((yx[60]-x0[1])^2+(yy[41]-x0[2])^2))

  end

@testset "PointForce(t) without spatial variation" begin
  Re = 500; # Reynolds number
  xlim = (-1.0,5.0)
  ylim = (-2.0,2.0)
  ReΔ = 16
  Δx = ReΔ/Re
  Δt = min(0.25*Δx,0.5*Δx^2*Re);
  sys = Systems.NavierStokes(Re, Δx, xlim, ylim, Δt)

  q = Edges(Primal,size(sys));
  qbuffer =  Edges(Primal,size(sys));
  xx, xy, yx, yy = coordinates(q,dx=sys.grid.Δx,I0=Systems.origin(sys))
  x0 = (xx[10], xy[30])
  f0 = (1.0, 2.0)
  σx = 2.0
  t0 = 2.0
  σt = 0.5
  pointf = PointForce(q, x0, f0, t0, σt, sys)

  # Test at the central location and time of the Gaussian distribution
  @test pointf(2.0).u[10, 30] ≈ (pointf.regop(pointf.ubuffer,deepcopy(pointf.fbuffer))).u[10, 30]

  # Test at the central time of the Gaussian distribution
  @test pointf(2.0).v.data  ≈  pointf.regop(pointf.ubuffer,deepcopy(pointf.fbuffer)).v.data

  @test pointf(2.5).u.data  ≈  pointf.regop(pointf.ubuffer,rmul!(deepcopy(pointf.fbuffer),
                              exp(-(2.5-pointf.t0)^2/pointf.σt^2))).u.data

  @test pointf(2.1).v.data  ≈  pointf.regop(pointf.ubuffer,rmul!(deepcopy(pointf.fbuffer),
                              exp(-(2.1-pointf.t0)^2/pointf.σt^2))).v.data

end

@testset "PointForce(t) with spatial variation" begin
        Re = 500; # Reynolds number
        xlim = (-1.0,5.0)
        ylim = (-2.0,2.0)
        ReΔ = 16
        Δx = ReΔ/Re
        Δt = min(0.25*Δx,0.5*Δx^2*Re);
        sys = Systems.NavierStokes(Re, Δx, xlim, ylim, Δt)

        q = Edges(Primal,size(sys));
        qbuffer =  Edges(Primal,size(sys));
        xx, xy, yx, yy = coordinates(q,dx=sys.grid.Δx,I0=Systems.origin(sys))
        x0 = (xx[10], xy[30])
        f0 = (1.0, 2.0)
        σx = 2.0
        t0 = 2.0
        σt = 0.5
        pointf = PointForce(q, x0, f0, t0, σt, sys, σx = σx)

        # Test at the central location and time of the Gaussian distribution
        @test pointf(2.0).u[10, 30] ≈  (pointf.regop(pointf.ubuffer,deepcopy(pointf.fbuffer))).u[10, 30]

        # Test at the central time of the Gaussian distribution
        @test pointf(2.0).v.data ≈  (pointf.xbuffer∘pointf.regop(pointf.ubuffer,deepcopy(pointf.fbuffer))).v.data

        @test pointf(2.5).u.data ≈  (pointf.xbuffer∘pointf.regop(pointf.ubuffer,rmul!(deepcopy(pointf.fbuffer),
                                    exp(-(2.5-pointf.t0)^2/pointf.σt^2)))).u.data

        @test pointf(2.1).v.data ≈  (pointf.xbuffer∘pointf.regop(pointf.ubuffer,rmul!(deepcopy(pointf.fbuffer),
                                    exp(-(2.1-pointf.t0)^2/pointf.σt^2)))).v.data

  end
end


text:
validate_buoyancy(::Nothing, tracers) = nothing

required_tracers(::Nothing) = ()

@inline buoyancy_perturbation(i, j, k, grid::AbstractGrid{FT}, ::Nothing, C) where FT = zero(FT)

@inline ∂x_b(i, j, k, grid::AbstractGrid{FT}, ::Nothing, C) where FT = zero(FT)
@inline ∂y_b(i, j, k, grid::AbstractGrid{FT}, ::Nothing, C) where FT = zero(FT)
@inline ∂z_b(i, j, k, grid::AbstractGrid{FT}, ::Nothing, C) where FT = zero(FT)

@inline x_dot_g_b(i, j, k, grid, ::Nothing, C) = 0
@inline y_dot_g_b(i, j, k, grid, ::Nothing, C) = 0
@inline z_dot_g_b(i, j, k, grid, ::Nothing, C) = 0


text:
using DataStructure101
const DS = DataStructure101

export bfs

function bfs(g::DS.AbstractGraph, v::Int64)
    traversed = Int[]
    queue = Int[]
    push!(queue, v)

    while length(queue) > 0
        current_v = queue[1]
        push!(traversed, popfirst!(queue))

        neighbors = DS.neighbor(g, current_v)
        for n in neighbors
            if !(n in traversed || n in queue) push!(queue, n) end
        end
    end

    return traversed
end


text:
abstract type AbstractIFParameter end
@with_kw struct IFParameter{VF} <: AbstractIFParameter
    τm::VF = 20ms
    τe::VF = 5ms
    τi::VF = 10ms
    Vt::VF = -50mV
    Vr::VF = -60mV
    El::VF = Vr
end

abstract type AbstractIF end
@with_kw mutable struct IF{VF,VB} <: AbstractIF
    param::IFParameter = IFParameter()
    N::SNNInt = 100
    v::VF = param.Vr .+ rand(eltype(VF), N) .* (param.Vt - param.Vr)
    ge::VF = zeros(eltype(VF), N)
    gi::VF = zeros(eltype(VF), N)
    fire::VB = zeros(eltype(VB), N)
    I::VF = zeros(eltype(VF), N)
    records::Dict = Dict()
end
IF(x;kwargs...) = IF{Vector{SNNFloat},Vector{Bool}}(;kwargs...)

function integrate!(p::IF, param::IFParameter, dt::SNNFloat)
    @unpack N, v, ge, gi, fire, I = p
    @unpack τm, τe, τi, Vt, Vr, El = param
    @inbounds for i = 1:N
        v[i] += dt * (ge[i] + gi[i] - (v[i] - El) + I[i]) / τm
        ge[i] += dt * -ge[i] / τe
        gi[i] += dt * -gi[i] / τi
    end
    @inbounds for i = 1:N
        fire[i] = v[i] > Vt
        v[i] = ifelse(fire[i], Vr, v[i])
    end
end


text:
using BEAST
using CompScienceMeshes
using Base.Test

m = readmesh(joinpath(dirname(@__FILE__),"sphere.in"))
X = raviartthomas(m)
κ = ω = c = 1.0

fn = Pkg.dir("BEAST","test","efie_solution.txt")
u = map(x->eval(parse(x)), readcsv(fn))

fcr = facecurrents(u,X)

T, P = linspace(0.0,π,7), 0.0
points = vec([point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for θ in T, ϕ in P])
utheta = vec([point(cos(ϕ)*cos(θ), sin(ϕ)*cos(θ), -sin(θ)) for θ in T, ϕ in P])
uphi   = vec([point(-sin(ϕ)*sin(θ), cos(ϕ)*sin(θ), 0) for θ in T, ϕ in P])
for i in eachindex(uphi)
    uphi[i] /= norm(uphi[i])
end

FF = BEAST.MWFarField3D(1.0im)
ff = potential(FF, points, u, X)
ff *= im*ω/(4π*c)

a = Float64[real(norm(f)) for f in ff]
b = Float64[sqrt(abs(dot(ff[i],utheta[i]))^2 + abs(dot(ff[i],uphi[i]))^2) for i in eachindex(ff)]
c = [dot(uphi[i], uphi[i]) for i in eachindex(ff)]

tol = eps() * 10000
@test a[1] ≈ 0.6426645162880508  atol=tol
@test a[2] ≈ 0.5222644345134998  atol=tol
@test a[3] ≈ 0.28287430186444407 atol=tol
@test a[4] ≈ 0.3873439919314161  atol=tol
@test a[5] ≈ 0.6788129041674158  atol=tol
@test a[6] ≈ 0.8800967217249486  atol=tol
@test a[7] ≈ 0.9487895695366841  atol=tol

# NF = BEAST.MWSingleLayerField3D(κ)
# xs = linspace(-8.0, 8.0, 100)
# ys = linspace(-8.0, 8.0, 100)
# points = [point(x,y,0) for x in xs, y in ys]
# nf = BEAST.potential(NF, points, u, X)
#
# pw = PlaneWaveMW(point(0,0,1),point(1,0,0),1.0,complex(1.0))
# for i in eachindex(points); nf[i] -= pw(points[i]); end
# a = Float64[abs(norm(f[1])) for f in nf]
# b = reshape(a,size(points)...)


text:
export pack, unpack

function pack(x::Var, dims::Vector{Int}, padding)
    @assert sum(dims) == size(x,ndims(x))
    s = Base.setindex(size(x), maximum(dims), ndims(x))
    ydata = similar(x.data, s..., length(dims))
    fill!(ydata, padding)

    xst = stride(x.data, ndims(x))
    yst = stride(ydata, ndims(ydata))
    xi = 1
    yi = 1
    for d in dims
        n = xst * d
        copyto!(ydata, yi, x.data, xi, n)
        xi += n
        yi += yst
    end
    Var(ydata, ∇pack!, (x,dims))
end
pack(x::Node, dims, padding) = Node(pack, (x,dims,padding))

function ∇pack!(y::Var, x::Var, dims::Vector{Int})
    isnothing(x.grad) && return
    xst = stride(x.data, ndims(x))
    yst = stride(y.data, ndims(y))
    xi = 1
    yi = 1
    for d in dims
        n = xst * d
        addto!(x.grad, xi, y.grad, yi, n)
        xi += n
        yi += yst
    end
end

function pack(x::UniArray{T,N}, dims::Vector{Int}, padding) where {T,N}
    @assert sum(dims) == size(x,N)
    s = Base.setindex(size(x), maximum(dims), N)
    y = similar(x, s..., length(dims))
    fill!(y, padding)

    xst = stride(x, N)
    yst = stride(y, N+1)
    xi = 1
    yi = 1
    for d in dims
        n = xst * d
        copyto!(y, yi, x, xi, n)
        xi += n
        yi += yst
    end
    y
end

function unpack(x::UniArray{T,N}, dims::Vector{Int}) where {T,N}
    @assert length(dims) == size(x,N)
    s = Base.setindex(Base.front(size(x)), sum(dims), N-1)
    y = similar(x, s...)

    xst = stride(x, N)
    yst = stride(y, N-1)
    xi = 1
    yi = 1
    for d in dims
        n = yst * d
        copyto!(y, yi, x, xi, n)
        xi += xst
        yi += n
    end
    y
end


text:
include(joinpath(WATERFALL_DIR,"src","includes.jl"))
include(joinpath(WATERFALL_DIR,"bin","io.jl"))

colorcycle = false

locals = (
    nsample = 3,
    permutation = perm,
    permute = true,
    ncor = true,
    correlate = true,
    colorcycle = false,
    ylabel = "Efficiency (%)",
    vmin = 12.0,
);

fun = Statistics.mean
args = (fun,)
cascade = define_from(Cascade{Data}, df; locals..., kwargs...)

# pos = scale_for(cascade, Geometry, args...; kwargs...)
# # all(length.(pos).==1) && (pos = getindex.(pos,1))

# sgn = sign(cascade)
# sgn[3][2] = -sgn[3][2]

# label = get_label.(collect_data(cascade))
# shape = _define_from(Box, Coloring, pos, sgn; style=:fill, kwargs...)

# ensurearray(x::AbstractArray) = x
# ensurearray(x) = [x;]

# annot = _define_annotation(cascade, Horizontal, Statistics.mean; locals..., kwargs...)
# cas = set_geometry(cascade, Horizontal, Statistics.mean; locals..., kwargs..., style=:stroke, alpha=1.0)
mutable struct Handle
    label::AbstractArray
    series::AbstractArray
end

mutable struct Legend
    handle::Vector{Handle}
end

Legend( ; handle) = Legend(handle)

#
#       TO-DO:
#
# Differentiate between blending/two labels and coloring/one label
# Make this possible: _define_from(Blending, [-1.,+1])

# function _define_from(::Type{Handle}, cascade::Cascade{T}; colorcycle, kwargs...) where T <: Geometry
# end

function _define_from(::Type{Handle}, cascade::Cascade{Data}, fun::Function, args...;
    x0 = WIDTH-100,
    y0 = 3*SEP,
    dx = 2*SEP,
    dy = SEP,
    space = SEP,
    idx = 1,
    kwargs...,
)
    label = _define_annotation(cascade, Horizontal, fun, args...; kwargs..., scale=0.9)
    cas = set_geometry(cascade, Horizontal, fun, args...; kwargs..., style=:stroke, alpha=1.0)

    # Need to update position.
    y0 = y0 + FONTSIZE*(idx-1)
    shape = first(cas.start.shape)
    setproperty!(shape, :position, (Luxor.Point(x0-dx, y0-dy), Luxor.Point(x0+dx,y0+dy)))

    description =_define_from(Label, string(fun), Luxor.Point(x0+dx+space, y0); halign=:left)

    return Handle([label, cas], [shape, description])
end

"""
This function defines the information required for a
"""
function _define_series(cascade::Cascade{T}, fun::Function, args...;
    x0 = WIDTH-100,
    y0 = 3*SEP,
    dx = 2*SEP,
    dy = SEP,
    space = SEP,
    idx = 1,
    kwargs...,
) where T <: Geometry

    y0 = y0 + FONTSIZE*(idx-1)

    shape = first(cascade.start.shape)
    setproperty!(shape, :position, (Luxor.Point(x0-dx, y0+dy), Luxor.Point(x0+dx,y0-dy)))

    label =_define_from(Label, string(fun), Luxor.Point(x0+dx+space, y0); halign=:left)
    return [shape, label]
end

function _define_from(::Type{Legend}, args...; kwargs...)
    return Legend([_define_from(Handle, args...; kwargs...)])
end

Base.push!(legend::Legend, h::Handle) = push!(legend.handle, h)

# Base.push!(legend::Legend, )

xshape = WIDTH-100;
yshape = 3*SEP

blend = Box(
    (Luxor.Point(xshape-SEP, yshape-0.5*SEP), Luxor.Point(xshape+SEP, yshape+0.5*SEP)),
    Blending(
        (Luxor.Point(xshape, yshape-0.5*SEP), Luxor.Point(xshape, yshape+0.5*SEP)),
        _define_gradient(HEX_GAIN, HEX_LOSS),
    ),
    :fill,
)

# # lab

# # Legend()

# series = [cas, annot]
# label = [shape, lab]

# vlims = vlim(cascade; locals..., kwargs...)
# pos = scale_for(calculate(copy(cascade), args...), Horizontal; kwargs...)

# lab = _define_annotation(cascade, Horizontal, args...)

# v = get_value(collect_data(cascade))


text:
module GraphBenchmarks

using BenchmarkTools

export benchmark, GraphGenerator, GraphAlgorithm, GraphBenchmarkSpec, GraphType,
       Graph500, BFS, AbstractBFS, construct, preparebench, picksources,
       preparebenchmark, generate

abstract GraphGenerator
abstract GraphAlgorithm
abstract GraphBenchmarkSpec
abstract GraphType

abstract AbstractBFS <: GraphAlgorithm
immutable BFS <: AbstractBFS end

immutable Graph500 <: GraphBenchmarkSpec end

include("generators/kronecker.jl")

function benchmark{B <: GraphBenchmarkSpec, T <: GraphType, A <: GraphAlgorithm, G <: GraphGenerator}(
        benchmark::B,
        t::Type{T},
        alg::A,
        generator::G
    )
    bench = preparebenchmark(benchmark, t, alg, generator)
    trial = run(bench)
    trial
end

function construct{T <: GraphType}(t::Type{T}, edges::Array)
    error("Construct not defined for $(T)")
end

function preparebench{B <: GraphBenchmarkSpec, A <: GraphAlgorithm, T <: GraphType}(
    benchmark::B,
    alg::A,
    g::T
    )
    error("preparebench not defined for $(typeof(g))")
end

function preparebenchmark{B <: GraphBenchmarkSpec, T <: GraphType, A <: GraphAlgorithm, G <: GraphGenerator}(
        benchmark::B,
        t::Type{T},
        alg::A,
        generator::G
    )
    #Generate the edges to be added based on the generator. All generators should return a 2D Array with
    #the edges.
    edges = generate(generator)

    #construct a graph `g` of type `t` with the edges in `edges`
	g = construct(t, edges)

    #Alternatively, we could have `generate(t, generator)` that does both of these steps together.
    #Allows for the generator to not have to generate the entire edge list
    bench = preparebench(benchmark, alg, g)
    bench
end

function picksources(::Graph500, nv::Int64)
    srand(0)
    sources = rand(1:nv, 64)
end

end


text:
using MacroTools: @q

"""
    @nograd(f...)

The output of the argument functions will be considered as constant
and will not contribute gradients.
"""
macro nograd(ex)
  isexpr(ex, :tuple) || (ex = Expr(:tuple, ex))
  blk = @q begin end
  for f in ex.args
    back = MacroTools.@q _ -> ($__source__; nothing)
    push!(blk.args, :(@inline Zygote._pullback(::Context, ::Core.Typeof($(esc(f))), args...) = $(esc(f))(args...), $back))
  end
  return blk
end

"""
    nograd(x)

Treats `x` as a constant. For example:

```julia-repl
julia> f(x) = sum(nograd(x) + x); gradient(f, randn(5))
([1.0, 1.0, 1.0, 1.0, 1.0],)
```

In the example above, the gradient is an array of ones even though
`f(x)` computes `sum(2x)`, because `nograd(x)` uses the value of `x`
as if it were a constant.
"""
nograd(x) = x
@nograd nograd

macro which(ex)
  @capture(ex, f_(args__)) || error("Zygote.@which f(args...)")
  :(InteractiveUtils.@which adjoint(Context(), $(esc(f)), $(esc.(args)...)))
end


max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
# auto reload edited files
#= @schedule begin =#
#=     sleep(0.1) =#
#=     @eval using Revise =#
#= end =#

max_stars_repo_name: cmey/dotfiles

max_issues_repo_name: cmey/dotfiles

max_forks_repo_name: cmey/dotfiles


text:
###########################################################################
#  Simplified manifold Metropolis adjusted Langevin algorithm (SMMALA)
#
#  Parameters :
#    - driftStep : drift step size (for scaling the jumps)
#    - tuner: for tuning the drift step size
#
###########################################################################

export SMMALA

println("Loading SMMALA(driftStep, tuner) sampler")

###########################################################################
# SMMALA specific tuners
###########################################################################
abstract SMMALATuner <: MCMCTuner

###########################################################################
# SMMALA type
###########################################################################
immutable SMMALA <: MCMCSampler
  driftStep::Float64
  tuner::Union(Nothing, MCMCTuner)

  function SMMALA(s::Real, t::Union(Nothing, MCMCTuner))
    @assert s>0 "SMMALA drift step should be > 0"
    new(s, t)
  end
end

SMMALA(s::Float64=1.0) = SMMALA(s, nothing)
SMMALA(s::MALATuner) = SMMALA(1.0, t)
SMMALA(;scale::Float64=1.0, tuner::Union(Nothing, MCMCTuner)=nothing) = SMMALA(scale, tuner)

###########################################################################
# SMMALA sampler
###########################################################################
function SamplerTask(model::MCMCModel, sampler::SMMALA, runner::MCMCRunner)
  # Define local variables
  local pars, proposedPars, parsMean
  local logTarget, proposedLogTarget
  local grad, proposedGrad
  local H, proposedH
  local G, invG, cholOfInvG, proposedG, proposedInvG
  local firstTerm, proposedFirstTerm
  local probNewGivenOld, probOldGivenNew
  local driftStep

  @assert hasgradient(model) "SMMALA sampler requires model with gradient function"
  @assert hastensor(model) "SMMALA sampler requires model with tensor function"

   # Task reset function
  function reset(resetPars::Vector{Float64})
    pars = copy(resetPars)
    logTarget, grad, G = model.evalallt(pars)
  end
  # Hook inside Task to allow remote resetting
  task_local_storage(:reset, reset)

  # Initialization
  pars = copy(model.init)
  logTarget, grad, G = model.evalallt(pars)
  @assert isfinite(logTarget) "Initial values out of model support, try other values"

  invG = inv(G)
  firstTerm = invG*grad

  if isa(sampler.tuner, EmpMCTuner); tune = EmpiricalMALATune(sampler.driftStep, 0, 0); end

  i = 1
  while true
    if isa(sampler.tuner, EmpMCTuner)
      tune.proposed += 1
      driftStep = tune.driftStep
    else
      driftStep = sampler.driftStep
    end

    # Calculate the drift term
    parsMean = pars+(driftStep/2)*firstTerm

    # Calculate proposed parameters
    cholOfInvG = chol(driftStep*invG)
    proposedPars = parsMean+cholOfInvG'*randn(model.size)

    # Update model based on the proposed parameters
    proposedLogTarget, proposedGrad, proposedG = model.evalallt(proposedPars)

    probNewGivenOld = (-sum(log(diag(cholOfInvG)))
      -(0.5*(parsMean-proposedPars)'*(G/driftStep)*(parsMean-proposedPars))[1])

    proposedInvG = inv(proposedG)
    proposedFirstTerm = proposedInvG*proposedGrad

    # Calculate the drift term
    parsMean = proposedPars+(driftStep/2)*proposedFirstTerm

    probOldGivenNew = (-sum(log(diag(chol(driftStep*eye(model.size)*proposedInvG))))
      -(0.5*(parsMean-pars)'*(proposedG/driftStep)*(parsMean-pars))[1])

    ratio = proposedLogTarget+probOldGivenNew-logTarget-probNewGivenOld

    if ratio > 0 || (ratio > log(rand()))  # i.e. if accepted
      produce(MCMCSample(proposedPars, proposedLogTarget, proposedGrad, pars, logTarget, grad, {"accept" => true}))
      pars, logTarget, grad = copy(proposedPars), copy(proposedLogTarget), copy(proposedGrad)
      G, invG, firstTerm = copy(proposedG), copy(proposedInvG), copy(proposedFirstTerm)
      if isa(sampler.tuner, EmpMCTuner); tune.accepted += 1; end
    else
      produce(MCMCSample(pars, logTarget, grad, pars, logTarget, grad, {"accept" => false}))
    end

    if isa(sampler.tuner, EmpMCTuner) && i<= runner.burnin && mod(i, sampler.tuner.adaptStep) == 0
      adapt!(tune, sampler.tuner)
      reset!(tune)
      if sampler.tuner.verbose
        println("Burn-in teration $i of $(runner.burnin): ", round(100*tune.rate, 2), " % acceptance rate")
      end
    end

    i += 1
  end
end


text:
using JuLIP: JVec
const BAI = BAInvariants
import Base: ==
using LinearAlgebra: dot, norm
using JuLIP: decode_dict

# -------------- IO -------------------

BondAngleDesc(transform::String, cutoff::Union{String, Tuple}) =
         BondAngleDesc(AnalyticTransform(transform), fcut_analyse(cutoff))

Dict(D::BondAngleDesc) = Dict( "__id__"    =>  "BondAngleDesc",
                                "transform" =>  Dict(D.transform),
                                "cutoff"    =>  Dict(D.cutoff) )

BondAngleDesc(D::Dict) = BondAngleDesc( decode_dict(D["transform"]),
                                        decode_dict(D["cutoff"]) )

==(D1::BondAngleDesc, D2::BondAngleDesc) =
      ( (D1.transform == D2.transform) && (D1.cutoff == D2.cutoff) )

Base.convert(::Val{:BondAngleDesc}, D::Dict) = BondAngleDesc(D)

# ------------- Interface Code ---------------

tdegrees(::BondAngleDesc, vN::Val{N}) where {N} = BAI.tdegrees(vN)

@inline ricoords(D::BondAngleDesc, Rs, J) = lengths_and_angles(Rs, J)

@inline gradri2gradR!(desc::BondAngleDesc, dVsite, dV_drθ, Rs, J, rθ) =
   _grad_rθ2pos!(dVsite, dV_drθ, Rs, J, rθ...)

# the cut-off for the bond-angle descriptor depends only on r but not on θ
@inline fcut(D::BondAngleDesc, rθ) = fcut(D.cutoff, rθ[1])

@inline function fcut_d(D::BondAngleDesc, rθ)
   fc, fc_d = fcut_d(D.cutoff, rθ[1])
   return fc, vcat(fc_d, zero(typeof(rθ[2])))
end

@inline skip_simplex(D::BondAngleDesc, rθ) = (maximum(rθ[1]) > cutoff(D.cutoff))

@inline _rθ2x(D::BondAngleDesc, r, θ) = vcat(transform.(Ref(D), r), θ)
@inline _rθ2x_d(D::BondAngleDesc, r, θ::SVector{K}) where {K} =
   vcat(transform_d.(Ref(D), r), @SVector ones(K))

@inline invariants(D::BondAngleDesc, rθ::Tuple) =
   BAI.invariants(_rθ2x(D, rθ...))

@inline invariants(D::BondAngleDesc, rθ::SVector{1}) =
   BAI.invariants(_rθ2x(D, rθ[1], (@SVector Float64[])))

@inline invariants(D::BondAngleDesc, rθ::SVector{3}) =
   BAI.invariants(_rθ2x(D, SVector(rθ[1], rθ[2]), SVector(rθ[3])))

@inline invariants(D::BondAngleDesc, rθ::SVector{6}) =
   BAI.invariants( _rθ2x(D, SVector(rθ[1], rθ[2], rθ[3]),
                            SVector(rθ[4], rθ[5], rθ[6])) )

@inline function invariants_ed(D::BondAngleDesc, rθ)
   x = _rθ2x(D, rθ...)
   I1, I2, DI1, DI2 = BAI.invariants_ed(x)
   x_d = _rθ2x_d(D, rθ...)
   return I1, I2, _sdot(x_d, DI1), _sdot(x_d, DI2)
end

# -------------- Kernel Functions --------------

"""
`_bondangles(Rs::AbstractVector)`

Take collection `Rs::AbstractVector{JVec}` representing a simplex
and convert into an ordered vector of bondlengths and bondangle
```
[ r1, r2, r3], [θ12, θ13, θ23 ]
[ r1, r2, r3, r4], [θ12, θ13, θ14, θ23, θ24, θ34 ]
```
"""
@generated function lengths_and_angles(Rs::AbstractVector{JVec{T}}, J::SVector{K}) where {T,K}
   # note K = N-1; eltype(Rs) == JVec{T}
   code = Expr[]
   idx = 0
   # bond lengths
   str_r = "r = @SVector T[ "
   for n = 1:K
      str_r *= "norm(Rs[J[$n]]), "
   end
   push_str!(code, str_r * "]")
   # bond angles
   str_θ = "θ = @SVector T[ "
   for n = 1:K-1, m = (n+1):K
      idx += 1
      str_θ *= "dot(Rs[J[$n]], Rs[J[$m]]) / (r[$n] * r[$m]), "
   end
   push_str!(code, str_θ * "]")
   # -----
   quote
      $(Expr(:meta, :inline))
      @inbounds $(Expr(:block, code...))
      return r, θ
   end
end

"""
convert ∇V (where ∇ is the gradient w.r.t. bond-lengths and bond-angle) into
forces, i.e., into ∇Vsite (where ∇ is the gradient w.r.t. positions)
"""
@generated function _grad_rθ2pos!(dVsite, dV_drθ, Rs, J::SVector{K, Int}, r, θ) where {K}
   # K is the number of neighbours, i.e. N = K+1 counting also the center atom
   # ------
   code = Expr[]
   idx = 0
   # the first K entries of dV, are the derivatives w.r.t. |Ri|
   for k = 1:K
      idx += 1   # idx == k of course
      push!(code, :( dVsite[J[$k]] += (dV_drθ[$idx]/r[$k]) * Rs[J[$k]] ))
   end
   # the remaining ones are dot(R̂_i,R̂_j) = dot(R_i,R_j)/(r_i r_j)
   for n = 1:K-1, m = (n+1):K
      idx += 1
      push!(code, :( theta = dot(Rs[J[$n]], Rs[J[$m]]) / (r[$n]*r[$m]) ) )
      push!(code, :( dVsite[J[$n]] += (dV_drθ[$idx]/r[$n]) *
                                      (Rs[J[$m]]/r[$m] - theta * Rs[J[$n]]/r[$n]) ))
      push!(code, :( dVsite[J[$m]] += (dV_drθ[$idx]/r[$m]) *
                                      (Rs[J[$n]]/r[$n] - theta * Rs[J[$m]]/r[$m]) ))
   end
   quote
      $(Expr(:meta, :inline))
      @inbounds $(Expr(:block, code...))
      return dVsite
   end
end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
# ---
# title: Generate your cover on the fly
# cover: assets/lena_color_256.png
# description: this demo shows you how to generate cover on the fly
# julia: 1.0.0
# ---

# There're many reasons that you don't want to mannually manage the cover image.
# DemoCards.jl allows you to generate the card cover on the fly for demos written
# in julia.

# Let's do this with a simple example

using TestImages, FileIO, ImageShow
## ImageIO backend such as ImageMagick is required
cover = testimage("lena_color_256")
save("assets/lena_color_256.png", cover)

# You can use this image later with `cover`, `# ![](assets/lena_color_256.png)`
# , or you can directly write it in your frontmatter

cover

# ## Advanced Literate & Documenter tricks
#
# The following tricks are what you'd probably need to work on a real-world demo.
#
# ### Hide asset generation
#
# There are cases that you want to hide the codes related to asset generation, e.g., the
# `save(...)` line. Fortunately, this is doable by combining the syntax of Documenter and Literate.
# For example, the following code can be inserted to whereever you need to generate the assets.
# The [`#src` flag](https://fredrikekre.github.io/Literate.jl/v2/fileformat/#Filtering-Lines) is
# used as a filter to tell Literate that this line is reserved only in the original source codes.
#
# ```julia
# save("assets/cover.png", img_noise) #src #md
# ```
#
# Gif image is also supported via `ImageMagick`:
#
# ```julia
# # --- save covers --- #src #md
# using ImageMagick #src #md
# imgs = map(10:5:50) do k #src #md
#     colorview(RGB, rank_approx.(svdfactors, k)...) #src #md
# end #src #md
# ImageMagick.save("assets/color_separations_svd.gif", cat(imgs...; dims=3); fps=2) #src #md
# ```
#
# ### Hide the output result
#
# It does not look good to show the `0` result in the above code block as it does not provide
# anything that the reader cares. There are two ways to fix this.
#
# The first way is to insert a [`#hide`
# flag](https://juliadocs.github.io/Documenter.jl/stable/man/syntax/#@example-block) to the line you
# don't want the user see. For example, if you write your source file in this way
#
# ```julia
# cover = testimage("lena_color_256")
# save("assets/lena_color_256.png", cover) #hide
# ```
#
# it will show up in the generated HTML page as

cover = testimage("lena_color_256")
save("assets/lena_color_256.png", cover) #hide

# The return value is still `0`. Sometimes it is wanted, sometime it is not. To also hide the return
# value `0`, you could insert a trivial `nothing #hide #md` to work it around.
#
# ```julia
# cover = testimage("lena_color_256")
# save("assets/lena_color_256.png", cover) #hide #md
# nothing #hide #md #md
# ```
#
# generates

cover = testimage("lena_color_256")
save("assets/lena_color_256.png", cover) #hide
nothing #hide #md

# No output at all. You could, of course, insert `cover #hide #md` or others you would like readers
# to see:
#
#
# ```julia
# cover = testimage("lena_color_256")
# save("assets/lena_color_256.png", cover) #hide #md #md
# cover #hide #md #md
# ```
#
# generates

cover = testimage("lena_color_256")
save("assets/lena_color_256.png", cover) #hide #md
cover #hide #md

# The `#md` flag is used to keep a clean `.jl` file provided by the download-julia badge:
# ![](https://img.shields.io/badge/download-julia-brightgreen.svg)

# ### Be careful about the spaces and identation
#
# Both Documenter and Literate has specifications on spaces, so you'd probably hit this issue sooner
# or later.
#
# Let's say, you want to create a `warning` box using the `!!! warning` syntax, and you'd probably
# ends up with something like this:
#
# ```text
# # !!! warning
# #    This is inside the warning box
# #
# # This is out side the warning box
# ```
#
# !!! warning
#    This is inside the warning box
#
# This is out side the warning box
#
# ---
#
# The reason that contents are not shown correctly inside the warning box is that you need to insert
# *five* (not four) spaces to the following lines: `Literate` consumes one space, and `Documenter`
# consumes the other four.
#
# ```diff
#  # !!! warning
# -#    This is inside the warning box
# +#     This is inside the warning box
#  #
#  # This is out side the warning box
# ```
#
# !!! warning
#     This is inside the warning box
#
# This is out side the warning box


text:
const DEFAULT_LEVEL = 1
const DEFAULT_WIDTH = 10

struct EntrySpec
    type::DataType
    fmt::String     # C-style format string
    uid::UInt16     # unique ID, corresponds to the order the entry was added
    level::UInt8      # verbosity level. 0 always prints.  Higher number -> lower priority
    width::UInt8    # column width in characters
    ccrayon::ConditionalCrayon
end
EntrySpec(T::DataType, fmt::String, eid, level=DEFAULT_LEVEL, width=DEFAULT_WIDTH;
    ccrayon=ConditionalCrayon()) = EntrySpec(T, fmt, UInt16(eid), UInt8(level), UInt8(width), ccrayon)

Base.@kwdef mutable struct LoggerOpts
    curlevel::UInt8 = DEFAULT_LEVEL
    freq::UInt16 = 10                 # how often header prints
    _count::UInt16 = 0                # internal counter
    headerstyle::Crayon = crayon"bold blue"
    linechar::Char = '—'
    enable::Bool = true
    autosize::Bool = true             # automatically expand columns
end

"""
    SolverLogger.Logger

A logger designed to print tabulated data that can be updated at any point.
It supports varying verbosity levels, each including all of the information
from previous levels.

# Constructor

    SolverLogger.Logger(io=stdout; opts...)
    SolverLogger.Logger(filename; opts...)

The constructor can take either an `IO` object or a filename, in which case it will
be opened with write permissions, replacing any existing contents.

The keyword arguments `opts` are one of the following:

* `curlevel` Current verbosity level of the solver. A non-negative integer.
* `freq` A non-negative integer specifying how often the header row should be printed.
* `headerstyle` A `Crayon` specifying the style of the header.
* `linechar` A `Char` that is used for the deliminating row underneath the header. Set to `\0` to not print a row below the header.
* `enable` Enabled/disabled state of the logger at construction.

# Typical Usage
The fields to be printed are specified before use via [`setentry`](@ref). Here the
user can specify properties of the field such as the width of the column, format
specifications (such as number of decimal places, numeric format, alignment, etc.),
column index, and even conditional formatting via a [`ConditionalCrayon`](@ref).
Each field is assigned a fixed verbosity level, and is only printed if the current
verbosity level of the logger, set via [`setlevel!`](@ref) is greater than or equal
to the level of the field.

Once all the fields for the logger have been specified, typical usage is via the
[`@log`](@ref) macro:

    @log logger "iter" 1
    @log logger "cost" cost * 10  # supports expressions
    @log logger "alpha" alpha
    @log logger alpha             # shortcut for the previous

All calls to `@log` overwrite any previous data. Data is only stored in the logger
if the field is active at the current verbosity.

To print the log, the easiest is via the [`printlog`](@ref) function, which will
automatically print the header rows for you, at the frequency specified by
`logger.opts.freq`. The period can be reset (printing a header at the next call to
`printlog`) via [`resetcount!`](@ref). For more control, the header and rows can
be printed directly via [`printheader`](@ref) and [`printrow`](@ref).

The logger can be completely reset via [`resetlogger!`](@ref).

# Enabling / Disabling
The logger can be enable/disabled via `SolverLogging.enable` and `SolverLogging.disable`.
This overwrites the verbosity level.

"""
struct Logger
    io::IO
    fmt::Dict{String,EntrySpec}  # Collection of entry specifications. UID for each entry is automatically assigned
    fmtfun::Dict{String,Function}
    idx::Vector{Int16}  # determines column order. idx[id] gives the column for entry with id.
    data::Vector{String}
    crayons::Vector{Crayon}
    defaults::Dict{DataType,String}
    opts::LoggerOpts
end
function Logger(io::IO=Base.stdout; opts...)
    fmt = Dict{String,EntrySpec}()
    fmtfun = Dict{String,Function}()
    idx = UInt16[]
    data = String[]
    crayons = Crayon[]
    defaults = _default_formats()
    Logger(io, fmt, fmtfun, idx, data, crayons, defaults, LoggerOpts(; opts...))
end
Logger(filename::AbstractString; kwargs...) = Logger(open(filename, "w"); kwargs...)

isenabled(log::Logger) = log.opts.enable
enable(log::Logger) = log.opts.enable = true
disable(log::Logger) = log.opts.enable = false

function _default_formats()
    Dict(
        AbstractFloat => "%.2e",
        AbstractString => "%s",
        Integer => "%d",
        Any => "%s"   # default to string printing
    )
end

"""
    empty!(log::Logger)

Clears out all data from logger and restores it to the default configuration.
Users should prefer to use [`resetlogger!`](@ref) which calls this function.
"""
function Base.empty!(log::Logger)
    empty!(log.fmt)
    empty!(log.fmtfun)
    empty!(log.idx)
    empty!(log.data)
    empty!(log.defaults)
    merge!(log.defaults, _default_formats())
    return log
end

function clear!(log::Logger)
    level = getlevel(log)
    for (k,v) in pairs(log.fmt)
        idx = log.idx[v.uid]
        if v.level > level
            log.data[idx] = ""
        else
            log.data[idx] = " "^v.width
        end
    end
end

"""
    resetcount!(logger)

Resets the row counter such that the subsequent call to [`printlog`](@ref) will
print a header row, and start the count from that point.
"""
resetcount!(log::Logger) = log.opts._count = 0

"""
    resetlogger!(logger)

Resets the logger to the default configuration. All current data will be lost, including
all fields and default formats.
"""
resetlogger!(log::Logger) = begin empty!(log); resetcount!(log) end

"""
    _log!(logger, name, val)

Internal method for logging a value with the logger. Users should prefer to use the
[`@log`](@ref) macro. If `name` is a registered field, `val` will be stored in the logger.

Internally, this method converts `val` to a string using the format specifications
and calculates the color using the [`ConditionalCrayon`](@ref) for the entry.
"""
function _log!(log::Logger, name::String, val, op::Symbol=:replace)
    if haskey(log.fmt, name)
        espec = log.fmt[name]
        if espec.level <= log.opts.curlevel
            idx = log.idx[espec.uid]
            fun = log.fmtfun[espec.fmt]
            crayon = espec.ccrayon(val)
            if op == :replace
                log.data[idx] = rpad(log.fmtfun[espec.fmt](val), espec.width)
            elseif op == :append
                if espec.type <: AbstractString
                    data0 = rstrip(log.data[idx])
                    newdata = log.fmtfun[espec.fmt](string(val))
                    if all(isspace, data0)
                        data = newdata
                    elseif data0[end] == '.'
                        data = data0 * " " * newdata
                    else
                        data = data0 * ". " * newdata
                    end
                    log.data[idx] = rpad( data, espec.width)
                else
                    @warn "Cannot append to a non-string entry."
                end
            elseif op == :add
                if espec.type <: Number
                    data0 = parse(espec.type, rstrip(log.data[idx]))
                    log.data[idx] = rpad(log.fmtfun[espec.fmt](val + data0), espec.width)
                else
                    @warn "Cannot add to a non-numeric field. Use :append for strings."
                end
            end
            log.crayons[idx] = crayon
            @debug "Logging $name with value $val at index $idx"
            if length(log.data[idx]) > espec.width
                if log.opts.autosize
                    newwidth = round(Int, length(log.data[idx])*1.5)
                    setentry(log, name, width = newwidth)
                    log.opts._count = 0
                    log.data[idx] = rpad(log.data[idx], newwidth)
                else
                    @warn "Entry for $name ($(log.data[idx])) is longer than field width ($(espec.width)). Alignment may be affected. Try increasing the field width."
                end
            end
        else
            @debug "Not logging $name, level not high enough"
        end
        return nothing
    else
        @debug "Rejecting log for $name with value $val"
    end
end

"""
    getlevel(logger)

Gets the current verbosity level for the logger.
"""
getlevel(logger::Logger) = logger.opts.curlevel

"""
    setlevel!(logger, level)

Set the verbosity level for the logger. High levels prints more information.
Returns the previous verbosity level.
"""
function setlevel!(log::Logger, level)
    prevlvl = log.opts.curlevel
    log.opts.curlevel = level

    # Reset all levels that are no longer active
    for (k,v) in pairs(log.fmt)
        idx = log.idx[v.uid]
        if v.level > level
            log.data[idx] = ""
        else
            log.data[idx] = rpad(log.data[idx], v.width)
        end
    end
    return prevlvl
end

"""
    printheader(logger)

Prints the header row(s) for the logger, including only the entries that are active
at the current verbosity level. The style of the header can be changed via
`logger.opts.headerstyle`, which is a Crayon that applies to the entire header.
The repeated character under the header can be specified via `logger.opts.linechar`.
This value can be set to the null character `\0` if this line should be excluded.
"""
function printheader(log::Logger)
    isenabled(log) || return
    _printheader(log.io, log)
    return nothing
end
function _printheader(io::IOStream, log::Logger)
    header = formheader(log)
    println(io, header)
    if log.opts.linechar != '\0'
        println(io, repeat(log.opts.linechar, length(header)))
    end
end
function _printheader(io::IO, log::Logger)
    header = formheader(log)
    println(log.opts.headerstyle(header))
    if log.opts.linechar != '\0'
        println(log.opts.headerstyle(repeat(log.opts.linechar, length(header))))
    end
end

"""
    formheader(logger)

Outputs the header as a string
"""
function formheader(log::Logger)
    names = fill("", length(log.idx))
    for (k,v) in pairs(log.fmt)
        idx = log.idx[v.uid]
        if v.level <= log.opts.curlevel
            names[idx] = rpad(k, v.width)
        end
    end
    header = ""
    for name in names
        header *= name
    end
    return header
end

"""
    printrow(logger)

Prints the data currently stored in the logger. Any entries with a
[`ConditionalCrayon`](@ref) will be printed in the specified color.
Only prints the data for the current verbosity level.
"""
function printrow(log::Logger)
    isenabled(log) || return
    _printrow(log.io, log)
    log.opts._count += 1
    return nothing
end
function _printrow(io::IOStream, log::Logger)
    for v in log.data
        print(io,v)
    end
    println(io)
    flush(io)
end
function _printrow(io::IO, log::Logger)
    for (c,v) in zip(log.crayons, log.data)
        print(c,v)
    end
    println()
    flush(stdout)
end

"""
    formrow(logger)

Outputs the data for the current verbosity level as a string.
"""
function formrow(log::Logger)
    row = ""
    for v in log.data
        row *= v
    end
    return row
end

"""
    printlog(logger)

Prints the data currently in the logger, automatically printing the header
at the frequency specified by `logger.opts.freq`.

The period of the header can be reset using [`resetcount!`](@ref).
"""
function printlog(log::Logger)
    cnt, freq = log.opts._count, log.opts.freq
    if cnt % freq == 0
        printheader(log)
        resetcount!(log)
    end
    printrow(log)
end

# Gets the index of the field `name`
@inline getidx(log::Logger, name::String) = log.idx[log.fmt[name].uid]

_getdata(log::Logger, name::String) = log.data[getidx(log, name)]


text:
using MomentBasedBayesianEstimators
using ModelsGenerators
using ProgressMeter
using Mamba
using Base.Test
srand(1113)
#=------------------------------------------------------------------------------
IV model with a first stage F of CP/m = 10 using ∑z'(y-xβ) = 0 as condition
------------------------------------------------------------------------------=#
srand(7875)
y, x, z = ModelsGenerators.sim_iv_d01(CP = 25);
y = vec(y)

data = (y, x, z);

X = Array(Float64, size(z)...)
f!(X, theta, data) = broadcast!(*, X, data[3], data[1]-data[2]*theta)
m = MinimumDivergenceProblem(X, zeros(5), div = ReverseKullbackLeibler());

bmm = BayesianMomentModel(f!, flatprior!, m, data)
#=
MCMC simulations
=#
n = 6000
burnin = 2000
sim = Chains(n + burnin, 1, names = ["Θ"]);
theta = RWMVariate([0.0]);

postf(theta) = posterior(bmm, theta)[1];

p = Progress(n, 1)
for i in 1:n + burnin
  rwm!(theta, [.5], postf, proposal = Normal)
  sim[i, :, 1] = [theta;]
  next!(p)
end
sim = sim[burnin+1:end,]
describe(sim)
p = plot(sim);
draw(p, filename="elpost_ivmom.svg")

#=------------------------------------------------------------------------------

------------------------------------------------------------------------------=#
function f!(X, θ, data)
    delta = θ[1:5]
    theta = θ[6]
    X[:] = [data[3].*(data[1]-data[3]*delta.*theta) data[3].*(data[2]-data[3]*delta)]
end

_logpostb(beta) = posterior(bmm, [delta; beta])[1]
_logpostd(delta) = posterior(bmm, [delta; beta])[1]

G = Array(Float64, 100, 10)
m = MinimumDivergenceProblem(G, zeros(10), div = ReverseKullbackLeibler())
bmm = BayesianMomentModel(f!, flatprior!, m, data)

n = 6000
burnin = 2000
sim_lik = Chains(n + burnin, 6, names = ["γ₁", "γ₂", "γ₃", "γ₄", "γ₅", "Θ"])
beta  = AMWGVariate([0.01])
delta = AMWGVariate(ones(5)*0.3)
p = Progress(n, 1)
for i in 1:n + burnin
    amwg!(beta, [.1], _logpostb, batchsize = 50, adapt = (i <= burnin))
    amwg!(delta, .1*ones(5), _logpostb, batchsize = 50, adapt = (i <= burnin))
    sim_lik[i, :, 1] = [delta; beta]
    next!(p)
end
describe(sim)

sim_lik2 = Chains(6000, 2, names = ["||δ||", "Θ"]);
sim_lik2[:,1,1] = mapslices(u -> norm(u), sim_lik.value[2001:8000,1:5,1], 2)
sim_lik2[:,2,1] = sim_lik.value[2001:8000,6,1]

p = plot(sim_lik2)
draw(p, filename="el_likmom.svg")

p = contourplot(sim_lik2, bins = 40);
draw(p, filename="el_likmom_contour.svg")

#=------------------------------------------------------------------------------

------------------------------------------------------------------------------=#

function jeffreyprior!(a, theta)
    a[1] = norm(theta[1:5])*(1+theta[6])^.5
    a
end

bmm = BayesianMomentModel(f!, jeffreyprior!, m, data)
_logpostb(beta) = posterior(bmm, [delta; beta])[1]
_logpostd(delta) = posterior(bmm, [delta; beta])[1]

n = 6000
burnin = 2000
sim_likj = Chains(n + burnin, 6, names = ["γ₁", "γ₂", "γ₃", "γ₄", "γ₅", "Θ"])
beta  = AMWGVariate([0.01])
delta = AMWGVariate(ones(5)*0.3)
p = Progress(n, 1)
for i in 1:n + burnin
    amwg!(beta, [.1], _logpostb, batchsize = 50, adapt = (i <= burnin))
    amwg!(delta, .1*ones(5), _logpostb, batchsize = 50, adapt = (i <= burnin))
    sim_likj[i, :, 1] = [delta; beta]
    next!(p)
end
describe(sim)

sim_likj2 = Chains(6000, 2, names = ["||δ||", "Θ"]);
sim_likj2[:,1,1] = mapslices(u -> norm(u), sim_likj.value[2001:8000,1:5,1], 2)
sim_likj2[:,2,1] = sim_likj.value[2001:8000,6,1]

p = plot(sim_likj2)
draw(p, filename="el_likmom_jeffrey.svg")

p = contourplot(sim_likj2, bins = 40);
draw(p, filename="el_likmom_jeffrey_contour.svg")


text:
import TriangleMesh, UnicodePlots, Plots
using Printf, LoopVectorization, LinearAlgebra, SparseArrays

include("CreateMeshFCFV.jl")
include("VisuFCFV.jl")
include("DiscretisationFCFV.jl")

#--------------------------------------------------------------------#

function SetUpProblem!( mesh, T, Tdir, Tneu, se, a, b, c, d, alp, bet )
    # Evaluate T analytic on cell faces
    @avx for in=1:mesh.nf
        x        = mesh.xf[in]
        y        = mesh.yf[in]
        Tdir[in] = exp(alp*sin(a*x + c*y) + bet*cos(b*x + d*y))
        dTdx     = Tdir[in] * (a*alp*cos(a*x + c*y) - b*bet*sin(b*x + d*y))
        dTdy     = Tdir[in] * (alp*c*cos(a*x + c*y) - bet*d*sin(b*x + d*y))
        Tneu[in] = -dTdy # nx*dTdx + nt*dTdy on SOUTH face
    end
    # Evaluate T analytic on barycentres
    @avx for iel=1:mesh.nel
        x       = mesh.xc[iel]
        y       = mesh.yc[iel]
        T       = exp(alp*sin(a*x + c*y) + bet*cos(b*x + d*y))
        T[iel]  = T
        se[iel] = T*(-a*alp*cos(a*x + c*y) + b*bet*sin(b*x + d*y))*(a*alp*cos(a*x + c*y) - b*bet*sin(b*x + d*y)) + T*(a^2*alp*sin(a*x + c*y) + b^2*bet*cos(b*x + d*y)) + T*(-alp*c*cos(a*x + c*y) + bet*d*sin(b*x + d*y))*(alp*c*cos(a*x + c*y) - bet*d*sin(b*x + d*y)) + T*(alp*c^2*sin(a*x + c*y) + bet*d^2*cos(b*x + d*y))
    end
    return
end

#--------------------------------------------------------------------#

function ComputeError( mesh, Te, qx, qy, a, b, c, d, alp, bet )
    eT  = zeros(mesh.nel)
    eqx = zeros(mesh.nel)
    eqy = zeros(mesh.nel)
    Ta  = zeros(mesh.nel)
    qxa = zeros(mesh.nel)
    qya = zeros(mesh.nel)
    @avx for iel=1:mesh.nel
        x        = mesh.xc[iel]
        y        = mesh.yc[iel]
        Ta[iel]  = exp(alp*sin(a*x + c*y) + bet*cos(b*x + d*y))
        qxa[iel] = -Ta[iel] * (a*alp*cos(a*x + c*y) - b*bet*sin(b*x + d*y))
        qya[iel] = -Ta[iel] * (alp*c*cos(a*x + c*y) - bet*d*sin(b*x + d*y))
        eT[iel]  = Te[iel] - Ta[iel]
        eqx[iel] = qx[iel] - qxa[iel]
        eqy[iel] = qy[iel] - qya[iel]
    end
    errT  = norm(eT)/norm(Ta)
    errqx = norm(eqx)/norm(qxa)
    errqy = norm(eqy)/norm(qya)
    return errT, errqx, errqy
end

#--------------------------------------------------------------------#

function StabParam(tau, dA, Vol, mesh_type)
    if mesh_type=="Quadrangles";        taui = tau; end
    if mesh_type=="UnstructTriangles";  taui = tau; end
    return taui
end

#--------------------------------------------------------------------#

@views function main( N, mesh_type )

    println("\n******** FCFV POISSON ********")

    # Create sides of mesh
    xmin, xmax = 0, 1
    ymin, ymax = 0, 1
    nx, ny     = N, N
    R          = 0.5
    inclusion  = 0

    # Generate mesh
    println("Mesh generation :")
    if mesh_type=="Quadrangles"
        tau  = 1
        @time mesh = MakeQuadMesh( nx, ny, xmin, xmax, ymin, ymax, inclusion, R )
    elseif mesh_type=="UnstructTriangles"
        tau  = 1
        @time  mesh = MakeTriangleMesh( nx, ny, xmin, xmax, ymin, ymax, inclusion, R )
    end
    println("Number of elements: ", mesh.nel, " number of dofs: ", mesh.nf)

    # Source term and BCs etc...
    Tanal  = zeros(mesh.nel)
    se     = zeros(mesh.nel)
    Tdir   = zeros(mesh.nf)
    Tneu   = zeros(mesh.nf)
    alp = 0.1; bet = 0.3; a = 5.1; b = 4.3; c = -6.2; d = 3.4;
    println("Model configuration :")
    @time SetUpProblem!(mesh , Tanal, Tdir, Tneu, se, a, b, c, d, alp, bet)

    # Compute some mesh vectors
    println("Compute FCFV vectors:")
    @time ae, be, ze = ComputeFCFV(mesh, se, Tdir, tau)

    # Assemble element matrices and RHS
    println("Compute element matrices:")
    @time Kv, fv = ElementAssemblyLoop(mesh, ae, be, ze, Tdir, Tneu, tau)

    # Assemble triplets and sparse
    println("Assemble triplets and sparse:")
    @time K, f = CreateTripletsSparse(mesh, Kv, fv)

    # Solve for hybrid variable
    println("Direct solve:")
    # @time Th   = K\f
    PC  = 0.5.*(K.+K')
    t = @elapsed PCc = cholesky(PC)
    @printf("Cholesky took = %02.2e s\n", t)
    Th  = zeros(mesh.nf)
    dTh = zeros(mesh.nf,1)
    r   = zeros(mesh.nf,1)
    r  .= f - K*Th
    # @time Th   = K\f
    for rit=1:5
        r    .= f - K*Th
        println("It. ", rit, " - Norm of residual: ", norm(r)/length(r))
        if norm(r)/length(r) < 1e-10
            break
        end
        dTh  .= PCc\r
        Th  .+= dTh[:]
    end

    # Reconstruct element values
    println("Compute element values:")
    @time Te, qx, qy = ComputeElementValues(mesh, Th, ae, be, ze, Tdir, tau)

    # Compute discretisation errors
    err_T, err_qx, err_qy = ComputeError( mesh, Te, qx, qy, a, b, c, d, alp, bet )
    println("Error in T:  ", err_T )
    println("Error in qx: ", err_qx)
    println("Error in qy: ", err_qy)

    return mesh.nf, err_T, err_qx, err_qy
end

N          = [8, 16, 32, 64 ]#, 128, 256, 512, 1024]
mesh_type  = "Quadrangles"
eT_quad    = zeros(size(N))
eqx_quad   = zeros(size(N))
eqy_quad   = zeros(size(N))
t_quad     = zeros(size(N))
ndof_quad  = zeros(size(N))
for k=1:length(N)
    t_quad[k]    = @elapsed ndof, err_T, err_qx, err_qy = main( N[k], mesh_type )
    eT_quad[k]   = err_T
    eqx_quad[k]  = err_qx
    eqy_quad[k]  = err_qy
    ndof_quad[k] = ndof
end

mesh_type  = "UnstructTriangles"
eT_tri     = zeros(size(N))
eqx_tri    = zeros(size(N))
eqy_tri    = zeros(size(N))
t_tri      = zeros(size(N))
ndof_tri   = zeros(size(N))
for k=1:length(N)
    t_tri[k]     = @elapsed ndof, err_T, err_qx, err_qy = main( N[k], mesh_type )
    eT_tri[k]    = err_T
    eqx_tri[k]   = err_qx
    eqy_tri[k]   = err_qy
    ndof_tri[k]  = ndof
end

p = Plots.plot(  log10.(1.0 ./ N) , log10.(eT_quad), markershape=:rect,      label="Quads"                          )
p = Plots.plot!( log10.(1.0 ./ N) , log10.(eT_tri),  markershape=:dtriangle, label="Triangles", legend=:bottomright, xlabel = "log_10(h_x)", ylabel = "log_10(err_T)" )
# p = Plots.plot(  ndof_quad[2:end], t_quad[2:end], markershape=:rect,      label="Quads"                          )
# p = Plots.plot!( ndof_tri[2:end],  t_tri[2:end],  markershape=:dtriangle, label="Triangles", legend=:bottomright, xlabel = "ndof", ylabel = "time" )
display(p)


text:
function penalty_α(α, pen_α)
    return - 0.5 * sum(α.*α) / pen_α^2
end

function penalty_params(model, pen_params)
    return sum(sum([[penalize(model[j, k], pen_params[j][k]) for k in 1:size(model)[2]] for j in 1:size(model)[1]]))
end

# tmp = [[2. 1. 3. 1.], [2. 3. 4. 5.], [1. 2. 3. 4.]],
#        [[2. 1. 3. 1.], [2. 3. 4. 5.], [1. 2. 3. 5.]]

# [[penalize(model[j, k], tmp[j][k]) for k in 1:size(model)[2]] for j in 1:size(model)[1]]

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
# ====================================================================
# Minimizers
#
abstract type AbstractMinimizer end
abstract type AbstractMinimizerStatus end

# --------------------------------------------------------------------
struct MinimizerStatusOK <: AbstractMinimizerStatus
    best::Vector{Float64}
    unc::Vector{Float64}
    specific
end

struct MinimizerStatusWarn <: AbstractMinimizerStatus
    best::Vector{Float64}
    unc::Vector{Float64}
    message::String
    specific
end

struct MinimizerStatusError <: AbstractMinimizerStatus
    message::String
    specific
end

as_string(::Nothing) = (crayon"bold red", "DRY", "")
function as_string(status::AbstractMinimizerStatus)
    if isa(status, MinimizerStatusOK)
        return (crayon"green", "OK", "")
    elseif isa(status, MinimizerStatusWarn)
        return (crayon"bold yellow", "WARNING!", status.message)
    elseif isa(status, MinimizerStatusWarn)
        return (crayon"bold red", "ERROR!", status.message)
    else
        error("Unsupported type: $(typeof(status))")
    end
end

# --------------------------------------------------------------------
import LsqFit

struct lsqfit <: AbstractMinimizer; end

function minimize(minimizer::lsqfit, func::Function, params::Vector{Parameter})
    ndata = length(func(getfield.(params, :val)))
    res = LsqFit.curve_fit((dummy, pvalues) -> func(pvalues),
                           1.:ndata, fill(0., ndata),
                           getfield.(params, :val),
                           lower=getfield.(params, :low),
                           upper=getfield.(params, :high))

    if !res.converged
        return MinimizerStatusError("Not converged", res)
    end
    return MinimizerStatusOK(getfield.(Ref(res), :param),
                             LsqFit.stderror(res),
                             res)
end

# --------------------------------------------------------------------
import CMPFit

#=
NOTE: using custom thresholds for ftol, gtol and xtol may lead to
unexpected behaviour.  E.g. settings ftol = 1.e-6 may lead to a
non-optimal fit with exit status 2 (possibly because in a single
iteration the improvement is particularly small).

The best approach is probably to use default tolerance values and
either increase the maximum allowed number of iterations
(config.maxiter) or set a threshold for relative fit statistic
improvements (Δfitstat_theshold) to be checked after the minimizer
iterated for the maximum allowed number of times.
=#

mutable struct cmpfit <: AbstractMinimizer
    config::CMPFit.Config
    Δfitstat_theshold::Float64
    cmpfit() = new(CMPFit.Config(), NaN)
end

function minimize(minimizer::cmpfit, func::Function, params::Vector{Parameter})
    guess = getfield.(params, :val)
    low   = getfield.(params, :low)
    high  = getfield.(params, :high)
    parinfo = CMPFit.Parinfo(length(guess))
    for i in 1:length(guess)
        llow  = isfinite(low[i])   ?  1  :  0
        lhigh = isfinite(high[i])  ?  1  :  0
        parinfo[i].limited = (llow, lhigh)
        parinfo[i].limits  = (low[i], high[i])
    end

    last_fitstat = sum(abs2, func(guess))
    while true
        res = CMPFit.cmpfit(func,
                            guess, parinfo=parinfo, config=minimizer.config)

        if res.status <= 0
            return MinimizerStatusError("Status = $(res.status)", res)
        end

        if (res.status == 5)  &&  isfinite(minimizer.Δfitstat_theshold)
            Δfitstat = (last_fitstat - res.bestnorm) / last_fitstat
            if Δfitstat > minimizer.Δfitstat_theshold
                println("\nReached max. number of iteration but relative Δfitstat = $(Δfitstat) > $(minimizer.Δfitstat_theshold), continue minimization...\n")
                last_fitstat = res.bestnorm
                guess = getfield.(Ref(res), :param)
                continue
            end
        end

        if res.status == 2
            return MinimizerStatusWarn(
                getfield.(Ref(res), :param),
                getfield.(Ref(res), :perror),
                "CMPFit status = 2 may imply one (or more) guess values are too far from optimum",
                res)
        end

        return MinimizerStatusOK(getfield.(Ref(res), :param),
                                 getfield.(Ref(res), :perror),
                                 res)
    end
end