text:
using Test

using TerminalClock
using TerminalClock: n2d, COLON_LARGE, setup_timer
using Dates
using TOML

# Taken from Prefernces.jl/test/runtests.jl
function with_temp_depot(f::Function)
    mktempdir() do dir
        saved_depot_path = copy(Base.DEPOT_PATH)
        empty!(Base.DEPOT_PATH)
        push!(Base.DEPOT_PATH, dir)
        try
            f()
        finally
            empty!(Base.DEPOT_PATH)
            append!(Base.DEPOT_PATH, saved_depot_path)
        end
    end
end

# Taken from Prefernces.jl/test/runtests.jl
function activate_and_run(project_dir::String, code::String; env::Dict=Dict())
    mktempdir() do dir
        open(joinpath(dir, "test_code.jl"), "w") do io
            write(io, code)
        end

        out = Pipe()
        cmd =
            setenv(`$(Base.julia_cmd()) --project=$(project_dir) $(dir)/test_code.jl`,
                env..., "JULIA_DEPOT_PATH" => Base.DEPOT_PATH[1])
        p = run(pipeline(cmd, stdout=out, stderr=out); wait=false)
        close(out.in)
        wait(p)
        output = String(read(out))
        if !success(p)
            println(output)
        end
        @test success(p)
        return output
    end
end

@testset "Dial" begin
    for n in 0:9
        @test length(split(n2d(n).str, "\n")) == 9
        @test length(split(n2d(n), "\n")) == 9
    end
    @test length(split(COLON_LARGE, "\n")) == 9
end

@testset "hcat" begin
    d1 = hcat(n2d(1), n2d(2))
    @test length(split(d1, "\n")) == 9

    d2 = hcat(n2d(1), n2d(2), COLON_LARGE, n2d(3), n2d(4))
    @test length(split(d2, "\n")) == 9
end

@testset "clock" begin
    dt = DateTime(2021, 11, 15, 12, 34, 56, 7)
    str = clock(dt)
    txt = joinpath("references", "ASCII", "clock.txt")
    @test str == join(readlines(txt), "\n")
end

@testset "stopwatch" begin
    t = Time(12, 34, 56, 789)
    str = stopwatch(t)
    txt = joinpath("references", "ASCII", "stopwatch.txt")
    @test str == join(readlines(txt), "\n")
end

@testset "setup_timer" begin
    @test setup_timer(hour=1) == Time(1, 0, 0)
    @test setup_timer(minute=2) == Time(0, 2, 0)
    @test setup_timer(second=3) == Time(0, 0, 3)
    @test setup_timer(hour=1, minute=2, second=3) == Time(1, 2, 3)
end

@testset "preference" begin
    local_prefs_toml =
        joinpath(dirname(dirname(pathof(TerminalClock))), "LocalPreferences.toml")
    rm(local_prefs_toml; force=true)
    with_temp_depot() do
        project_dir = joinpath(dirname(@__DIR__))

        # test for set_dials
        activate_and_run(
            project_dir,
            """
    using Pkg; Pkg.instantiate()
    using TerminalClock
    unicodebox = joinpath(dirname(pathof(TerminalClock)), "dials", "UnicodeBox.toml")
    TerminalClock.set_dials(unicodebox)
""",
        )
        prefs = local_prefs_toml |> TOML.parsefile
        @test haskey(prefs, "TerminalClock")
        @test basename(prefs["TerminalClock"]["tomlfile"]) == "UnicodeBox.toml"

        activate_and_run(
            project_dir,
            """
    using Test
    using TerminalClock, Dates
    dt = DateTime(2021, 11, 15, 12, 34, 56, 7)
    str = clock(dt)
    txt = joinpath("references", "UnicodeBox", "clock.txt")
    @test str == join(readlines(txt), "\n")

    t = Time(12, 34, 56, 789)
    str = stopwatch(t)
    txt = joinpath("references", "UnicodeBox", "stopwatch.txt")
    @test str == join(readlines(txt), "\n")
""",
        )

        # test for clear_dials
        activate_and_run(
            project_dir,
            """
    using TerminalClock
    TerminalClock.clear_dials()
""",
        )

        activate_and_run(
            project_dir,
            """
    using Test
    using TerminalClock, Dates
    dt = DateTime(2021, 11, 15, 12, 34, 56, 7)
    str = clock(dt)
    txt = joinpath("references", "ASCII", "clock.txt")
    @test str == join(readlines(txt), "\n")

    t = Time(12, 34, 56, 789)
    str = stopwatch(t)
    txt = joinpath("references", "ASCII", "stopwatch.txt")
    @test str == join(readlines(txt), "\n")
""",
        )
    end # with_temp_depot
end


text:
# Import packages
using Test
using SafeTestsets
using ModelCheck

using Distributions
using Random
using DataFrames
using StatsBase
using GLM
using Gadfly
using Statistics
using Compose
using Colors
using SafeTestsets

# Import simulated data and Poisson general linear model to test functions
include("test_simulated_model.jl")

##############################################
# - Run unit tests
##############################################

@testset "ModelCheck Tests" begin

# Test #1: Test that plot_qq returns an object of type 'Plot' when using Poisson GLM
@test   typeof(plot_qqplot(model = modPoisson,
                           df = df,
                           y = "y")) == Plot

# Test #2: Test that plot_residual_fitted returns an object of type 'Plot' when using Poisson GLM
@test   typeof(plot_residual_fitted(model = modPoisson,
                                    df = df,
                                    y = "y")) == Plot

# Test #3: Test that plot_scale_location returns an object of type 'Plot' when using Poisson GLM
@test   typeof(plot_scale_location(model = modPoisson,
                                   df = df,
                                   y = "y")) == Plot

# Test #4: Test that plot_model_check returns an object of type 'Context' when using Poisson GLM
@test   typeof(plot_model_check(model = modPoisson,
                                df = df,
                                y = "y")) == Context

# Test #5: Test that plot_model_check returns an object of type 'Vector{Float64}' when using Poisson GLM
@test   typeof(calc_deviance_resids(model = modPoisson,
                                    df = df,
                                    y = "y")) == Vector{Float64}

end

ext: jl


text:
"""
pairzip((a, b), (c, d)) gives (a=>c, b=>d)

This function was written by Fengyang Wang and shared on the Julia discourse
forum: https://discourse.julialang.org/t/type-stable-zip-to-pairs/3390/2
"""
pairzip(::Tuple{}, ::Tuple{}) = ()
pairzip(::Tuple{}, ::Tuple) = throw(ArgumentError("args must be equal in length"))
pairzip(::Tuple, ::Tuple{}) = throw(ArgumentError("args must be equal in length"))
pairzip(t::Tuple, u::Tuple) = (t[1] => u[1], pairzip(Base.tail(t), Base.tail(u))...)
pairzip(p::Pair{<:Tuple, <:Tuple}) = pairzip(p.first, p.second)
# inefficient but convenient method to allow subs(p, (x, y)=>[1, 2])
pairzip(p::Pair) = pairzip(Tuple(p.first), Tuple(p.second))

tuplezip(::Tuple{}, ::Tuple{}) = ()
tuplezip(::Tuple{}, ::Tuple) = throw(ArgumentError("args must be equal in length"))
tuplezip(::Tuple, ::Tuple{}) = throw(ArgumentError("args must be equal in length"))
tuplezip(t::Tuple, u::Tuple) = ((t[1], u[1]), tuplezip(Base.tail(t), Base.tail(u))...)

_zip(t::Tuple, u::Tuple) = tuplezip(t, u)
# `tuplezip` would be type unstable
_zip(t::Vector, u::Vector) = zip(t, u)


text:
@doc raw"""
    sigmoid(x)

Elementwise evaluation of the [sigmoid function](https://en.wikipedia.org/wiki/Sigmoid_function)

```math
    \sigma(t) = \frac{1}{1+\exp(t)}
```
"""
sigmoid(x) = @. 1 / (1 + exp(-x))

"""
    allequal(x)::Bool

Inspired by a [Stackoverflow post](https://stackoverflow.com/questions/47564825/check-if-all-the-elements-of-a-julia-array-are-equal/47578613)
"""
@inline function allequal(x)::Bool
    length(x) < 2 && return true
    x1 = x[1]
    @inbounds for i in eachindex(x)
        x[i] == x1 || return false
    end
    return true
end

"""
    normalize_columns!(X::Matrix, p::Real=2)

using `LA.normalize!(x, p)` for each column `x` of `X`.
"""
function normalize_columns!(X::Matrix, p::Real=2)
    for x in eachcol(X)
        LA.normalize!(x, p)
    end
    return X
end

"""
    pairwise_distances(X, Y) = Distances.pairwise(Distances.Euclidean(1e-8), X, Y; dims=2)

Pairwise euclidean distance matrix between columns of `X` and `Y`.
Equivalent to `[norm(x - y) for x in X, y in Y]`.
"""
pairwise_distances(X, Y) = Distances.pairwise(Distances.Euclidean(1e-8), X, Y; dims=2)

"""
    pairwise_distances(X; diag_coeff=Inf) = Distances.pairwise(Distances.Euclidean(1e-8), X; dims=2)

Pairwise euclidean distance matrix between columns of `X`
Equivalent to [`PRS.pairwise_distances`](@ref)`(X, X)` and setting the diagonal elements to `diag_coeff`.
"""
function pairwise_distances(X; diag_coeff=Inf)
    dist = Distances.pairwise(Distances.Euclidean(1e-8), X; dims=2)
    dist[LA.diagind(dist)] .= diag_coeff
    return dist
end

## Graph functions

edgemap(g::LG.AbstractGraph) = Dict(LG.edges(g) .=> 1:LG.ne(g))

sink_nodes(g::LG.SimpleDiGraph) = [v for v in LG.vertices(g) if LG.outdegree(g, v) == 0]

"""
    strong_product(g::G, h::G)::G where {G<:LG.AbstractGraph}

Return the [strong product](https://en.wikipedia.org/wiki/Strong_product_of_graphs) of two graphs

- Section 5 of [BNBLPR19](@cite)
"""
function strong_product(g::G, h::G)::G where {G<:LG.AbstractGraph}
    sp_gh = G(LG.nv(g) * LG.nv(h))
    id(i, j) = (i - 1) * LG.nv(h) + j
    for e1 in LG.edges(g)
        i1, i2 = Tuple(e1)
        for j = 1:LG.nv(h)
            LG.add_edge!(sp_gh, id(i1, j), id(i2, j))
        end
        for e2 in LG.edges(h)
            j1, j2 = Tuple(e2)
            LG.add_edge!(sp_gh, id(i1, j1), id(i2, j2))
            LG.add_edge!(sp_gh, id(i1, j2), id(i2, j1))
        end
    end
    for e in LG.edges(h)
        j1, j2 = Tuple(e)
        for i in LG.vertices(g)
            LG.add_edge!(sp_gh, id(i, j1), id(i, j2))
        end
    end
    return sp_gh
end

"""
    strong_product(g::G)::G where {G<:LG.AbstractGraph}

Return the [strong product](https://en.wikipedia.org/wiki/Strong_product_of_graphs) of a graph with itself

- Section 5 of [BNBLPR19](@cite)
"""
function strong_product(g::G)::G where {G<:LG.AbstractGraph}
    n = LG.nv(g)
    sp_gg = G(n * n)  # nv = n^2, ne = 2n(2n+1)
    id(i, j) = (i - 1) * n + j
    for e1 in LG.edges(g)
        i1, i2 = Tuple(e1)
        for j in LG.vertices(g)
            LG.add_edge!(sp_gg, id(i1, j), id(i2, j))
            LG.add_edge!(sp_gg, id(j, i1), id(j, i2))
        end
        for e2 in LG.edges(g)
            j1, j2 = Tuple(e2)
            LG.add_edge!(sp_gg, id(i1, j1), id(i2, j2))
            LG.add_edge!(sp_gg, id(i1, j2), id(i2, j1))
        end
    end
    return sp_gg
end

king_graph(k::Integer, l::Integer) = strong_product(LG.path_graph(l), LG.path_graph(k))
king_graph(k::Integer) = strong_product(LG.path_graph(k))

"""
    uniform_weighted_graph(
        [rng::Random.AbstractRNG,]
        graph::LG.AbstractGraph
    )::SWG.SimpleWeightedGraph

Return a weighted version of `graph` where each edge is attached an independent uniform random variable.
"""
function uniform_weighted_graph(
    rng::Random.AbstractRNG,
    graph::LG.AbstractGraph
)::SWG.SimpleWeightedGraph
    g = SWG.SimpleWeightedGraph(graph)
    for e in LG.edges(g)
        i, j = Tuple(e)
        @inbounds g.weights[i, j] = g.weights[j, i] = rand(rng, weighttype(g))
    end
    return g
end

uniform_weighted_graph(graph) = uniform_weighted_graph(Random.default_rng(), graph)

"""
    random_edge_orientation(
        rng::Random.AbstractRNG,
        graph::LG.SimpleGraph{T},
        p::Real=0.5
    )::LG.SimpleDiGraph{T} where {T}

Construct an directed version of ``
"""
function random_edge_orientation(
    rng::Random.AbstractRNG,
    graph::LG.SimpleGraph{T},
    p::Real=0.5
)::LG.SimpleDiGraph{T} where {T}
    @assert 0 <= p <= 1
    return LG.SimpleDiGraphFromIterator(flip_edge(rng, e, p) for e in LG.edges(graph))
end

function random_edge_orientation(
    graph::LG.SimpleGraph,
    p::Real=0.5
)
    return random_edge_orientation(Random.default_rng(), graph, p)
end

"""
    flip_edge(
        [rng::Random.AbstractRNG,]
        edge::LG.Edge,
        p::Real=0.5
    )::LG.Edge

return `edge` with probability `p` or `reverse(edge)` with probability `1-p`.
"""
function flip_edge(
    rng::Random.AbstractRNG,
    edge::LG.Edge,
    p::Real=0.5
)::LG.Edge
    return rand(rng) < p ? edge : reverse(edge)
end

flip_edge(edge::LG.Edge, p::Real=0.5) = flip_edge(Random.default_rng(), edge, p)

"""
    random_neighbor_assignment(
        [rng::Random.AbstractRNG,]
        graph::LG.SimpleGraph{T},
        roots=T[]
    )::LG.SimpleDiGraph{T} where {T}

Return a oriented subgraph of `graph` where each vertex except the `roots` is connected to a unique neighbor, i.e., each vertex has outdegree equal to one.
"""
function random_neighbor_assignment(
    rng::Random.AbstractRNG,
    graph::LG.SimpleGraph{T},
    roots=T[]
)::LG.SimpleDiGraph{T} where {T}
    g = LG.SimpleDiGraph(LG.nv(graph))
    for v in LG.vertices(g)
        if v ∉ roots
            w = rand(rng, LG.neighbors(graph, v))
            LG.add_edge!(g, v, w)
        end
    end
    return g
end

function random_neighbor_assignment(
    graph::LG.SimpleGraph{T},
    roots=T[]
) where {T}
    return random_neighbor_assignment(Random.default_rng(), graph, roots)
end

## String methods

function findall_overlap(pat::AbstractString, str::AbstractString)
    found = UnitRange{Int64}[]
    isempty(pat) && return found

    rng = findfirst(pat, str)
    isnothing(rng) && return found

    i, j = first(rng), last(rng)
    k, e = first(rng), lastindex(str)
    while k <= e
        @inbounds k = nextind(str, k)
        rng = findnext(pat, str, k)
        if isnothing(rng)
            push!(found, i:j)
            break
        end
        k = first(rng)
        if j < k
            push!(found, i:j)
            i = k
        end
        j = last(rng)
    end

    return found
end

function has_common_prefix_suffix(s::String)
    n = length(s)
    n <= 1 && return false
    for i in 1:div(n, 2)
        @inbounds s[1:i] == s[n-i+1:n] && return true
    end
    return false
end


text:
@time itercount(100_000_000)
@time itercountcombos(8)
@time itercountcombos(17)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# ---
# title: 25. Reverse Nodes in k-Group
# id: problem25
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Linked List
# link: <https://leetcode.com/problems/reverse-nodes-in-k-group/description/>
# hidden: true
# ---
#
# Given a linked list, reverse the nodes of a linked list _k_ at a time and
# return its modified list.
#
# _k_ is a positive integer and is less than or equal to the length of the
# linked list. If the number of nodes is not a multiple of _k_ then left-out
# nodes, in the end, should remain as it is.
#
# **Follow up:**
#
#   * Could you solve the problem in `O(1)` extra memory space?
#   * You may not alter the values in the list's nodes, only nodes itself may be changed.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/03/reverse_ex1.jpg)
#
#
#
#     Input: head = [1,2,3,4,5], k = 2
#     Output: [2,1,4,3,5]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/03/reverse_ex2.jpg)
#
#
#
#     Input: head = [1,2,3,4,5], k = 3
#     Output: [3,2,1,4,5]
#
#
# **Example 3:**
#
#
#
#     Input: head = [1,2,3,4,5], k = 1
#     Output: [1,2,3,4,5]
#
#
# **Example 4:**
#
#
#
#     Input: head = [1], k = 1
#     Output: [1]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the list is in the range `sz`.
#   * `1 <= sz <= 5000`
#   * `0 <= Node.val <= 1000`
#   * `1 <= k <= sz`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Clang
using Clang.LibClang
using Test

include("cbasic.jl")
include("typedef.jl")
include("struct.jl")
include("wrap_c.jl")
include("macro.jl")
include("union.jl")
# include("generator.jl")


text:
#
# Wrapper for the paralel/structured PIPS-NLP interface
#

module PipsNlpSolver

using StructJuMPSolverInterface
import MPI

try
  sharedLib=ENV["PIPS_NLP_PAR_SHARED_LIB"]

  # Explicitly check if the file exists. dlopen sometimes does not throw an
  # error for invalid filenames, resulting in a seg fault)
  if(!isfile(sharedLib))
    error(string("The specified shared library ([", sharedLib, "]) does not exist"))
  end
  global const libparpipsnlp=Libdl.dlopen(get(ENV,"PIPS_NLP_PAR_SHARED_LIB",""))
catch
  warn("Could not load PIPS-NLP shared library. Make sure the ENV variable 'PIPS_NLP_PAR_SHARED_LIB' points to its location, usually in the PIPS repo at PIPS/build_pips/PIPS-NLP/libparpipsnlp.so")
  rethrow()
end

#######################
type FakeModel <: ModelInterface
    sense::Symbol
    status::Int
    nscen::Int
    rowmap::Dict{Int,Int}
    colmap::Dict{Int,Int}
    eq_rowmap::Dict{Int,Int}
    inq_rowmap::Dict{Int,Int}

    get_num_scen::Function
    get_sense::Function

    get_status::Function
    get_num_rows::Function
    get_num_cols::Function
    get_num_eq_cons::Function
    get_num_ineq_cons::Function

    set_status::Function
    set_num_rows::Function
    set_num_cols::Function
    set_num_eq_cons::Function
    set_num_ineq_cons::Function

    str_init_x0::Function
    str_prob_info::Function
    str_eval_f::Function
    str_eval_g::Function
    str_eval_grad_f::Function
    str_eval_jac_g::Function
    str_eval_h::Function
    str_write_solution::Function

    function FakeModel(sense::Symbol,status::Int,nscen::Int, str_init_x0, str_prob_info, str_eval_f, str_eval_g,str_eval_grad_f,str_eval_jac_g, str_eval_h)
        instance = new(sense,status,nscen,Dict{Int,Int}(),Dict{Int,Int}(),Dict{Int,Int}(),Dict{Int,Int}())
        instance.str_init_x0 = str_init_x0
        instance.str_prob_info = str_prob_info
        instance.str_eval_f = str_eval_f
        instance.str_eval_g = str_eval_g
        instance.str_eval_grad_f = str_eval_grad_f
        instance.str_eval_jac_g = str_eval_jac_g
        instance.str_eval_h = str_eval_h

        instance.str_write_solution = function(id::Integer, x::Vector{Float64}, y_eq::Vector{Float64}, y_ieq::Vector{Float64})
            @show id
            @show x, y_eq, y_ieq
        end

        instance.get_num_scen = function()
            return instance.nscen
        end
        instance.get_sense = function()
            return instance.sense
        end
        instance.get_status = function()
            return instance.status
        end
        instance.get_num_rows = function(id::Integer)
            return instance.rowmap[id]
        end
        instance.get_num_cols = function(id::Integer)
            return instance.colmap[id]
        end
        instance.get_num_eq_cons = function(id::Integer)
            return instance.eq_rowmap[id]
        end
        instance.get_num_ineq_cons = function(id::Integer)
            return instance.inq_rowmap[id]
        end
        instance.set_status = function(s::Integer)
            instance.status = s
        end
        instance.set_num_rows = function(id::Integer, v::Integer)
            return instance.rowmap[id] = v
        end
        instance.set_num_cols = function(id::Integer, v::Integer)
            return instance.colmap[id] = v
        end
        instance.set_num_eq_cons = function(id::Integer, v::Integer)
            return instance.eq_rowmap[id] = v
        end
        instance.set_num_ineq_cons = function(id::Integer, v::Integer)
            return instance.inq_rowmap[id] = v
        end
        return instance
    end
end

#######################

type PipsNlpProblemStruct
    ref::Ptr{Void}
    model::ModelInterface
    comm::MPI.Comm
    prof::Bool

    n_iter::Int
    t_jl_init_x0::Float64
    t_jl_str_prob_info::Float64
    t_jl_eval_f::Float64
    t_jl_eval_g::Float64
    t_jl_eval_grad_f::Float64

    t_jl_eval_jac_g::Float64
    t_jl_str_eval_jac_g::Float64
    t_jl_eval_h::Float64
    t_jl_str_eval_h::Float64
    t_jl_write_solution::Float64

    t_jl_str_total::Float64
    t_jl_eval_total::Float64

    function PipsNlpProblemStruct(comm, model, prof)
        prob = new(C_NULL, model, comm, prof,-3
            ,0.0,0.0,0.0,0.0,0.0
            ,0.0,0.0,0.0,0.0,0.0
            ,0.0,0.0
            )
        finalizer(prob, freeProblemStruct)

        return prob
    end
end

immutable CallBackData
	prob::Ptr{Void}
	row_node_id::Cint
    col_node_id::Cint
    flag::Cint
end

export  ModelInterface, FakeModel,
        createProblemStruct, solveProblemStruct

###########################################################################
# Callback wrappers
###########################################################################

function str_init_x0_wrapper(x0_ptr::Ptr{Float64}, cbd::Ptr{CallBackData})
	# @show " julia - str_init_x0_wrapper "
    # @show cbd
    data = unsafe_load(cbd)
    # @show data
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    # @show prob
    # data = unsafe_pointer_to_objref(cbd)::CallBackData
    # out = Array(Ptr{CallBackData},1)
    rowid = Int(Int(data.row_node_id))
    colid = Int(Int(data.col_node_id))
    assert(rowid == colid)
    n0 = prob.model.get_num_cols(colid)
    x0 = unsafe_wrap(Array, x0_ptr, n0)

    @timing prob.prof tic()
    prob.model.str_init_x0(colid,x0)
    @timing prob.prof prob.t_jl_init_x0 += toq()

    return Int32(1)
end

# prob info (prob_info)
function str_prob_info_wrapper(n_ptr::Ptr{Cint}, col_lb_ptr::Ptr{Float64}, col_ub_ptr::Ptr{Float64}, m_ptr::Ptr{Cint}, row_lb_ptr::Ptr{Float64}, row_ub_ptr::Ptr{Float64}, cbd::Ptr{CallBackData})
    # @show " julia - str_prob_info_wrapper "
    # @show cbd
    data = unsafe_load(cbd)
    # @show data
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    # @show prob
    # data = unsafe_pointer_to_objref(cbd)::CallBackData
    # out = Array(Ptr{CallBackData},1)
    rowid = Int(Int(data.row_node_id))
    colid = Int(data.col_node_id)
    flag = Int(data.flag)
    assert(rowid == colid)

	mode = (col_lb_ptr == C_NULL) ? (:Structure) : (:Values)
    # @show flag
    if flag == 0
        # @show mode
    	if(mode==:Structure)
            col_lb = unsafe_wrap(Array,col_lb_ptr,0)
    		col_ub = unsafe_wrap(Array,col_ub_ptr,0)
    		row_lb = unsafe_wrap(Array,row_lb_ptr,0)
    		row_ub = unsafe_wrap(Array,row_ub_ptr,0)
            @timing prob.prof tic()
    		(n,m) = prob.model.str_prob_info(colid,flag,mode,col_lb,col_ub,row_lb,row_ub)
            @timing prob.prof prob.t_jl_str_prob_info += toq()

    		unsafe_store!(n_ptr,convert(Cint,n)::Cint)
    		unsafe_store!(m_ptr,convert(Cint,m)::Cint)
            # @show typeof(colid), typeof(m)
    		prob.model.set_num_rows(colid, m)
    		prob.model.set_num_cols(colid, n)
    	else
    		n = unsafe_load(n_ptr)
    		m = unsafe_load(m_ptr)
    		col_lb = unsafe_wrap(Array,col_lb_ptr,n)
    		col_ub = unsafe_wrap(Array,col_ub_ptr,n)
    		row_lb = unsafe_wrap(Array,row_lb_ptr,m)
    		row_ub = unsafe_wrap(Array,row_ub_ptr,m)

            @timing prob.prof tic()
    		prob.model.str_prob_info(colid,flag,mode,col_lb,col_ub,row_lb,row_ub)
            @timing prob.prof prob.t_jl_str_prob_info += toq()

    		neq = 0
    		nineq = 0
    		for i = 1:length(row_lb)
    			if row_lb[i] == row_ub[i]
    				neq += 1
    			else
    				nineq += 1
    			end
    		end
    		assert(neq+nineq == length(row_lb) == m)
    		prob.model.set_num_eq_cons(colid,neq)
    		prob.model.set_num_ineq_cons(colid,nineq)
    	end
    else
        @assert flag ==1
        if mode == :Structure
            col_lb = unsafe_wrap(Array,col_lb_ptr,0)
            col_ub = unsafe_wrap(Array,col_ub_ptr,0)
            row_lb = unsafe_wrap(Array,row_lb_ptr,0)
            row_ub = unsafe_wrap(Array,row_ub_ptr,0)
            (n,m) = prob.model.str_prob_info(colid,flag,mode,col_lb,col_ub,row_lb,row_ub)
            # @show n,m
            unsafe_store!(m_ptr,convert(Cint,m)::Cint)
        else
            n = unsafe_load(n_ptr)
            m = unsafe_load(m_ptr)
            @assert m==0
        end
    end
	return Int32(1)
end
# Objective (eval_f)
function str_eval_f_wrapper(x0_ptr::Ptr{Float64}, x1_ptr::Ptr{Float64}, obj_ptr::Ptr{Float64}, cbd::Ptr{CallBackData})
    # @show " julia - eval_f_wrapper "
    data = unsafe_load(cbd)
    # @show data
    # @show data
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    rowid = Int(Int(data.row_node_id))
    colid = Int(data.col_node_id)
    assert(rowid == colid)
    n0 = prob.model.get_num_cols(0)
    n1 = prob.model.get_num_cols(colid)
    # Calculate the new objective
    x0 = unsafe_wrap(Array, x0_ptr, n0)
    x1 = unsafe_wrap(Array, x1_ptr, n1)

    @timing prob.prof tic()
    new_obj = convert(Float64, prob.model.str_eval_f(colid,x0,x1))::Float64
    @timing prob.prof prob.t_jl_eval_f += toq()

    # Fill out the pointer
    unsafe_store!(obj_ptr, new_obj)
    # Done
    return Int32(1)
end

# Constraints (eval_g)
function str_eval_g_wrapper(x0_ptr::Ptr{Float64}, x1_ptr::Ptr{Float64}, eq_g_ptr::Ptr{Float64}, inq_g_ptr::Ptr{Float64}, cbd::Ptr{CallBackData})
    # @show " julia - eval_g_wrapper "
    data = unsafe_load(cbd)
    # @show data
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    rowid = Int(data.row_node_id)
    colid = Int(data.col_node_id)
    assert(rowid == colid)
    n0 = prob.model.get_num_cols(0)
    n1 = prob.model.get_num_cols(colid)
    x0 = unsafe_wrap(Array, x0_ptr, n0)
    x1 = unsafe_wrap(Array, x1_ptr, n1)
    # Calculate the new constraint values
    neq = prob.model.get_num_eq_cons(rowid)
    nineq = prob.model.get_num_ineq_cons(rowid)
    new_eq_g = unsafe_wrap(Array, eq_g_ptr, neq)
    new_inq_g = unsafe_wrap(Array, inq_g_ptr, nineq)

    @timing prob.prof tic()
    prob.model.str_eval_g(colid,x0,x1,new_eq_g,new_inq_g)
    @timing prob.prof prob.t_jl_eval_g += toq()

    # Done
    return Int32(1)
end

# Objective gradient (eval_grad_f)
function str_eval_grad_f_wrapper(x0_ptr::Ptr{Float64}, x1_ptr::Ptr{Float64}, grad_f_ptr::Ptr{Float64}, cbd::Ptr{CallBackData})
    # @show " julia -  eval_grad_f_wrapper "
    # Extract Julia the problem from the pointer
    # @show cbd
    data = unsafe_load(cbd)
    # @show data
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    rowid = Int(data.row_node_id)
    colid = Int(data.col_node_id)
    n0 = prob.model.get_num_cols(0)
    n1 = prob.model.get_num_cols(rowid)
    # @show n0,n1
    x0 = unsafe_wrap(Array, x0_ptr, n0)
    x1 = unsafe_wrap(Array, x1_ptr, n1)
    # Calculate the gradient
    grad_len = prob.model.get_num_cols(colid)
    new_grad_f = unsafe_wrap(Array, grad_f_ptr, grad_len)

    @timing prob.prof tic()
    prob.model.str_eval_grad_f(rowid,colid,x0,x1,new_grad_f)
    @timing prob.prof prob.t_jl_eval_grad_f += toq()

    if prob.model.get_sense() == :Max
        new_grad_f *= -1.0
    end
    # Done
    return Int32(1)
end

# Jacobian (eval_jac_g)
function str_eval_jac_g_wrapper(x0_ptr::Ptr{Float64}, x1_ptr::Ptr{Float64},
	e_nz_ptr::Ptr{Cint}, e_values_ptr::Ptr{Float64}, e_row_ptr::Ptr{Cint}, e_col_ptr::Ptr{Cint},
	i_nz_ptr::Ptr{Cint}, i_values_ptr::Ptr{Float64}, i_row_ptr::Ptr{Cint}, i_col_ptr::Ptr{Cint},
	cbd::Ptr{CallBackData}
	)
    # @show " julia -  eval_jac_g_wrapper "
    # Extract Julia the problem from the pointer
    data = unsafe_load(cbd)
    # @show data
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    rowid = Int(Int(data.row_node_id))
    colid = Int(Int(data.col_node_id))
    flag = Int(data.flag)
    n0 = prob.model.get_num_cols(0)
    n1 = prob.model.get_num_cols(rowid) #we can do this because of 2-level and no linking constraint
    # @show n0, n1
    x0 = unsafe_wrap(Array, x0_ptr, n0)
    x1 = unsafe_wrap(Array, x1_ptr, n1)
    # @show x0
    # @show x1
    nrow = prob.model.get_num_rows(rowid)
    ncol = prob.model.get_num_cols(colid)
    #@show prob
    # Determine mode
    mode = (e_col_ptr == C_NULL && i_col_ptr == C_NULL) ? (:Structure) : (:Values)
    if flag != 1
        if(mode == :Structure)
        	e_values = unsafe_wrap(Array,e_values_ptr,0)
    		e_colptr = unsafe_wrap(Array,e_col_ptr,0)
    		e_rowidx = unsafe_wrap(Array,e_row_ptr,0)
    		i_values = unsafe_wrap(Array,i_values_ptr,0)
    		i_colptr = unsafe_wrap(Array,i_col_ptr,0)
    		i_rowidx = unsafe_wrap(Array,i_row_ptr,0)

            @timing prob.prof tic()
            (e_nz,i_nz) = prob.model.str_eval_jac_g(rowid,colid,flag,x0,x1,mode,e_rowidx,e_colptr,e_values,i_rowidx,i_colptr,i_values)
            @timing prob.prof prob.t_jl_str_eval_jac_g += toq()

            unsafe_store!(e_nz_ptr,convert(Cint,e_nz)::Cint)
    		unsafe_store!(i_nz_ptr,convert(Cint,i_nz)::Cint)
    		# @show "structure - ",(e_nz,i_nz)
        else
        	e_nz = unsafe_load(e_nz_ptr)
        	e_values = unsafe_wrap(Array,e_values_ptr,e_nz)
        	e_rowidx = unsafe_wrap(Array,e_row_ptr,e_nz)
        	e_colptr = unsafe_wrap(Array,e_col_ptr,ncol+1)
        	i_nz = unsafe_load(i_nz_ptr)
        	# @show "values - ",(e_nz,i_nz), ncol
        	i_values = unsafe_wrap(Array,i_values_ptr,i_nz)
        	i_rowidx = unsafe_wrap(Array,i_row_ptr,i_nz)
        	i_colptr = unsafe_wrap(Array,i_col_ptr,ncol+1)
            # @show x0
            # @show x1
            @timing prob.prof tic()
        	prob.model.str_eval_jac_g(rowid,colid,flag,x0,x1,mode,e_rowidx,e_colptr,e_values,i_rowidx,i_colptr,i_values)
            @timing prob.prof prob.t_jl_eval_jac_g += toq()
        end
    else
        @assert flag == 1
        if mode == :Structure
            e_nz = 0
            i_nz = 0
            unsafe_store!(e_nz_ptr,convert(Cint,e_nz)::Cint)
            unsafe_store!(i_nz_ptr,convert(Cint,i_nz)::Cint)
        else
            e_nz = unsafe_load(e_nz_ptr)
            i_nz = unsafe_load(i_nz_ptr)
            @assert e_nz == i_nz == 0
        end
    end
    # Done
    return Int32(1)
end

# Hessian
function str_eval_h_wrapper(x0_ptr::Ptr{Float64}, x1_ptr::Ptr{Float64}, lambda_ptr::Ptr{Float64}, nz_ptr::Ptr{Cint}, values_ptr::Ptr{Float64}, row_ptr::Ptr{Cint}, col_ptr::Ptr{Cint}, cbd::Ptr{CallBackData})
    # @show " julia - eval_h_wrapper "
    # Extract Julia the problem from the pointer
    data = unsafe_load(cbd)
    # @show data
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    rowid = Int(data.row_node_id)
    colid = Int(data.col_node_id)
    flag = Int(data.flag)
    # @message @sprintf(" julia - eval_h_wrapper - %d %d %d", rowid, colid, flag)
    @timing prob.prof begin
        if rowid == colid ==0
            prob.n_iter += 1
            if prob.n_iter == 0
                prob.t_jl_str_total = t_reset(prob)
            end
        end
    end
    # @show prob.n_iter

    high = max(rowid,colid)
    low  = min(rowid,colid)
    n0 = prob.model.get_num_cols(0)
    n1 = prob.model.get_num_cols(high)
    x0 = unsafe_wrap(Array,x0_ptr,n0)
    x1 = unsafe_wrap(Array,x1_ptr,n1)
    # @show x0
    # @show x1
    ncol = prob.model.get_num_cols(low)
    g0 = prob.model.get_num_rows(high)
    # @show g0
    # @show ncol
    lambda = unsafe_wrap(Array,lambda_ptr,g0)
    obj_factor = 1.0
    if prob.model.get_sense() == :Max
        obj_factor *= -1.0
    end
    # Did the user specify a Hessian
    mode = (col_ptr == C_NULL) ? (:Structure) : (:Values)
    if(mode == :Structure)
    	values = unsafe_wrap(Array,values_ptr,0)
		colptr = unsafe_wrap(Array,col_ptr,0)
		rowidx = unsafe_wrap(Array,row_ptr,0)

        @timing prob.prof tic()
		nz = prob.model.str_eval_h(rowid,colid,flag, x0,x1,obj_factor,lambda,mode,rowidx,colptr,values)
        @timing prob.prof prob.t_jl_str_eval_h += toq()

		unsafe_store!(nz_ptr,convert(Cint,nz)::Cint)
		# @show "structure - ", nz
    else
    	nz = unsafe_load(nz_ptr)
    	values = unsafe_wrap(Array, values_ptr, nz)
    	rowidx = unsafe_wrap(Array, row_ptr, nz)
    	colptr = unsafe_wrap(Array, col_ptr, ncol+1)
    	# @show "value - ", nz
        @timing prob.prof tic()
    	prob.model.str_eval_h(rowid,colid,flag,x0,x1,obj_factor,lambda,mode,rowidx,colptr,values)
        @timing prob.prof prob.t_jl_eval_h += toq()
    end
    # Done

    return Int32(1)
end

#write solution
function str_write_solution_wrapper(x_ptr::Ptr{Float64}, y_eq_ptr::Ptr{Float64}, y_ieq_ptr::Ptr{Float64}, cbd::Ptr{CallBackData})
    data = unsafe_load(cbd)
    # @show data
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    rowid = Int(data.row_node_id)
    colid = Int(data.col_node_id)
    @assert rowid == colid

    nx = prob.model.get_num_cols(rowid)
    neq = prob.model.get_num_eq_cons(rowid)
    nieq = prob.model.get_num_ineq_cons(rowid)
    x = unsafe_wrap(Array,x_ptr,nx)
    y_eq = unsafe_wrap(Array,y_eq_ptr,neq)
    y_ieq = unsafe_wrap(Array,y_ieq_ptr,nieq)
    @timing prob.prof tic()
    prob.model.str_write_solution(rowid,x,y_eq,y_ieq)
    @timing prob.prof prob.t_jl_write_solution += toq()

    return Int32(1)
end
###########################################################################
# C function wrappers
###########################################################################

function createProblemStruct(comm::MPI.Comm, model::ModelInterface, prof::Bool)
	# println(" createProblemStruct  -- julia")
	str_init_x0_cb = cfunction(str_init_x0_wrapper, Cint, (Ptr{Float64}, Ptr{CallBackData}) )
    str_prob_info_cb = cfunction(str_prob_info_wrapper, Cint, (Ptr{Cint}, Ptr{Float64}, Ptr{Float64}, Ptr{Cint}, Ptr{Float64}, Ptr{Float64}, Ptr{CallBackData}) )
    str_eval_f_cb = cfunction(str_eval_f_wrapper,Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{CallBackData}) )
    str_eval_g_cb = cfunction(str_eval_g_wrapper,Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{CallBackData}) )
    str_eval_grad_f_cb = cfunction(str_eval_grad_f_wrapper, Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{CallBackData}) )
    str_eval_jac_g_cb = cfunction(str_eval_jac_g_wrapper, Cint, (Ptr{Float64}, Ptr{Float64},
    	Ptr{Cint}, Ptr{Float64}, Ptr{Cint}, Ptr{Cint},
    	Ptr{Cint}, Ptr{Float64}, Ptr{Cint}, Ptr{Cint},
    	Ptr{CallBackData}))
    str_eval_h_cb = cfunction(str_eval_h_wrapper, Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Cint}, Ptr{Float64}, Ptr{Cint}, Ptr{Cint}, Ptr{CallBackData}))
    str_write_solution_cb = cfunction(str_write_solution_wrapper, Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{CallBackData}))

    # println(" callback created ")
    prob = PipsNlpProblemStruct(comm, model, prof)
    # @show prob
    ret = ccall(Libdl.dlsym(libparpipsnlp,:CreatePipsNlpProblemStruct),Ptr{Void},
            (MPI.Comm,
            Cint, Ptr{Void}, Ptr{Void},
	    Ptr{Void}, Ptr{Void}, Ptr{Void},
	    Ptr{Void}, Ptr{Void}, Ptr{Void},Any
            ),
            comm,
            model.get_num_scen(),
            str_init_x0_cb,
            str_prob_info_cb,
            str_eval_f_cb,
            str_eval_g_cb,
            str_eval_grad_f_cb,
            str_eval_jac_g_cb,
            str_eval_h_cb,
            str_write_solution_cb,
            prob
            )
    # println(" ccall CreatePipsNlpProblemStruct done ")
    # @show ret

    if ret == C_NULL
        error("PIPS-NLP: Failed to construct problem.")
    else
        prob.ref = ret
    end
    # @show prob
    # println("end createProblemStruct - julia")
    return prob
end

function solveProblemStruct(prob::PipsNlpProblemStruct)
    # println("solveProblemStruct - julia")
    # @show prob

    ret = ccall(Libdl.dlsym(libparpipsnlp,:PipsNlpSolveStruct), Cint,
            (Ptr{Void},),
            prob.ref)
    # @show ret
    prob.model.set_status(Int(ret))

    prob.t_jl_eval_total = report_total_now(prob)
    # @show prob
    return prob.model.get_status()
end

function freeProblemStruct(prob::PipsNlpProblemStruct)
    # @show "freeProblemStruct"
    ret = ccall(Libdl.dlsym(libparpipsnlp,:FreePipsNlpProblemStruct),
            Void, (Ptr{Void},),
            prob.ref)
    # @show ret
    return ret
end

function report_total_now(prob::PipsNlpProblemStruct)
    total = 0.0
    total += prob.t_jl_init_x0
    total += prob.t_jl_str_prob_info
    total += prob.t_jl_eval_f
    total += prob.t_jl_eval_g
    total += prob.t_jl_eval_grad_f
    total += prob.t_jl_eval_jac_g
    total += prob.t_jl_str_eval_jac_g
    total += prob.t_jl_eval_h
    total += prob.t_jl_str_eval_h
    total += prob.t_jl_write_solution
    return total
end

function t_reset(prob::PipsNlpProblemStruct)
    total = report_total_now(prob)
    prob.t_jl_init_x0  = 0.0
    prob.t_jl_str_prob_info  = 0.0
    prob.t_jl_eval_f = 0.0
    prob.t_jl_eval_g = 0.0
    prob.t_jl_eval_grad_f = 0.0
    prob.t_jl_eval_jac_g = 0.0
    prob.t_jl_str_eval_jac_g = 0.0
    prob.t_jl_eval_h = 0.0
    prob.t_jl_str_eval_h = 0.0
    prob.t_jl_write_solution = 0.0
    return total
end

end


text:
using DiffEqDevTools
using Test

# write your own tests here
@time @testset "Benchmark Tests" begin include("benchmark_tests.jl") end
@time @testset "ODE AppxTrue Tests" begin include("ode_appxtrue_tests.jl") end
@time @testset "Analyticless Convergence Tests" begin include("analyticless_convergence_tests.jl") end
@time @testset "ODE Tableau Convergence Tests" begin include("ode_tableau_convergence_tests.jl") end ## Windows 32-bit fails on Butcher62 convergence test
@time @testset "Analyticless Stochastic WP" begin include("analyticless_stochastic_wp.jl") end
@time @testset "Stability Region Tests" begin include("stability_region_test.jl") end


text:
using HarmonicOrthogonalPolynomials, Plot
plotly()

G = HarmonicOrthogonalPolynomials.grid(SphericalHarmonic()[:,Block.(Base.OneTo(10))])
n,m = size(G)
g = append!(vec(G), [SphericalCoordinate(0,0), SphericalCoordinate(π,0)])
x,y,z = ntuple(k ->getindex.(g,k), 3);
N = length(g)

i = [range(0; length=m, step=n); ]
j = [range(n; length=m-1, step=n); 0]
k = fill(N-2, m)

for ν = range(0; length=m-1, step=n)
    append!(i, range(ν; length=n-1))
    append!(i, range(ν+n; length=n-1))
    append!(j, range(ν+1; length=n-1))
    append!(j, range(ν+1; length=n-1))
    append!(k, range(ν+n; length=n-1))
    append!(k, range(ν+n+1; length=n-1))
end

append!(i, range(n*(m-1); length=n-1))
append!(i, range(0; length=n-1))
append!(j, range(n*(m-1)+1; length=n-1))
append!(j, range(n*(m-1)+1; length=n-1))
append!(k, range(0; length=n-1))
append!(k, range(1; length=n-1))

append!(i, [range(2n-1; length=m-1, step=n); n-1])
append!(j, range(n-1; length=m, step=n))
append!(k, fill(N-1, m))

mesh3d(x,y,z; connections = (i,j,k))


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -2.1106331886844685 <= q <= -0.28632222807890817)

                     add_NL_constraint(m, :(softplus(-0.901869745284984 + -0.09697470954184428*softplus(-0.8276381457841007 + 0.36077566442210207*$(x[1]) + 0.4291228542230603*$(x[2]) + 0.43678775223595334*$(x[3])) + 0.09177414115023685*softplus(0.9224048735895125 + -0.5729901252145462*$(x[1]) + -0.6177157080444822*$(x[2]) + -0.37503852018425166*$(x[3]))) + softplus(-0.5704938352741911 + -0.31710381800885923*softplus(-0.8276381457841007 + 0.36077566442210207*$(x[1]) + 0.4291228542230603*$(x[2]) + 0.43678775223595334*$(x[3])) + -0.16638488372553573*softplus(0.9224048735895125 + -0.5729901252145462*$(x[1]) + -0.6177157080444822*$(x[2]) + -0.37503852018425166*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#
# Model TV for a given pixelManifold if
# * pixelM
# * and given data f
# * weight α
#
rep(d) = (d > 1) ? [ones(Int, d)..., d] : d
res(s::Int) = res([s])
res(s) = (length(s) > 1) ? [s..., length(s)] : s
d = length(size(f))
s = size(f)
M = PowerManifold(pixelM, NestedPowerRepresentation(), s...)
M2 = PowerManifold(pixelM, NestedPowerRepresentation(), res(s)...)
N = TangentBundle(M2)
m2(m) = repeat(m; inner=rep(length(size(m))))
#
# Build TV functionals
#
fidelity(x) = 1 / 2 * distance(M, x, f)^2
function Λ(x)
    return ProductRepr(m2(x), forward_logs(M, x)) # on N=TM, namely in T_xM
end
function prior(x)
    # inner 2-norm over logs, 1-norm over the pixel
    return norm(norm.(Ref(pixelM), x, submanifold_component(N, Λ(x), 2)), 1)
end
cost(x) = (1 / α) * fidelity(x) + prior(x)

proxFidelity(M, m, λ, x) = prox_distance(M, λ / α, f, x, 2)
function proxPriorDual(N, n, λ, ξ)
    return ProductRepr(
        submanifold_component(N, ξ, 1),
        project_collaborative_TV(
            base_manifold(N),
            λ,
            submanifold_component(N, n, 1),
            submanifold_component(N, ξ, 2),
            Inf,
            Inf,
        ),
    )
end
function DΛ(m, ξm)
    return ProductRepr(
        repeat(ξm; inner=rep(length(size(m)))), differential_forward_logs(M, m, ξm)
    )
end
AdjDΛ(m, ξn) = adjoint_differential_forward_logs(M, m, submanifold_component(N, ξn, 2));


text:
#=
  Inner product (dot algorithm)
  Outer product (cross algorithm)
  Generate identity matrix or a given size
  Generate zeroes matrix of a given size
  calculate determinant of a triangular matrix
  calculate determinant via minors
  calculate determinant via Cramers Rule
=#

function identityMatrix(m)
  A = spzeroes(m,m)
  for i = 1:m
    A[i,i] = 1
  end
  return A

function zeroesMatrix(m,n)
  A = spzeroes(m,n)
  return A

function checkCompatability(A,B)
  if lenth(A)[2] == length(B)[1]
    return true
  else
    return false
  end

function dot(a,b)
  #equal size
  try
    for i=1:size(a)
      x = a[i] * b[i]
    end
  finally
    return e
  end
  #print error: dimensions

function cross(a,b)
  #produces a by b matrix
  A = zeroesMatrix(length(a),length(b))
  for i = 1:length(b)
    for j = 1:length(a)
      A[j,i] = b[i] * a[j]
    end
  end
  return A

function multiply(A,B)
  if checkCompatability(A,B) == true
    continue
  else
    return error e
  end
  C = new zeroesMatrix(length(A)[1],length(B)[2])
  #columns of B
  for i = 1:length(B)[2]
    #rows of A
    for j = 1:length(A)[1]
      C[i,j] = dot(A[j,:],B[:,i])
    end
  end
  return C

function determinantTriangular(A)
  a = 1
  i = 1
  while i !> length(A[:,1]) & i !> length(A[1,:])
    a *= A[i,i]
    i + 1
  end
  return a

function determinantCofactors(A)
  a = 0
  #down column one
  for i = 1:length(A)[:,1]
    if length(B)[1,:] != 2
      a += B[] * (-1)^ * determinantCofactors(reshape(B, ))
    else
      return (B[1,1] * B[2,2]) - (B[2,1] * b[1,2])
    end
  end
  return a

max_stars_repo_name: K-P-S/Julia-Scientific-Computing

max_issues_repo_name: K-P-S/Julia-Scientific-Computing

max_forks_repo_name: K-P-S/Julia-Scientific-Computing


text:
function test_true_rates(rate, val, val_eps)
    if iszero(rate)
        val >= rate >= val_eps
    else
        val >= rate > val_eps
    end
end

function test_false_rates(rate, val, val_eps)
    if isone(rate)
        val <= rate <= val_eps
    else
        val <= rate < val_eps
    end
end

n = 1000
targets = rand(0:1, n)
scores = rand(n)

# make non-unique scores
scores[rand(1:1000, 10)] .= 0
scores[rand(1:1000, 10)] .= 1
for k in 1:18
    scores[rand(1:1000, 10)] .= rand()
end

scores_sorted = sort(scores)
rates = collect(0:0.01:1)
rates = vcat(rates, [rates])
ks = 1:10

set_encoding(OneZero())

encs = [
    OneZero(),
    OneMinusOne(),
    OneTwo(),
    OneVsOne(3,4),
    OneVsOne(:three,:four),
    OneVsOne("three","four"),
    OneVsRest(1, [2,3,4]),
    OneVsRest(:one, [:two, :three]),
    OneVsRest("one", ["two", "three"]),
    RestVsOne([1,2,3], 4),
    RestVsOne([:one, :two], :three),
    RestVsOne(["one", "two"], "three")
]

@testset "threshold_at_k: k = $k" for k in ks
    @test threshold_at_k(scores, k) == scores_sorted[end - k + 1]
    @test threshold_at_k(scores, k; rev = true) == scores_sorted[end - k + 1]
    @test threshold_at_k(scores, k; rev = false) == scores_sorted[k]
end

@testset "Thresholds for $enc encoding" for enc in encs
    global targets

    targets = recode.(current_encoding(), enc, targets)
    set_encoding(enc)

    @testset "rate = $(rate)" for rate in rates
        @testset "threshold_at_tpr" begin
            thres = threshold_at_tpr(targets, scores, rate)
            val = true_positive_rate(targets, scores, thres)
            val_eps = true_positive_rate(targets, scores, thres .+ eps.(thres))

            @test all(test_true_rates.(rate, val, val_eps))

            thres = threshold_at_tpr(enc, targets, scores, rate)
            tpr = true_positive_rate(enc, targets, scores, thres)
            tpr_eps = true_positive_rate(enc, targets, scores, thres .+ eps.(thres))

            @test all(test_true_rates.(rate, tpr, tpr_eps))
        end

        @testset "threshold_at_tnr" begin
            thres = threshold_at_tnr(targets, scores, rate)
            val = true_negative_rate(targets, scores, thres)
            val_eps = true_negative_rate(targets, scores, thres .- eps.(thres))

            @test all(test_true_rates.(rate, val, val_eps))

            thres = threshold_at_tnr(enc, targets, scores, rate)
            val = true_negative_rate(enc, targets, scores, thres)
            val_eps = true_negative_rate(enc, targets, scores, thres .- eps.(thres))

            @test all(test_true_rates.(rate, val, val_eps))
        end

        @testset "threshold_at_fpr" begin
            thres = threshold_at_fpr(targets, scores, rate)
            val = false_positive_rate(targets, scores, thres)
            val_eps = false_positive_rate(targets, scores, thres .- eps.(thres))

            @test all(test_false_rates.(rate, val, val_eps))

            thres = threshold_at_fpr(enc, targets, scores, rate)
            val = false_positive_rate(enc, targets, scores, thres)
            val_eps = false_positive_rate(enc, targets, scores, thres .- eps.(thres))

            @test all(test_false_rates.(rate, val, val_eps))
        end

        @testset "threshold_at_fnr" begin
            thres = threshold_at_fnr(targets, scores, rate)
            val = false_negative_rate(targets, scores, thres)
            val_eps = false_negative_rate(targets, scores, thres .+ eps.(thres))

            @test all(test_false_rates.(rate, val, val_eps))

            thres = threshold_at_fnr(enc, targets, scores, rate)
            val = false_negative_rate(enc, targets, scores, thres)
            val_eps = false_negative_rate(enc, targets, scores, thres .+ eps.(thres))

            @test all(test_false_rates.(rate, val, val_eps))
        end
    end
end


text:
using Dates, Gadfly

set_default_plot_size(6inch, 3inch)

# single_datetime, single_boxplot, single_violin

p1 = plot(x=[unix2datetime(100)], y=[10], Geom.point)
p2 = plot(y=randn(100), Geom.boxplot)
p3 = plot(y=randn(100), Geom.violin)

hstack(p1, p2, p3)


text:
#@TODO create nice ticks when ks=:relative
"""
Recipe to plot a DFBand. If ks is set to relative it calculates the relative offset of the k_points to the middle of the k-point path and puts that on the x-axis.

Input:    band::DFBand
Optional: ks=nothing
Kwargs:   fermi=0, -> applies fermi level to band eigenvalues before plotting.
linewidth=2
"""
@recipe function f(band::DFBand, ks=nothing; fermi=0, linewidth=2)
    if ks == :relative_cart
        ks = []
        k_m = band.k_points_cart[div(size(band.k_points_cart)[1] + 1, 2)]
        for k in band.k_points_cart
            push!(ks, norm(k - k_m))
        end
        ks[1:div(length(ks), 2)] = -ks[1:div(length(ks), 2)]
    elseif ks == :relative_cryst
        ks = []
        k_m = band.k_points_cryst[div(size(band.k_points_cryst)[1] + 1, 2)]
        for k in band.k_points_cryst
            push!(ks, norm(k - k_m))
        end
        ks[1:div(length(ks), 2)] = -ks[1:div(length(ks), 2)]
    else
        ks = collect(1:length(band.k_points_cart))
    end
    if fermi != 0
        band = apply_fermi_level(band, fermi)
    end
    linewidth --> linewidth
    title     --> "Eigenvalues"
    yguide    --> "Energy (eV)"
    legend    --> false
    out = band.eigvals
    ks, out
end

"""
To plot multiple bands on one plot.
"""
@recipe function f(bands::Vector{<:DFBand}, ks=nothing)
    for (i, band) in enumerate(bands)
        @series begin
            band, ks
        end
    end
end


text:
using Pkg;Pkg.activate(".")
using Test, Plots, BenchmarkTools, Revise
using VortexDistributions
gr(transpose=true)

## field with 2 points per healing length
Nx = 400; Ny = 400
Lx = 200; Ly = Lx
x = LinRange(-Lx/2,Lx/2, Nx+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, Ny+1)[1:end-1];
psi0 = one.(x*y') |> complex;

## make a field with a dipole
# Periodic boundary conditions
psi = Torus(copy(psi0),x,y)

xp = 99.5
yp = 0
xn = xp
yn = -50
vp = PointVortex(xp,yp,1)
vn = PointVortex(xp,yn,-1)
dip = ScalarVortex([vp;vn])

periodic_dipole!(psi,dip)
heatmap(x,y,angle.(psi.ψ))

## import for debugging
using Parameters, Interpolations
import VortexDistributions:findvortices, zoom_interp, zoom_grid

# @btime vort = findvortices(psi)

## test !periodic detection
vort = findvortices(psi)
@test length(vort) == 0

# @btime findvortices(psi)

## test periodic detection
vort = findvortices(psi,periodic=true)
@test length(vort) == 2

# @btime findvortices(psi,periodic=true)

## single vortex test
psi = Torus(copy(psi0),x,y)

xp = 1.133
yp = 2.787
vp = PointVortex(xp,yp,1)

sp = ScalarVortex(vp)

vortex!(psi,sp)
vfound = findvortices(psi);@show vortex_array(vfound)

## test zoom window requirements and stability for single vortex
# round one:

vort = findvortices_grid(psi)
vort = remove_vortices_edge(vort,psi)

psi_int,xint,yint = zoom_interp(psi.ψ,psi.x,psi.y,vort[1].xv,vort[1].yv,periodic=true)

v1 = findvortices_grid(Torus(psi_int,xint,yint))
vint = remove_vortices_edge(v1,Torus(psi_int,xint,yint))[1]

@show vint
heatmap(xint,yint,angle.(psi_int))
scatter!([vint.xv],[vint.yv],label=false)

## round two
psi_int,xint,yint = zoom_interp(psi_int,xint,yint,vint.xv,vint.yv)
v2 = findvortices_grid(Torus(psi_int,xint,yint))
vint = remove_vortices_edge(v2,Torus(psi_int,xint,yint))[1]
@show vint
heatmap(xint,yint,angle.(psi_int))
scatter!([vint.xv],[vint.yv],label=false)

## test zoom stages on failing cases in real data
using JLD2
loadpath="/Users/abradley/Dropbox/Julia/Vortices - simple results"
f1 = joinpath(loadpath,"nv20-t_0.jld2")
f2 = joinpath(loadpath,"nv20-t_200.jld2")

## load good data
@load f2 𝛹200 x y
psi200 = Torus(𝛹200,x,y)
vfound = findvortices(psi200,periodic=true)
p1=heatmap(x,y,angle.(psi200.ψ))
vdata = vortex_array(vfound)

for j in 1:size(vdata)[1]
    scatter!(p1,[vdata[j,1]],[vdata[j,2]],label=false,color=(vdata[j,3]==1 ? :green : :blue),ms=4,alpha=.6,markerstrokecolor=(vdata[j,3]==1 ? :green : :blue))
end
p1


text:
# Functions for area of various geometric shapes
# Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)

"""
    surfarea_cube(side)

Finds surface area of a cube

# Example

```julia
surfarea_cube(1)  # returns 6
surfarea_cube(3)  # returns 54
surfarea_cube(-1) # returns DomainError
```
"""
function surfarea_cube(side)
    if side < 0
        throw(
            DomainError(
                side,
                "surfarea_cube() only works for non-negative values",
            ),
        )
    end
    return 6(side^2)
end

"""
    surfarea_sphere(side)

Finds surface area of a sphere

# Example

```julia
surfarea_sphere(5)  # returns 314.1592653589793
surfarea_sphere(1)  # returns 12.566370614359172
surfarea_sphere(-1) # returns DomainError
```
"""
function surfarea_sphere(radius)
    if radius < 0
        throw(
            DomainError(
                radius,
                "surfarea_sphere() only works for non-negative values",
            ),
        )
    end
    return 4π * (radius^2)
end

"""
    area_rectangle(length, width)

Finds area of the rectangle

# Example

```julia
area_rectangle(10,20) # returns 200
area_rectangle(-1,-2) # returns DomainError
area_rectangle(1,-2)  # returns DomainError
area_rectangle(-1,2)  # returns DomainError
```
"""
function area_rectangle(length, width)
    if length < 0 || width < 0
        throw(
            DomainError("area_rectangle() only works for non-negative values"),
        )
    end
    return length * width
end

"""
    area_square(side)

Finds area of the area_square

# Example

```julia
area_square(10) # returns 100
area_square(-1) # returns DomainError
```
"""
function area_square(side)
    if side < 0
        throw(DomainError("area_square() only works for non-negative values"))
    end
    return side^2
end

"""
    area_triangle(base, height)

Finds area of the right angled triangle with base height

# Example

```julia
area_triangle(10,10) # returns 50.0
area_triangle(-1,-2) # returns DomainError
area_triangle(1,-2)  # returns DomainError
area_triangle(-1,2)  # returns DomainError
```
"""
function area_triangle(base, height)
    if base < 0 || height < 0
        throw(DomainError("area_triangle() only accepts non-negative values"))
    end
    return (base * height) / 2
end

"""
    area_heron_triangle(side1, side2, side3)

Finds area of a triangle using heron's formula

# Example

```julia
area_heron_triangle(5,12,13) # returns 30.0
area_heron_triangle(-1,-2,1) # returns DomainError
area_heron_triangle(1,-2,1)  # returns DomainError
area_heron_triangle(-1,2,1)  # returns DomainError
```
"""
function area_heron_triangle(side1, side2, side3)
    if side1 < 0 || side2 < 0 || side3 < 0
        throw(
            DomainError(
                "area_triangle_three_sides() only accepts non-negative values",
            ),
        )
    elseif side1 + side2 < side3 ||
           side1 + side3 < side2 ||
           side2 + side3 < side1
        throw(DomainError("Given three sides do not form a triangle"))
    end
    semi_perimeter = (side1 + side2 + side3) / 2
    area = sqrt(
        semi_perimeter *
        (semi_perimeter - side1) *
        (semi_perimeter - side2) *
        (semi_perimeter - side3),
    )
    return area
end

"""
    area_parallelogram(base, height)

Finds area of the parallelogram

# Example

```julia
area_parallelogram(10,20) # returns 200
area_parallelogram(-1,-2) # returns DomainError
area_parallelogram(1,-2)  # returns DomainError
area_parallelogram(-1,2)  # returns DomainError
```
"""
function area_parallelogram(base, height)
    if base < 0 || height < 0
        throw(
            DomainError(
                "area_parallelogram() only accepts non-negative values",
            ),
        )
    end
    return base * height
end

"""
area_trapezium(base1,base2,height)

Finds area of the traπzium

# Example

```julia
area_trapezium(10, 20, 30)  # returns 450.0
area_trapezium(-1, -2, -3)  # returns DomainError
area_trapezium(-1, 2, 3)    # returns DomainError
area_trapezium(1, -2, 3)    # returns DomainError
area_trapezium(1, 2, -3)    # returns DomainError
area_trapezium(-1, -2, 3)   # returns DomainError
area_trapezium(1, -2, -3)   # returns DomainError
area_trapezium(-1, 2, -3)   # returns DomainError
```
"""
function area_trapezium(base1, base2, height)
    if base1 < 0 || base2 < 0 || height < 0
        throw(DomainError("area_trapezium() only accepts non-negative values"))
    end
    return 1 / 2 * (base1 + base2) * height
end

"""
    area_circle(radius)

Finds area of the circle

# Example

```julia
area_circle(20) # returns 1256.6370614359173
area_circle(-1) # returns DomainError
```
"""
function area_circle(radius)
    if radius < 0
        throw(DomainError("area_circle() only accepts non-negative values"))
    end
    return π * radius^2
end

"""
    area_ellipse(radius_x, radius_y)

Finds area of the ellipse

# Example

```julia
area_ellipse(10, 10) # returns 314.1592653589793
area_ellipse(10, 20) # returns 628.3185307179587
area_ellipse(1, -2) # returns DomainError
area_ellipse(-1, 2) # returns DomainError
area_ellipse(-1, -2) # returns DomainError
```
"""
function area_ellipse(radius_x, radius_y)
    if radius_x < 0 || radius_y < 0
        throw(DomainError("area_ellipse() only accepts non-negative values"))
    end
    return π * radius_x * radius_y
end

"""
    area_rhombus(diagonal_1, diagonal_2)

Finds area of the rhombus

# Example

```julia
area_rhombus(10, 20) # returns 100.0
area_rhombus(-1,-2) # returns DomainError
area_rhombus(1,-2)  # returns DomainError
area_rhombus(-1,2)  # returns DomainError
```
"""
function area_rhombus(diagonal_1, diagonal_2)
    if diagonal_1 < 0 || diagonal_2 < 0
        throw(DomainError("area_rhombus() only accepts non-negative values"))
    end
    return 1 / 2 * diagonal_1 * diagonal_2
end

"""
    trapazoidal_area(f, x_start, x_end, steps)

Approximates the area under the curve using the trapezoidal rule
Arguments:
 - f: function for the
 - x_start: starting value for x
 - x_end: ending value for x
 - steps: steps taken while integrating.
"""
function trapazoidal_area(f, x_start, x_end, steps)
    x1 = x_start
    fx1 = f(x_start)
    area = 0.0

    for _ in 1:steps
        x2 = ((x_end - x_start) / steps) + x1
        fx2 = f(x2)
        area += abs(fx2 + fx1) * (x2 - x1) / 2
        x1 = x2
        fx1 = fx2
    end

    return area
end

"""
    area_regular_polygon(sides, side_len)

Finds area of any regular Polygon

# Example

```julia
area_regular_polygon(1, 5)  # returns DomainError
area_regular_polygon(3, 5)  # returns 10.825317547305486
area_regular_polygon(7, 15)  # returns 817.6302999003576
area_regular_polygon(-1, 4) # returns DomainError
area_regular_polygon(4, -3) # returns DomainError
area_regular_polygon(-12, -4) # returns DomainError
```
"""
function area_regular_polygon(sides, side_len)
    sides < 3 && throw(DomainError("A polygon should have atleast 3 sides"))
    (sides < 0 || side_len < 0) && throw(
        DomainError(
            "area_regular_polygon() only works for non-negative values",
        ),
    )
    return sides * side_len^2 * cot(pi / sides) / 4
end

"""
    area_polygon(V)

Finds area of any Polygon given by continuous sequence of vertex coordinates
Arguments:
- coords: x,y co-ordinates of the vertices
    [Vector of Tuples / Matrix with 2 rows or 2 columns]

# Example

```julia
area_polygon([(0, 0), (100, 0), (0, 100)])  # returns 5000.0
area_polygon([0 0;100 0;100 100;0 100])  # returns 10000.0
area_polygon([(6, 4.5), (5, 4.5), (4.5, 5.5), (5, 6.5)]) # returns 1.5
area_polygon([0 0;100 0]) # returns DomainError
area_polygon([(6, 4.63), (5, 4.63)]) # returns DomainError
```
"""
function area_polygon(coords::Matrix{T}) where {T<:Real}
    r, c = size(coords)
    (~((c == 2 && r > 2) || (r == 2 && c > 2))) &&
        throw(DomainError("A polygon should have at least 3 vertices"))

    reshape(coords, (:, 2))
    if ~(coords[1, :] == coords[end, :])
        coords = vcat(coords, coords[1, :]') #to make a loop of the vertices
    end
    return @views abs(sum(coords[1:end-1, :] .* (-coords[2:end, 2:-1:1]))) / 2
end

function area_polygon(coords::Vector{Tuple{T,N}}) where {T,N<:Real}
    length(coords) < 3 &&
        throw(DomainError("A polygon should have at least 3 vertices"))
    if ~(first(coords) == last(coords))
        push!(coords, first(coords)) #to make a loop of the vertices
    end
    res = 0
    for i in 1:length(coords)-1
        res =
            res + coords[i][1] * coords[i+1][2] - coords[i+1][1] * coords[i][2]
    end
    return abs(res) / 2
end

max_stars_repo_name: AugustoCL/Julia

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


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

mutable struct Ipv6ExpressRouteCircuitPeeringConfig <: SwaggerModel
    primaryPeerAddressPrefix::Any # spec type: Union{ Nothing, String } # spec name: primaryPeerAddressPrefix
    secondaryPeerAddressPrefix::Any # spec type: Union{ Nothing, String } # spec name: secondaryPeerAddressPrefix
    microsoftPeeringConfig::Any # spec type: Union{ Nothing, ExpressRouteCircuitPeeringConfig } # spec name: microsoftPeeringConfig
    routeFilter::Any # spec type: Union{ Nothing, RouteFilter } # spec name: routeFilter
    state::Any # spec type: Union{ Nothing, String } # spec name: state

    function Ipv6ExpressRouteCircuitPeeringConfig(;primaryPeerAddressPrefix=nothing, secondaryPeerAddressPrefix=nothing, microsoftPeeringConfig=nothing, routeFilter=nothing, state=nothing)
        o = new()
        validate_property(Ipv6ExpressRouteCircuitPeeringConfig, Symbol("primaryPeerAddressPrefix"), primaryPeerAddressPrefix)
        setfield!(o, Symbol("primaryPeerAddressPrefix"), primaryPeerAddressPrefix)
        validate_property(Ipv6ExpressRouteCircuitPeeringConfig, Symbol("secondaryPeerAddressPrefix"), secondaryPeerAddressPrefix)
        setfield!(o, Symbol("secondaryPeerAddressPrefix"), secondaryPeerAddressPrefix)
        validate_property(Ipv6ExpressRouteCircuitPeeringConfig, Symbol("microsoftPeeringConfig"), microsoftPeeringConfig)
        setfield!(o, Symbol("microsoftPeeringConfig"), microsoftPeeringConfig)
        validate_property(Ipv6ExpressRouteCircuitPeeringConfig, Symbol("routeFilter"), routeFilter)
        setfield!(o, Symbol("routeFilter"), routeFilter)
        validate_property(Ipv6ExpressRouteCircuitPeeringConfig, Symbol("state"), state)
        setfield!(o, Symbol("state"), state)
        o
    end
end # type Ipv6ExpressRouteCircuitPeeringConfig

const _property_map_Ipv6ExpressRouteCircuitPeeringConfig = Dict{Symbol,Symbol}(Symbol("primaryPeerAddressPrefix")=>Symbol("primaryPeerAddressPrefix"), Symbol("secondaryPeerAddressPrefix")=>Symbol("secondaryPeerAddressPrefix"), Symbol("microsoftPeeringConfig")=>Symbol("microsoftPeeringConfig"), Symbol("routeFilter")=>Symbol("routeFilter"), Symbol("state")=>Symbol("state"))
const _property_types_Ipv6ExpressRouteCircuitPeeringConfig = Dict{Symbol,String}(Symbol("primaryPeerAddressPrefix")=>"String", Symbol("secondaryPeerAddressPrefix")=>"String", Symbol("microsoftPeeringConfig")=>"ExpressRouteCircuitPeeringConfig", Symbol("routeFilter")=>"RouteFilter", Symbol("state")=>"String")
Base.propertynames(::Type{ Ipv6ExpressRouteCircuitPeeringConfig }) = collect(keys(_property_map_Ipv6ExpressRouteCircuitPeeringConfig))
Swagger.property_type(::Type{ Ipv6ExpressRouteCircuitPeeringConfig }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_Ipv6ExpressRouteCircuitPeeringConfig[name]))}
Swagger.field_name(::Type{ Ipv6ExpressRouteCircuitPeeringConfig }, property_name::Symbol) =  _property_map_Ipv6ExpressRouteCircuitPeeringConfig[property_name]

const _allowed_Ipv6ExpressRouteCircuitPeeringConfig_state = ["Disabled", "Enabled"]

function check_required(o::Ipv6ExpressRouteCircuitPeeringConfig)
    true
end

function validate_property(::Type{ Ipv6ExpressRouteCircuitPeeringConfig }, name::Symbol, val)
    if name === Symbol("state")
        Swagger.validate_param(name, "Ipv6ExpressRouteCircuitPeeringConfig", :enum, val, _allowed_Ipv6ExpressRouteCircuitPeeringConfig_state)
    end
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _gal_uvw{T<:AbstractFloat}(ra::T, dec::T, pmra::T, pmdec::T,
                                    vrad::T, plx::T, lsr::Bool)
    cosdec = cosd(dec)
    sindec = sind(dec)
    cosra  = cosd(ra)
    sinra  = sind(ra)
    k = 4.740470463533348 # = 149597870.7/(86400*365.25) = 1 AU/year in km/s
    A_G = [[ 0.0548755604, +0.4941094279, -0.8676661490];
           [ 0.8734370902, -0.4448296300, -0.1980763734];
           [ 0.4838350155,  0.7469822445, +0.4559837762]]
    vec2 = k*pmra/plx
    vec3 = k*pmdec/plx
    u = ( A_G[1]*cosra*cosdec + A_G[4]*sinra*cosdec + A_G[7]*sindec)*vrad +
        (-A_G[1]*sinra        + A_G[4]*cosra                       )*vec2 +
        (-A_G[1]*cosra*sindec - A_G[4]*sinra*sindec + A_G[7]*cosdec)*vec3
    v = ( A_G[2]*cosra*cosdec + A_G[5]*sinra*cosdec + A_G[8]*sindec)*vrad +
        (-A_G[2]*sinra        + A_G[5]*cosra                       )*vec2 +
        (-A_G[2]*cosra*sindec - A_G[5]*sinra*sindec + A_G[8]*cosdec)*vec3
    w = ( A_G[3]*cosra*cosdec + A_G[6]*sinra*cosdec + A_G[9]*sindec)*vrad +
        (-A_G[3]*sinra        + A_G[6]*cosra                       )*vec2 +
        (-A_G[3]*cosra*sindec - A_G[6]*sinra*sindec + A_G[9]*cosdec)*vec3
    if lsr
        u += -8.5
        v += 13.38
        w +=  6.49
    end
    return u, v, w
end

"""
    gal_uvw(ra, dec, pmra, pmdec, vrad, plx[, lsr=true]) -> u, v, w

### Purpose ###

Calculate the Galactic space velocity \$(u, v, w)\$ of a star.

### Explanation ###

Calculates the Galactic space velocity \$(u, v, w)\$ of a star given its (1)
coordinates, (2) proper motion, (3) parallax, and (4) radial velocity.

### Arguments ###

User must supply a position, proper motion, radial velocity and parallax.
Either scalars or arrays all of the same length can be supplied.

(1) Position:

* `ra`: right ascension, in degrees
* `dec`: declination, in degrees

(2) Proper Motion

* `pmra`: proper motion in right ascension in arc units (typically
  milli-arcseconds/yr).  If given \$\\mu_\\alpha\$ -- proper motion in seconds of
  time/year -- then this is equal to \$15 \\mu_\\alpha \cos(\\text{dec})\$.
* `pmdec`: proper motion in declination (typically mas/yr).

(3) Radial Velocity

* `vrad`: velocity in km/s

(4) Parallax

* `plx`: parallax with same distance units as proper motion measurements
  typically milliarcseconds (mas)

If you know the distance in parsecs, then set `plx` to \$1000/\\text{distance}\$,
if proper motion measurements are given in milli-arcseconds/yr.

There is an additional optional keyword:

* `lsr` (optional boolean keyword): if this keyword is set to `true`, then the
  output velocities will be corrected for the solar motion
  \$(u, v, w)_\\odot = (-8.5, 13.38, 6.49)\$ (Coşkunoǧlu et al. 2011 MNRAS, 412,
  1237;
  DOI:[10.1111/j.1365-2966.2010.17983.x](http://dx.doi.org/10.1111/j.1365-2966.2010.17983.x))
  to the local standard of rest (LSR).  Note that the value of the solar motion
  through the LSR remains poorly determined.

### Output ###

The 3-tuple \$(u, v, w)\$

* \$u\$: velocity (km/s) positive toward the Galactic *anti*center
* \$v\$: velocity (km/s) positive in the direction of Galactic rotation
* \$w\$: velocity (km/s) positive toward the North Galactic Pole

### Method ###

Follows the general outline of Johnson & Soderblom (1987, AJ, 93, 864;
DOI:[10.1086/114370](http://dx.doi.org/10.1086/114370)) except that \$u\$ is
positive outward toward the Galactic *anti*center, and the J2000 transformation
matrix to Galactic coordinates is taken from the introduction to the Hipparcos
catalog.

### Example ###

Compute the U,V,W coordinates for the halo star HD 6755.  Use values from
Hipparcos catalog, and correct to the LSR.

``` julia
ra=ten(1,9,42.3)*15.; dec = ten(61,32,49.5);
pmra = 627.89;  pmdec = 77.84; # mas/yr
vrad = -321.4; dis = 129; # distance in parsecs
u, v, w = gal_uvw(ra, dec, pmra, pmdec, vrad, 1e3/dis, lsr=true)
# => (118.2110474553902,-466.4828898385057,88.16573278565097)
```

### Notes ###

This function does not take distance as input.  See "Arguments" section above
for how to provide it using parallax argument.

Code of this function is based on IDL Astronomy User's Library.
"""
gal_uvw(ra::Real, dec::Real, pmra::Real, pmdec::Real,
        vrad::Real, plx::Real; lsr::Bool=false) =
            _gal_uvw(promote(float(ra), float(dec), float(pmra), float(pmdec),
                             float(vrad), float(plx))..., lsr)

function gal_uvw{R<:Real,
                 D<:Real,
                 PR<:Real,
                 PD<:Real,
                 VR<:Real,
                 PL<:Real}(ra::AbstractArray{R}, dec::AbstractArray{D},
                           pmra::AbstractArray{PR}, pmdec::AbstractArray{PD},
                           vrad::AbstractArray{VR}, plx::AbstractArray{PL};
                           lsr::Bool=false)
    @assert length(ra) == length(dec) == length(pmra) ==
        length(pmdec) == length(vrad) == length(plx)
    typer = typeof(float(one(R)))
    u = similar(ra, typer)
    v = similar(ra, typer)
    w = similar(ra, typer)
    for i in eachindex(ra)
        u[i], v[i], w[i] = gal_uvw(ra[i], dec[i], pmra[i], pmdec[i],
                                   vrad[i], plx[i], lsr=lsr)
    end
    return u, v, w
end


text:
export  OrthoPoly,
        MultiOrthoPoly,
        Quad,
        OrthoPolyQ,
        Measure,
        MultiMeasure,
        Tensor

struct Measure
    name::String
    w::Function
    dom::Tuple{Float64,Float64}
    symmetric::Bool
    pars::Dict
    function Measure(name::String,w::Function,dom::Tuple{Real,Real},symm::Bool,d::Dict=Dict())
        @assert dom[1] < dom[2] "Invalid support."
        new(lowercase(name),w,Float64.(dom),symm,d)
    end
end

struct MultiMeasure
    name::Vector{String}
    w::Function
    w_uni::Vector{Function}
    dom::Vector{Tuple{Float64,Float64}}
    symmetric::Vector{Bool}
    pars::Vector{Dict}
end

# constructor for classic distributions
function Measure(name::String,d::Dict=Dict())
    name = lowercase(name)
    name == "legendre"         &&  return Measure(name,w_legendre,(-1,1),true,d)
    name == "jacobi"           &&  d[:shape_a] > -1 && d[:shape_b] > -1 && return Measure(name,build_w_jacobi(d[:shape_a], d[:shape_b]),(-1,1),d[:shape_a] == d[:shape_b],d)
    name == "laguerre"         &&  return Measure(name,w_laguerre,(0,Inf),false,d)
    name == "genlaguerre"      &&  return Measure(name,build_w_genlaguerre(d[:shape]),(0,Inf),false,d)
    name == "hermite"          &&  return Measure(name,w_hermite,(-Inf,Inf),true,d)
    name == "genhermite"       &&  return Measure(name,build_w_genhermite(Float64(d[:mu])),(-Inf,Inf),true,d)
    name == "meixnerpollaczek" &&  return Measure(name,build_w_meixner_pollaczek(d[:lambda], d[:phi]),(-Inf,Inf),false,d)
    # measures corresponding to probability density functions:
    name == "gaussian"  &&  return Measure(name,w_gaussian,(-Inf,Inf),true,d)
    name == "uniform01" &&  return Measure(name,w_uniform01,(0,1),true,d)
    name == "beta01"    &&  d[:shape_a] > 0 && d[:shape_b] > 0 &&  return Measure(name,build_w_beta(d[:shape_a],d[:shape_b]),(0,1),d[:shape_a] == d[:shape_b],d)
    name == "gamma"     &&  d[:rate] == 1. && return Measure(name,build_w_gamma(d[:shape]),(0, Inf),false,d)
    name == "logistic"  &&  return Measure(name,w_logistic,(-Inf,Inf),true,d)
    # error handling
    name == "jacobi" && !(d[:shape_a] > -1 && d[:shape_b] > -1) &&  throw(AssertionError("Invalid shape parameters."))
    name == "beta01" && !(d[:shape_a] >  0 && d[:shape_b] >  0) &&  throw(AssertionError("Invalid shape parameters."))
    name == "gamma"  && !(d[:rate] == 1.) &&  error("Rates different from one not supported currently.")
    throw(error("A measure of name $name is not implemented."))
end

struct OrthoPoly
    name::String
    deg::Int64          # maximum degree
    α::Vector{Float64}  # recurrence coefficients
    β::Vector{Float64}  # recurrence coefficients
    meas::Measure
    # inner constructor
    function OrthoPoly(name::String,deg::Int64,α::Vector{Float64},β::Vector{Float64},m::Measure)
        @assert deg >= 0 "Degree has to be non-negative."
        @assert length(α) == length(β) "Different number of recursion coefficients α and β supplied."
        new(lowercase(name),deg,α,β,m)
    end
end

# constructor for classic distributions
function OrthoPoly(name::String,deg::Int64,d::Dict=Dict();Nrec::Int64=deg+1)
  @assert Nrec >= deg + 1 "Not enough recurrence coefficients specified"
  name = lowercase(name)
  name == "legendre"         && return OrthoPoly(name,deg,rm_legendre(Nrec)...,Measure(name,Dict()))
  name == "jacobi"           && return OrthoPoly(name,deg,rm_jacobi(Nrec,Float64(d[:shape_a]), Float64(d[:shape_b]))...,Measure(name,d))
  name == "laguerre"         && return OrthoPoly(name,deg,rm_laguerre(Nrec)...,Measure(name,Dict()))
  name == "genlaguerre"      && return OrthoPoly(name,deg,rm_laguerre(Nrec,Float64(d[:shape]))...,Measure(name,d))
  name == "hermite"          && return OrthoPoly(name,deg,rm_hermite(Nrec)...,Measure(name,Dict()))
  name == "genhermite"       && return OrthoPoly(name,deg,rm_hermite(Nrec,Float64(d[:mu]))...,Measure(name,d))
  name == "meixnerpollaczek" && return OrthoPoly(name,deg,rm_meixner_pollaczek(Nrec,Float64(d[:lambda]),Float64(d[:phi]))...,Measure(name,d))
  # measures corresponding to probability density functions:
  name == "gaussian"  && return OrthoPoly(name,deg,r_scale(1/sqrt(2pi),rm_hermite_prob(Nrec)...)...,Measure(name,Dict()))
  name == "uniform01" && return OrthoPoly(name,deg,r_scale(1.,rm_legendre01(Nrec)...)...,Measure(name,Dict()))
  name == "beta01"    && return OrthoPoly(name,deg,r_scale(1/beta(d[:shape_a],d[:shape_b]),rm_jacobi01(Nrec,Float64(d[:shape_b])-1,Float64(d[:shape_a])-1)...)...,Measure(name,d))
  name == "gamma"     && return OrthoPoly(name,deg,r_scale((d[:rate]^d[:shape])/gamma(d[:shape]),rm_laguerre(Nrec,Float64(d[:shape]-1.))...)...,Measure(name,d))
  name == "logistic"  && return OrthoPoly(name,deg,r_scale(1.,rm_logistic(Nrec)...)...,Measure(name,Dict()))
  # error handling
  error("$name is not yet implemented.")
end

# constructor for known Measure
function OrthoPoly(name::String,deg::Int64,m::Measure;Nrec=deg+1,Nquad=10*Nrec,quadrature::Function=clenshaw_curtis,discretization::Function=stieltjes)
  @assert Nrec >= deg + 1 "Not enough recurrence coefficients specified"
  name = lowercase(name)
  a,b = rm_compute(m,Nrec,Nquad,quadrature=quadrature,discretization=discretization)
  return OrthoPoly(name,deg,a,b,m)
end

# general constructor
function OrthoPoly(name::String,deg::Int64,w::Function,s::Tuple{Real,Real},symm::Bool,d::Dict=Dict();Nrec=deg+1,Nquad=10*Nrec,quadrature::Function=clenshaw_curtis,discretization::Function=stieltjes)
  name = lowercase(name)
  m = Measure(name,w,s,symm,d)
  OrthoPoly(name,deg,m;Nrec=Nrec,Nquad=Nquad,quadrature=quadrature,discretization=discretization)
end

struct Quad
    name::String              # name of quadrature
    Nquad::Int64              # number of qudrature points
    nodes::Vector{Float64}
    weights::Vector{Float64}
    meas::Measure
    function Quad(name::String,N::Int64,nodes::Vector{Float64},weights::Vector{Float64},m::Measure)
        @assert N >= 1 "Number of qudrature points has to be positive"
        @assert length(nodes) == length(weights) "Inconsistent number of nodes and weights inconsistent."
        new(lowercase(name),N,nodes,weights,m)
    end
end

# general constructor
function Quad(N::Int,α::Vector{Float64},β::Vector{Float64},m::Measure)
    @assert length(α) == length(β) "Inconsistent length of recurrence coefficients."
    @assert N <= length(α) - 1 "Requested number of quadrature points $N cannot be provided with $(length(α)) recurrence coefficients"
    n,w = gauss(N,α,β)
    Quad("golubwelsch",N,n,w,m)
end
Quad(N::Int,op::OrthoPoly) = Quad(N,op.α,op.β,op.meas)

function Quad(N::Int64,weight::Function,α::Vector{Float64},β::Vector{Float64},supp::Tuple{Float64,Float64},symm::Bool,d::Dict=Dict())
    m = Measure("fun_"*String(nameof(weight)),weight,supp,symm,d)
    Quad(N,α,β,m)
end

# all-purpose constructor (last resort!)
function Quad(N::Int64,m::Measure;quadrature::Function=clenshaw_curtis)
  n, w = quadgp(m.w,m.dom[1],m.dom[2],N;quadrature=quadrature)
  Quad("quadgp",N,n,w,m)
end

function Quad(N::Int64,weight::Function,supp::Tuple{Real,Real},symm::Bool,d::Dict=Dict();quadrature::Function=clenshaw_curtis)
    @assert N >= 1 "Number of qudrature points has to be positive"
    m = Measure("fun_"*String(nameof(weight)),weight,supp,symm,d)
    Quad(N,m;quadrature=quadrature)
end

# Struct that contains pre-computed nodes and weights
struct OrthoPolyQ
    op::OrthoPoly
    quad::Quad
end

function OrthoPolyQ(op::OrthoPoly,N::Int64)
    q = Quad(N,op.α,op.β,op.meas)
    return OrthoPolyQ(op,q)
end
OrthoPolyQ(op::OrthoPoly) = OrthoPolyQ(op,length(op.α)-1)

function OrthoPolyQ(name::String,N::Int64,d::Dict=Dict();Nrec::Int64=N+1)
    op = OrthoPoly(name,N,d;Nrec=Nrec)
    OrthoPolyQ(op)
end

struct MultiOrthoPoly
name::Vector{String}
deg::Int64
dim::Int64
ind::Matrix{Int64} # multi-index
meas::MultiMeasure
uni::Union{Vector{OrthoPoly},Vector{OrthoPolyQ}}
    function MultiOrthoPoly(uni::Union{Vector{OrthoPoly},Vector{OrthoPolyQ}},deg::Int64)
      t = typeof(uni) == Vector{OrthoPoly}
      degs = [ t ? u.deg : u.op.deg for u in uni ]
      @assert deg <= minimum(degs) "Requested degree $deg is greater than smallest univariate degree $(minimum(degs))."
      Nuni::Int64 = length(uni)

      name_meas = [ t ? u.meas.name : u.op.meas.name for u in uni ]
      supp      = [ t ? u.meas.dom : u.op.meas.dom for u in uni       ]
      symm      = [ t ? u.meas.symmetric : u.op.meas.symmetric for u in uni ]
      pars      = [ t ? u.meas.pars : u.op.meas.pars for u in uni      ]
      w_uni     = [ t ? u.meas.w : u.op.meas.w for u in uni         ]
      w(t) = t ? prod([uni[i].meas.w(t[i]) for i=1:Nuni]) : prod([uni[i].op.meas.w(t[i]) for i=1:Nuni])
      m = MultiMeasure(name_meas,w,w_uni,supp,symm,pars)

      name = [ t ? u.name : u.op.name for u in uni ]
      ind = calculateMultiIndices(Nuni,deg)
      dim = size(ind,1)

      new(name,deg,dim,ind,m,uni)
    end
end

struct Tensor
dim::Int64          # "dimension"
T::SparseVector{Float64,Int64}
get::Function
op::Union{OrthoPolyQ,MultiOrthoPoly}
  function Tensor(dim::Int64,mop::MultiOrthoPoly)
    T = computeTensorizedSP(dim,mop)
    g(a) = getentry(a,T,mop.ind,dim)
    new(dim,T,g,mop)
  end
  function Tensor(dim::Int64,opq::OrthoPolyQ)
    T = computeTensorizedSP(dim,opq)
    g(a) = getentry(a,T,calculateMultiIndices(1,opq.op.deg),dim)
    new(dim,T,g,opq)
  end
end


text:
function Max(K)
	M = -Inf;
	len=length(K)
	for i=1:len
		if K[i]>M
			M=K[i]
		end
	end
	return M
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
using Test, GpABC, DifferentialEquations, Distributions, Distances, LinearAlgebra

@testset "LNA Test" begin

    params = [2.0, 1.0, 15.0, 1.0, 1.0, 1.0, 100.0, 1.0, 1.0, 1.0]
    Tspan = (0.0, 10.0)
    volume = 100.0
    n_samples = 10
    n_design_points = 200
    priors = [Uniform(0., 5.), Uniform(0., 5.), Uniform(10., 20.)]
    distance_metric = euclidean
    saveat = 0.1

    #stochiometry matrix
    S = [1.0 -1.0 0.0 0.0 0.0 0.0;
        0.0 0.0 1.0 -1.0 0.0 0.0;
        0.0 0.0 0.0 0.0 1.0 -1.0]

    # manual input of rates of 3 three example
    reaction_rate_function = function(x,params)
        f = [params[1]/(1+params[7]*x[3]),
            params[4]*x[1],
            params[2]*params[8]*x[1]/(1+params[8]*x[1]),
            params[5]*x[2],
            params[3]*params[9]*x[1]*params[10]*x[2]/(1+params[9]*x[1])/(1+params[10]*x[2]),
            params[6]*x[3]]
    end

    #LNA Mean Var example
    x0 = ([3.0, 2.0, 1.0], 0.4*Matrix(1.0I, 3,3))
    input = LNAInput(params, S, reaction_rate_function, volume)
    lna = compute_LNA(input, x0, Tspan, saveat)
    @test size(lna.traj_means,2) == length(lna.time_points)
    @test length(lna.traj_covars) == length(lna.time_points)
    @test size(lna.traj_means,1) == size(input.S,1)

    lna_trajectories = sample_LNA_trajectories(lna, n_samples)
    @test size(lna_trajectories,1) == size(input.S,1)
    @test size(lna_trajectories,2) == length(lna.time_points)

    #ABC-SMC simulation with LNA
    threshold_schedule = [4.0, 3.0, 2.0]
    n_samples = 100
    n_particles = 500
    n_var_params = length(priors)
    reference_data = sample_LNA_trajectories(lna, n_samples)

    function simulator_function(var_params)
        input = LNAInput(vcat(var_params, params[n_var_params+1:end]), S, reaction_rate_function, volume)
        return get_LNA_trajectories(input, n_samples, x0, Tspan, saveat)
    end

    sim_abcsmc_res = SimulatedABCSMC(reference_data,
    simulator_function,
    priors,
    threshold_schedule,
    n_particles,
    )

    @test size(sim_abcsmc_res.population, 1) > 0

    #ABC-SMC emulation with LNA

    emu_abcsmc_res = EmulatedABCSMC(reference_data,
    simulator_function,
    priors,
    threshold_schedule,
    n_particles,
    n_design_points,
    )

    @test size(emu_abcsmc_res.population, 1) > 0

end


text:
function xxz_longrange(lx::Int, delta::Float64=1.0, r::Float64=0.5)
    (lx < 2) && error("The system size is too small : ", lx)

    sz, sp, sm = sparse.(spinoperators(1/2))

    heis_term = 0.5 * (sp ⊗ sm + sm ⊗ sp) + delta * sz ⊗ sz

    Hmat = heis_term
    # recursive generation of Hamiltonian
    if lx > 2
        for n=3:lx
            Hmat = (Hmat ⊗ I(2)) + (I(2^(n-2)) ⊗ heis_term)
            for l in n-2:-1:1
                d = n-l-1
                Hmat += (r^d) * I(2^(l-1)) ⊗ (
                    delta * (sz ⊗ I(2^d) ⊗ sz) +
                    0.5 * (sp ⊗ I(2^d) ⊗ sm + sm ⊗I(2^d) ⊗ sp)
                )
            end
        end
    end
    Hmat
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
# ---
# title: 322. Coin Change
# id: problem322
# author: Indigo
# date: 2021-06-10
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/coin-change/description/>
# hidden: true
# ---
#
# You are given coins of different denominations and a total amount of money
# _amount_. Write a function to compute the fewest number of coins that you need
# to make up that amount. If that amount of money cannot be made up by any
# combination of the coins, return `-1`.
#
# You may assume that you have an infinite number of each kind of coin.
#
#
#
# **Example 1:**
#
#
#
#     Input: coins = [1,2,5], amount = 11
#     Output: 3
#     Explanation: 11 = 5 + 5 + 1
#
#
# **Example 2:**
#
#
#
#     Input: coins = [2], amount = 3
#     Output: -1
#
#
# **Example 3:**
#
#
#
#     Input: coins = [1], amount = 0
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: coins = [1], amount = 1
#     Output: 1
#
#
# **Example 5:**
#
#
#
#     Input: coins = [1], amount = 2
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= coins.length <= 12`
#   * `1 <= coins[i] <= 231 - 1`
#   * `0 <= amount <= 104`
#
#
## @lc code=start
using LeetCode

function coin_change(coins::Vector{Int}, amount::Int)
    res = OffsetArray(fill(amount + 1, 1 + amount), -1)
    res[0] = 0
    @inbounds for coin in coins, i in 0:(amount - coin)
        res[i + coin] = min(res[i + coin], res[i] + 1)
    end
    return res[end] == amount + 1 ? -1 : res[end]
end
## @lc code=end


text:
@inline function ExaTronKernels.dtrqsol(n::Int,x,
                         p,delta::Float64,
                         tx)
    zero = 0.0
    sigma = zero

    ptx = ddot(n, p, 1, x, 1, tx)
    ptp = ddot(n, p, 1, p, 1, tx)
    xtx = ddot(n, x, 1, x, 1, tx)
    dsq = delta^2

    # Guard against abnormal cases.
    rad = ptx^2 + ptp*(dsq - xtx)
    rad = sqrt(max(rad, zero))

    if ptx > zero
        sigma = (dsq - xtx)/(ptx + rad)
    elseif rad > zero
        sigma = (rad - ptx)/ptp
    else
        sigma = zero
    end
    @synchronize

    return sigma
end


text:
using DistributedArrays

include("kernels_naive.jl") #For validation only
include("kernels_parallel.jl")

################################################################################

#Helper routines that initialize, execute, time, and validate the kernels

function runhplp(n, k, l)
    #Initialize
    hpl(distribute(rand(1, 1)), rand(1), 1, 1) #Precompile
    #A = drandn(n, n, workers(), [nworkers(), 1])
    A = distribute(randn(n, n), procs=workers(), dist=[nworkers(), 1])
    b = randn(n)
    A′= Array(A)

    #Run
    #success = false
    #local t, x
    #while !success
    #    try
    #        success = true
    #    end
    #end
    t = @elapsed try x=hpl(A, b, k, l) end
    #Validate
    #r = b - A′*x
    #ϵ = eps()
    #r₀= norm(r,Inf)
    #nrmA1 = norm(A′, 1)
    #r₁= r₀ / (ϵ*nrmA1*n)
    #r₂= r₀ / (ϵ*nrmA1*norm(x,1))
    #r₃= r₀ / (ϵ*norm(A′,Inf)*norm(x,Inf)*n)
    #err = max(r₁, r₂, r₃)
    # if err ≥ 16
    #     warn("Error $err exceeds allowed value of 16")
    # end
    return t
end

function runrandomupdatep(m, nupdate=4m, seed=1)
    #Initialize
    randomupdate!(dzeros(UInt64,1), 1) #Precompile
    srand(seed)
    T = @DArray [UInt64(i) for i in 1:m]

    @assert T.pids==workers() "Not enough work to distribute across all cores"
    #Run
    t = @elapsed randomupdate!(T, nupdate)

    #Validate
    srand(seed)
    T′ = Array(T)
    randomupdate!(T′, nupdate)

    err = 0
    for i = 1:m
        err += (T′[i] != i)
    end
    err /= m
    if err > 0.01
        warn("Error rate = ", 100err, "%, exceeding 1% allowance")
    end
    return t
end

function runstreamtriadp(m, ntrials=10, validate=true)
    #Initialize
    streamtriad!(distribute([0.0]), distribute([1.0]), 2.0, distribute([3.0]))
    a = dzeros(m)
    b = drandn(m)
    c = drandn(m)
    α = randn()

    @assert a.pids==workers() "Not enough work to parallelize over all cores"
    @assert a.pids==b.pids==c.pids "Cuts are not aligned"

    #Run
    t = Inf
    for i in 1:ntrials
        t0 = @elapsed streamtriad!(a, b, α, c)
        t = min(t, t0)
    end

    #Validate
    if validate
        a′ = Array(a)
        d  = zeros(m)
        streamtriad!(d, Array(b), α, Array(c))

        err = 0.0
        for i in 1:m
            err += abs(d[i]-a′[i])
        end
        if err > 1e-6
            warn("Error = $err exceeds threshold")
        end
    end
    return t
end

function runfftp(m)
    #Initialize
    FFT(distribute([1.0+0.0im]))
    z = @DArray [randn() + im*randn() for i in 1:m]
    z′ = Array(z)

    #Run
    t = @elapsed FFT(z)

    #Validate
    ifft!(Array(z))
    err = maxabs(z′ - z)
    errtol = 16*eps()*log2(m)
    if err > errtol
        warn("Error = $err exceeds tolerance of $errtol")
    end

    return t
end


text:
using LinearAlgebra
using Revise
using Statistics
using Test
using MFGnet

@testset "Float64 3D Array" begin
    A = randn(32,16,2)
    At = MFGnet.linInter1D(0.5,1.,A)
    @test At ≈ mean(A,dims=3)
end

@testset "Float64 2D Array" begin
    B = randn(32,2)
    Bt = MFGnet.linInter1D(0.5,1.,B)
    @test Bt ≈ mean(B,dims=2)
end

@testset "tuple" begin
    A = randn(32,16,2)
    B = randn(32,2)
    AB = MFGnet.linInter1D(0.5,1.,(A,B))
    @test AB[1] ≈ mean(A,dims=3)
    @test AB[2] ≈ mean(B,dims=2)
end

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
using Compose
import Cairo, Fontconfig

set_default_graphic_size(14cm, 10cm)

y = [0.26, 0.5, missing, 0.4, NaN, 0.48, 0.58, 0.83]
p1 = collect(Tuple, zip(1:8, y))
p2 = collect(Tuple, zip(1:9, vcat(NaN, y)))

img = compose(context(units=UnitBox(0, 0, 10, 1)),
    (context(), line([p1]), stroke("black")),
    (context(), line([p2]), stroke("red"))
)

imgs = [SVG("forms_and_nans.svg"), PDF("forms_and_nans.pdf")]

draw.(imgs, [img])


text:
abstract type AbstractBuildToolchain{C} end

struct CMake{C} <: AbstractBuildToolchain{C} end
struct Meson{C} <: AbstractBuildToolchain{C} end

c_compiler(::AbstractBuildToolchain{:clang}) = "clang"
cxx_compiler(::AbstractBuildToolchain{:clang}) = "clang++"
c_compiler(::AbstractBuildToolchain{:gcc}) = "gcc"
cxx_compiler(::AbstractBuildToolchain{:gcc}) = "g++"
fortran_compiler(::AbstractBuildToolchain) = "gfortran"

function cmake_arch(p::AbstractPlatform)
    if arch(p) == "powerpc64le"
        return "ppc64le"
    else
        return arch(p)
    end
end

function cmake_os(p::AbstractPlatform)
    if Sys.islinux(p)
        return "Linux"
    elseif Sys.isfreebsd(p)
        return "FreeBSD"
    elseif Sys.isapple(p)
        return "Darwin"
    elseif Sys.iswindows(p)
        return "Windows"
    end
end

function toolchain_file(bt::CMake, p::AbstractPlatform; is_host::Bool=false)
    target = triplet(p)
    aatarget = aatriplet(p)

    # CMake uses the setting of `HOST_SYSTEM_NAME` and `SYSTEM_NAME` to decide
    # whether the current build is a cross-compilation or not:
    # <https://cmake.org/cmake/help/latest/variable/CMAKE_CROSSCOMPILING.html>.
    # We want to have the host toolchain always setting `HOST_SYSTEM_NAME`, and
    # the target toolchain always setting `SYSTEM_NAME`.
    system_name_var = if is_host
        "CMAKE_HOST_SYSTEM_NAME"
    else
        "CMAKE_SYSTEM_NAME"
    end

    if Sys.isapple(p)
        darwin_ver = something(os_version(p), v"14.5.0")
        maj_ver = darwin_ver.major
        min_ver = darwin_ver.minor
        return """
        # CMake toolchain file for $(c_compiler(bt)) running on $(target)
        set($(system_name_var) $(cmake_os(p)))
        set(CMAKE_SYSTEM_PROCESSOR $(cmake_arch(p)))
        set(CMAKE_SYSTEM_VERSION $(maj_ver).$(min_ver))
        set(DARWIN_MAJOR_VERSION $(maj_ver))
        set(DARWIN_MINOR_VERSION $(min_ver))

        # Enable rpath support
        set(CMAKE_SHARED_LIBRARY_RUNTIME_C_FLAG "-Wl,-rpath,")

        set(CMAKE_SYSROOT /opt/$(aatarget)/$(aatarget)/sys-root)
        set(CMAKE_SYSTEM_FRAMEWORK_PATH
            \${CMAKE_SYSROOT}/System/Library/Frameworks
            \${CMAKE_SYSROOT}/System/Library/PrivateFrameworks
        )
        set(CMAKE_INSTALL_PREFIX \$ENV{prefix})

        set(CMAKE_C_COMPILER   /opt/bin/$(target)/$(aatarget)-$(c_compiler(bt)))
        set(CMAKE_CXX_COMPILER /opt/bin/$(target)/$(aatarget)-$(cxx_compiler(bt)))
        set(CMAKE_Fortran_COMPILER /opt/bin/$(target)/$(aatarget)-$(fortran_compiler(bt)))

        set(CMAKE_LINKER  /opt/bin/$(target)/$(aatarget)-ld)
        set(CMAKE_OBJCOPY /opt/bin/$(target)/$(aatarget)-objcopy)

        set(CMAKE_AR     /opt/bin/$(target)/$(aatarget)-ar)
        set(CMAKE_NM     /opt/bin/$(target)/$(aatarget)-nm)
        set(CMAKE_RANLIB /opt/bin/$(target)/$(aatarget)-ranlib)

        if( \$ENV{CC} MATCHES ccache )
            set_property(GLOBAL PROPERTY RULE_LAUNCH_COMPILE ccache)
        endif()
        """
    else
        return """
        # CMake toolchain file for $(c_compiler(bt)) running on $(target)
        set($(system_name_var) $(cmake_os(p)))
        set(CMAKE_SYSTEM_PROCESSOR $(cmake_arch(p)))

        set(CMAKE_SYSROOT /opt/$(aatarget)/$(aatarget)/sys-root/)
        set(CMAKE_INSTALL_PREFIX \$ENV{prefix})

        set(CMAKE_C_COMPILER   /opt/bin/$(target)/$(aatarget)-$(c_compiler(bt)))
        set(CMAKE_CXX_COMPILER /opt/bin/$(target)/$(aatarget)-$(cxx_compiler(bt)))
        set(CMAKE_Fortran_COMPILER /opt/bin/$(target)/$(aatarget)-$(fortran_compiler(bt)))

        set(CMAKE_LINKER  /opt/bin/$(target)/$(aatarget)-ld)
        set(CMAKE_OBJCOPY /opt/bin/$(target)/$(aatarget)-objcopy)

        set(CMAKE_AR     /opt/bin/$(target)/$(aatarget)-ar)
        set(CMAKE_NM     /opt/bin/$(target)/$(aatarget)-nm)
        set(CMAKE_RANLIB /opt/bin/$(target)/$(aatarget)-ranlib)

        if( \$ENV{CC} MATCHES ccache )
            set_property(GLOBAL PROPERTY RULE_LAUNCH_COMPILE ccache)
        endif()
        """
    end
end

function meson_c_link_args(p::AbstractPlatform)
    if arch(p) == "powerpc64le" && Sys.islinux(p)
        return "'-Wl,-rpath-link,/workspace/destdir/lib64'"
    else
        return ""
    end
end
meson_cxx_link_args(p::AbstractPlatform) = meson_c_link_args(p)
meson_fortran_link_args(p::AbstractPlatform) = meson_c_link_args(p)

# We can run native programs only on
# * i686-linux-gnu
# * x86_64-linux-gnu
# * x86_64-linux-musl
function meson_is_foreign(p::AbstractPlatform)
    if Sys.islinux(p) && proc_family(p) == "intel" && (libc(p) == "glibc" || (libc(p) == "musl" && arch(p) == "x86_64"))
        # Better to explicitly return the string we expect rather than
        # relying on the representation of the boolean values (even though
        # the result is the same)
        return "false"
    else
        return "true"
    end
end

# https://mesonbuild.com/Reference-tables.html#operating-system-names
meson_system(p::AbstractPlatform) = lowercase(cmake_os(p))

# https://github.com/mesonbuild/meson/blob/6e39dcad2fbd8d1c739e262b0e7b7d901cf1ce08/mesonbuild/environment.py#L412-L440
meson_cpu(p::AbstractPlatform) = cmake_arch(p)

# https://mesonbuild.com/Reference-tables.html#cpu-families
function meson_cpu_family(p::AbstractPlatform)
    if arch(p) == "powerpc64le"
        return "ppc64"
    elseif arch(p) == "i686"
        return "x86"
    elseif arch(p) == "x86_64"
        return "x86_64"
    elseif arch(p) == "aarch64"
        return "aarch64"
    elseif startswith(arch(p), "arm")
        return "arm"
    end
end

function toolchain_file(bt::Meson, p::AbstractPlatform)
    target = triplet(p)
    aatarget = aatriplet(p)

    return """
    [binaries]
    c = '/opt/bin/$(target)/$(aatarget)-$(c_compiler(bt))'
    cpp = '/opt/bin/$(target)/$(aatarget)-$(cxx_compiler(bt))'
    fortran = '/opt/bin/$(target)/$(aatarget)-$(fortran_compiler(bt))'
    objc = '/opt/bin/$(target)/$(aatarget)-objc'
    ar = '/opt/bin/$(target)/$(aatarget)-ar'
    ld = '/opt/bin/$(target)/$(aatarget)-ld'
    nm = '/opt/bin/$(target)/$(aatarget)-nm'
    strip = '/opt/bin/$(target)/$(aatarget)-strip'
    pkgconfig = '/usr/bin/pkg-config'

    [properties]
    c_args = []
    cpp_args = []
    fortran_args = []
    c_link_args = [$(meson_c_link_args(p))]
    cpp_link_args = [$(meson_cxx_link_args(p))]
    fortran_link_args = [$(meson_fortran_link_args(p))]
    needs_exe_wrapper = $(meson_is_foreign(p))

    [build_machine]
    system = 'linux'
    cpu_family = 'x86_64'
    cpu = 'x86_64'
    endian = 'little'

    [host_machine]
    system = '$(meson_system(p))'
    cpu_family = '$(meson_cpu_family(p))'
    cpu = '$(meson_cpu(p))'
    endian = 'little'

    [paths]
    prefix = '/workspace/destdir'
    """
end

function generate_toolchain_files!(platform::AbstractPlatform;
                                   toolchains_path::AbstractString,
                                   host_platform::AbstractPlatform = default_host_platform,
                                   )

    # Generate the files fot bot the host and the target platforms
    for p in unique((platform, host_platform))
        dir = joinpath(toolchains_path, triplet(p))
        mkpath(dir)

        for compiler in (:clang, :gcc)
            # Target CMake toolchain
            if platforms_match(p, platform)
                write(joinpath(dir, "target_$(aatriplet(p))_$(compiler).cmake"), toolchain_file(CMake{compiler}(), p; is_host=false))
            end
            # Host CMake toolchain
            if platforms_match(p, host_platform)
                write(joinpath(dir, "host_$(aatriplet(p))_$(compiler).cmake"), toolchain_file(CMake{compiler}(), p; is_host=true))
            end
        end
        write(joinpath(dir, "$(aatriplet(p))_clang.meson"), toolchain_file(Meson{:clang}(), p))
        write(joinpath(dir, "$(aatriplet(p))_gcc.meson"), toolchain_file(Meson{:gcc}(), p))

        symlink_if_exists(target, link) = ispath(joinpath(dir, target)) && symlink(target, link)

        # On FreeBSD and MacOS we actually want to default to clang, otherwise gcc
        if Sys.isbsd(p)
            symlink_if_exists("host_$(aatriplet(p))_clang.cmake", joinpath(dir, "host_$(aatriplet(p)).cmake"))
            symlink_if_exists("target_$(aatriplet(p))_clang.cmake", joinpath(dir, "target_$(aatriplet(p)).cmake"))
            symlink("$(aatriplet(p))_clang.meson", joinpath(dir, "$(aatriplet(p)).meson"))
        else
            symlink_if_exists("host_$(aatriplet(p))_gcc.cmake", joinpath(dir, "host_$(aatriplet(p)).cmake"))
            symlink_if_exists("target_$(aatriplet(p))_gcc.cmake", joinpath(dir, "target_$(aatriplet(p)).cmake"))
            symlink("$(aatriplet(p))_gcc.meson", joinpath(dir, "$(aatriplet(p)).meson"))
        end
    end
end


text:
# Batched matrix multiplication
# 1st argument is produced by NNlib.storage_type(A)
NNlib._batched_gemm!(::Type{<:CuArray}, transA::Char, transB::Char, α::Number, A, B, β::Number, C) =
     CUBLAS.gemm_strided_batched!(transA, transB, α, A, B, β, C)

Base.unsafe_convert(::Type{CuPtr{T}}, A::NNlib.BatchedAdjOrTrans{T}) where {T} =
    Base.unsafe_convert(CuPtr{T}, parent(A))


text:
#################
# DUMMY SOLVERS
#################
import GeoStatsBase: solve, solvesingle

@estimsolver DummyEstimSolver begin end
function solve(problem::EstimationProblem, solver::DummyEstimSolver)
  sdat = data(problem)
  sdom = domain(problem)
  npts = nelms(sdom)
  vars = [v for (v,V) in variables(problem)]
  μvar = Dict(v => fill(mean(sdat[v]),npts) for v in vars)
  σvar = Dict(v => fill(var(sdat[v]),npts) for v in vars)
  EstimationSolution(sdom, μvar, σvar)
end

@simsolver DummySimSolver begin end
function solvesingle(problem::SimulationProblem, covars::NamedTuple,
                     solver::DummySimSolver, preproc)
  reals = map(covars.names) do var
    npts = nelms(domain(problem))
    V    = variables(problem)[var]
    real = vcat(fill(zero(V), npts÷2), fill(one(V), npts÷2))
    var => real
  end
  Dict(reals)
end

###################
# DUMMY ESTIMATOR
###################
import GeoStatsBase: fit, predict, status

struct DummyEstimator end
struct FittedDummyEstimator
  z
end

fit(estimator::DummyEstimator, X, z) = FittedDummyEstimator(z)
function predict(fitted::FittedDummyEstimator, xₒ)
  z  = fitted.z
  μ  = sum(z) / length(z)
  σ² = sum((v - μ)^2 for v in z) / length(z)
  μ, σ²
end
status(fitted::FittedDummyEstimator) = true

#########################
# OTHER EXAMPLE SOLVERS
#########################
@estimsolver ESolver begin
  @param A=1.0
  @param B=2
  @jparam J="foo"
  @global C=true
end

@simsolver SSolver begin
  @param A=1.0
  @param B=2
  @jparam J="foo"
  @global C=true
end

########################
# DUMMY LEARNING MODEL
########################
import MLJModelInterface
const MI = MLJModelInterface

struct DummyModel <: MI.Supervised end
MI.fit(m::DummyModel, v, X, y) = rand(unique(y),length(y)), 0, 0
MI.predict(m::DummyModel, θ, X) = rand(θ, Tables.rowcount(X))
MI.target_scitype(m::DummyModel) = AbstractVector{<:MI.Finite}