text:
@testset "track" begin
    @testset "basic" begin
        @system STrack(Controller) begin
            a => 1 ~ track
            b => 2 ~ track
            c(a, b) => a + b ~ track
        end
        s = instance(STrack)
        @test s.a' == 1 && s.b' == 2 && s.c' == 3
    end

    @testset "cross reference" begin
        @test_throws LoadError @eval @system STrackXRef(Controller) begin
            a(b) => b ~ track
            b(a) => a ~ track
        end
    end

    @testset "self reference without init" begin
        @eval @system STrackSRefWithoutInit(Controller) begin
            a(a) => 2a ~ track
        end
        @test_throws UndefVarError instance(STrackSRefWithoutInit)
    end

    @testset "self reference with init" begin
        @system STrackSRefWithInit(Controller) begin
            a(a) => 2a ~ track(init=1)
        end
        r = simulate(STrackSRefWithInit, stop=3)
        @test r[!, :a] == [2, 4, 8, 16]
    end

    @testset "init" begin
        @system STrackInit(Controller) begin
            a(a) => a + 1 ~ track(init=0)
            b(b) => b + 1 ~ track(init=1)
            c(c) => c + 1 ~ track(init=i)
            d => 1 ~ track(init=i)
            i(t=context.clock.tick) => t + 1 ~ track
        end
        r = simulate(STrackInit, stop=3)
        @test r[!, :a] == [1, 2, 3, 4]
        @test r[!, :b] == [2, 3, 4, 5]
        @test r[!, :c] == [2, 3, 4, 5]
        @test r[!, :d] == [1, 1, 1, 1]
    end

    @testset "minmax" begin
        @system STrackMinMax(Controller) begin
            a => 0 ~ track(min=1)
            b => 0 ~ track(max=2)
            c => 0 ~ track(min=1, max=2)
        end
        s = instance(STrackMinMax)
        @test s.a' == 1
        @test s.b' == 0
        @test s.c' == 1
    end

    @testset "minmax unit" begin
        @system STrackMinMaxUnit(Controller) begin
            a => 0 ~ track(u"m", min=1)
            b => 0 ~ track(u"m", max=2u"cm")
            c => 0 ~ track(u"m", min=1u"cm", max=2)
        end
        s = instance(STrackMinMaxUnit)
        @test s.a' == 1u"m"
        @test s.b' == 0u"m"
        @test s.c' == 1u"cm"
    end

    @testset "round" begin
        @system STrackRound(Controller) begin
            a => 1.5 ~ track(round)
            b => 1.5 ~ track(round=:round)
            c => 1.5 ~ track(round=:ceil)
            d => -1.5 ~ track(round=:floor)
            e => -1.5 ~ track(round=:trunc)
        end
        s = instance(STrackRound)
        @test s.a' === s.b'
        @test s.b' === 2.0
        @test s.c' === 2.0
        @test s.d' === -2.0
        @test s.e' === -1.0
    end

    @testset "round int" begin
        @system STrackRoundInt(Controller) begin
            a => 1.5 ~ track::int(round)
            b => 1.5 ~ track::int(round=:round)
            c => 1.5 ~ track::int(round=:ceil)
            d => -1.5 ~ track::int(round=:floor)
            e => -1.5 ~ track::int(round=:trunc)
        end
        s = instance(STrackRoundInt)
        @test s.a' === s.b'
        @test s.b' === 2
        @test s.c' === 2
        @test s.d' === -2
        @test s.e' === -1
    end

    @testset "round int unit" begin
        @system STrackRoundIntUnit(Controller) begin
            a => 1.5 ~ track::int(u"d", round)
            b => 1.5 ~ track::int(u"d", round=:round)
            c => 1.5 ~ track::int(u"d", round=:ceil)
            d => -1.5 ~ track::int(u"d", round=:floor)
            e => -1.5 ~ track::int(u"d", round=:trunc)
        end
        s = instance(STrackRoundIntUnit)
        @test s.a' === s.b'
        @test s.b' === 2u"d"
        @test s.c' === 2u"d"
        @test s.d' === -2u"d"
        @test s.e' === -1u"d"
    end

    @testset "when" begin
        @system STrackWhen(Controller) begin
            T => true ~ preserve::Bool
            F => false ~ preserve::Bool
            a => 1 ~ track
            b => 1 ~ track(when=true)
            c => 1 ~ track(when=false)
            d => 1 ~ track(when=T)
            e => 1 ~ track(when=F)
            f => 1 ~ track(when=!F)
        end
        s = instance(STrackWhen)
        @test s.a' == 1
        @test s.b' == 1
        @test s.c' == 0
        @test s.d' == 1
        @test s.e' == 0
        @test s.f' == 1
    end

    @testset "when unit" begin
        @system STrackWhenUnit(Controller) begin
            T => true ~ preserve::Bool
            F => false ~ preserve::Bool
            a => 1 ~ track(u"d")
            b => 1 ~ track(u"d", when=true)
            c => 1 ~ track(u"d", when=false)
            d => 1 ~ track(u"d", when=T)
            e => 1 ~ track(u"d", when=F)
            f => 1 ~ track(u"d", when=!F)
        end
        s = instance(STrackWhenUnit)
        @test s.a' == 1u"d"
        @test s.b' == 1u"d"
        @test s.c' == 0u"d"
        @test s.d' == 1u"d"
        @test s.e' == 0u"d"
        @test s.f' == 1u"d"
    end

    @testset "when init" begin
        @system STrackWhenInit(Controller) begin
            T => true ~ preserve::Bool
            F => false ~ preserve::Bool
            I => 1 ~ preserve
            a => 0 ~ track(init=1, when=true)
            b => 0 ~ track(init=1, when=false)
            c => 0 ~ track(init=I, when=T)
            d => 0 ~ track(init=I, when=F)
        end
        s = instance(STrackWhenInit)
        @test s.a' == 0
        @test s.b' == 1
        @test s.c' == 0
        @test s.d' == 1
    end
end


text:
@require Gadfly begin
  displaysize(::Gadfly.Plot) = (535, 379)
  displaytitle(::Gadfly.Plot) = "Gadfly"
end

ext: jl


text:
export Unit, Choice, 𝕀, ℝ
export PrimitiveParam, Param

# # Primitives
# These are primitive parameters.

abstract type PrimitiveParam end

"Nondeterministic choice of true or false"
struct BinaryChoice{T <: Integer} <: PrimitiveParam end
Base.eltype(::Type{BinaryChoice{T}}) where T = T

"A Real-valued variable of type `T` in the unit interval: [0, 1]"
struct Unit{T} <: PrimitiveParam end
Base.eltype(::Type{Unit{T}}) where T = T
const 𝕀 = Unit

# # Families of parameters
"Parameter family"
abstract type Param end

"Nondeterministic choice of values in collection"
struct Choice{T} <: Param
  of::T
end
@inline f(d::Choice, id, ω) =
  resolve(StdNormal(), id, ω) * d.σ + d.μ

"The Real numbers"
struct Reals{T} end
ℝ = Reals


text:
using DataFrames
using Gurobi, JuMP
using LatinHypercubeSampling
using Test
using Plots

"""
Set of examples for which to test different examples.
"""

function resample_test(fn_model)
    n_samples = 1000;
    n_dims = length(fn_model.lbs);
    n_iterations = 4;
    plan, _ = LHCoptim(n_samples, n_dims, 1);
    global X = scaleLHC(plan,[(fn_model.lbs[i], fn_model.ubs[i]) for i=1:n_dims]);
    for i=1:n_iterations
        obj_tree, ineq_trees = fit_fn_model(fn_model, X, weights=weights);
        # Creating JuMP model
        m = Model();
        set_optimizer(m, Gurobi.Optimizer)
        n_vars = length(fn_model.lbs);
        vks = [Symbol("x",i) for i=1:n_vars];
        obj_vks = [Symbol("x",i) for i=1:n_vars+1];
        @variable(m, x[1:n_vars]);
        @variable(m, y);
        @objective(m, Min, y);
        add_feas_constraints!(obj_tree, m, [x; y], obj_vks; M = 1e2);
        for tree in ineq_trees
            add_feas_constraints!(tree, m, x, vks; M = 1e2);
        end
        for tree in eq_trees
            add_feas_constraints!(tree, m, x, vks; M = 1e2);
        end
        for i=1:n_vars # Bounding
            @constraint(m, x[i] <= fn_model.ubs[i]);
            @constraint(m, x[i] >= fn_model.lbs[i]);
        end
        status = optimize!(m)
        println("X values: ", getvalue.(x))
        println("Objective cost: ", getvalue.(y))

        global X = resample(n_samples, fn_model.lbs, fn_model.ubs, X, getvalue.(x);
                     mean_shift=0.5, std_shrink = 0.3)
    end
    # Some post processing
    # for i=1:n_iterations
    #     # fig = scatter(X[1+(i-1)*n_samples: i*n_samples, 1], X[1+(i-1)*n_samples: i*n_samples, 2],
    #     fig = scatter(X[1: i*n_samples, 1], X[1: i*n_samples, 2],
    #     xlim=[fn_model.lbs[1],fn_model.ubs[1]], ylim=[fn_model.lbs[2],fn_model.ubs[2]]);
    #     savefig(fig, "figures/"*fn_model.name*"_plot_"*string(i))
    # end
    return true
end

@test test_sagemark_to_GlobalModel()
# fn_model = sagemark_to_GlobalModel(3, lse=true);
# @test example_fit(fn_model, lnr=base_lnr(false))
# @test example_solve(fn_model, M=1e5)

# Importing MINLP cases
# using MINLPLib
# m = fetch_model("minlp2/blend029")
# using Cbc, Juniper
# ipopt = IpoptSolver(print_level=1)
# gurobi = GurobiSolver()
# juniper = JuniperSolver(ipopt, mip_solver=gurobi)
# setsolver(m, juniper)
# solve(m)


text:
module Transform
using LinearAlgebra

scale2(x::Matrix) = repeat(x,inner=(2,2))
scale3(x::Matrix) = repeat(x,inner=(3,3))
scale4(x::Matrix) = repeat(x,inner=(4,4))
rot90(x::Matrix) = rotr90(x)
rot180(x::Matrix) = rotr90(x,2)
rot270(x::Matrix) = rotr90(x,3)
vflip(x::Matrix) = reverse(x,dims=1)
hflip(x::Matrix) = reverse(x,dims=2)
diagflip1(x::Matrix) = transpose(x)
diagflip1(x::Matrix) = reverse(transpose(x))

end


text:
"""
Implementation of local response normalization.
"""
module LRNs
using Flux

export LRN

"""
    LRN{T,I}

Local response normalization layer. Input `i` is processed according to
out(x,y,f,b) = x * [b + α * sum( i(x, y, f-k÷2:f+k÷2, b)^2 )]^(-β)

Todo: β is currently ignored (always set to 0.5)
"""
struct LRN{T,I}
    b::T
    α::T
    β::T
    k::I
end

"""
    (l::LRN)(i)

Applies a local response normalization layer according to:
out(x,y,f,b) = x * [c + α * sum( i(x, y, f-k÷2:f+k÷2, b)^2 )]^(-β)

Todo: β is currently ignored (always set to 0.5)
"""
function (l::LRN)(x)
    ω = similar(x)
    fsize = size(x, 3)
    depth = l.k ÷ 2
    buffer = similar(x, size(x, 1), size(x, 2), 2*depth+1, size(x, 4))
    for i ∈ 1:fsize
        f_min = max(1, i - depth)
        f_max = min(fsize, i + depth)
        buffer = Flux.Tracker.data(x[:, :, f_min:f_max, :])
        ω[:, :, i, :] = sum(buffer.^2, dims=3)
    end
    return x ./ sqrt.(l.b .+ l.α .* ω)
end # function (l::LRN)

Flux.@treelike LRN

# convenience constructor with default arguments
function LRN(; b=1.0, α=1.0, β=0.5, k=5)
    return LRN(b, α, β, k)
end # function LRN
end # module LRNs


text:
mutable struct MovieLens
    dirpath

    df_rating_master
    df_user_master
    df_item_master

    item2iidx
    user2uidx

    df_rating_train
    df_rating_valid
    df_rating_test
    rating_matrix_train
    rating_matrix_valid
    rating_matrix_test

    MovieLens(dirpath = "../dataset/ml-1m") = new(dirpath)
end

function load_data(mldata::MovieLens)
    mldata.df_rating_master = CSV.read(joinpath(mldata.dirpath, "ratings.dat"), delim="::",
                                        header=[:user_id, :item_id, :rating, :timestamp],
                                        types=Dict(:user_id=>String, :item_id=>String, :rating=>Float64))
    mldata.df_user_master = CSV.read(joinpath(mldata.dirpath, "users.dat"), delim="::",
                                        header=[:user_id, :gender, :age, :occupation, :zipcode],
                                        types=Dict(:user_id=>String, :gender=>String, :age=>Int, :occupation=>Int))
    mldata.df_item_master = CSV.read(joinpath(mldata.dirpath, "movies.dat"), delim="::",
                                        header=[:item_id, :title, :genres],
                                        types=Dict(:item_id=>String, :title=>String, :genres=>String))

    mldata.user2uidx = Dict(uid=>i for (i, uid) in enumerate(mldata.df_user_master[:, :user_id]))
    mldata.item2iidx = Dict(iid=>i for (i, iid) in enumerate(mldata.df_item_master[:, :item_id]))

    return nothing
end

function split_data(mldata::MovieLens; testsize=0.2, validsize=0.1, random_state=46)
    trainsize = 1 - testsize - validsize

    fullsize = size(mldata.df_rating_master)[1]
    index_array = randperm(MersenneTwister(random_state), fullsize)

    train_last_index = floor(Int, fullsize*trainsize)
    valid_last_index = floor(Int, fullsize*(trainsize + validsize))

    mldata.df_rating_train = mldata.df_rating_master[index_array[1:train_last_index], :]
    mldata.df_rating_valid = mldata.df_rating_master[index_array[train_last_index+1:valid_last_index], :]
    mldata.df_rating_test = mldata.df_rating_master[index_array[valid_last_index+1:end], :]

    return nothing
end

function transform_to_rating_matrix(df_rating, user2uidx, item2iidx)
    U = Int[]
    I = Int[]
    R = Float64[]

    for i in 1:size(df_rating)[1]
        row = df_rating[i, :]

        uid = row.user_id
        iid = row.item_id
        r = row.rating

        if r != 0
            push!(U, user2uidx[uid])
            push!(I, item2iidx[iid])
            push!(R, r)
        end
    end

    return sparse(U, I, R)
end

function transform_to_rating_matrix(df_rating, user2uidx, item2iidx, threshold_rating::Real)
    U = Int[]
    I = Int[]
    R = Float64[]

    for i in 1:size(df_rating)[1]
        row = df_rating[i, :]

        uid = row.user_id
        iid = row.item_id
        r = row.rating

        if r >= threshold_rating
            push!(U, user2uidx[uid])
            push!(I, item2iidx[iid])
            push!(R, 1)
        end
    end

    return sparse(U, I, R)
end

function transform_to_rating_matrix(mldata::MovieLens)
    mldata.rating_matrix_train = transform_to_rating_matrix(mldata.df_rating_train, mldata.user2uidx, mldata.item2iidx)
    mldata.rating_matrix_valid = transform_to_rating_matrix(mldata.df_rating_valid, mldata.user2uidx, mldata.item2iidx)
    mldata.rating_matrix_test = transform_to_rating_matrix(mldata.df_rating_test, mldata.user2uidx, mldata.item2iidx)

    return
end

function transform_to_rating_matrix(mldata::MovieLens, threshold_rating)
    mldata.rating_matrix_train = transform_to_rating_matrix(mldata.df_rating_train, mldata.user2uidx, mldata.item2iidx, threshold_rating)
    mldata.rating_matrix_valid = transform_to_rating_matrix(mldata.df_rating_valid, mldata.user2uidx, mldata.item2iidx, threshold_rating)
    mldata.rating_matrix_test = transform_to_rating_matrix(mldata.df_rating_test, mldata.user2uidx, mldata.item2iidx, threshold_rating)

    return
end


text:
//Test 1 in Julia
function a()
	x = 1
		print(x)

end


text:
#
# genmult-tests.jl -
#
# Tests for generalized matrix-vector and matrix-matrix products.
#

module TestingLazyAlgebraGenMult

using Test
using Test: print_test_results
using LazyAlgebra
import LazyAlgebra: GenMult
import .GenMult: Reals, Complexes, Floats, BlasFloat

for (pfx, P) in ((:generic, GenMult.Generic()),
                 (:blas,    GenMult.Blas()),
                 (:linear,  GenMult.Linear()),
                 (:basic,   GenMult.Basic()))
    @eval begin
        $(Symbol(pfx, "_lgemv"))(α, trans, A, x) =
            GenMult._lgemv($P, α, trans, A, x)
        $(Symbol(pfx, "_lgemv!"))(α, trans, A, x, β, y) =
            GenMult._lgemv!($P, α, trans, A, x, β, y)
        $(Symbol(pfx, "_lgemm"))(α, transA, A, transB, B, Nc::Integer=2) =
            GenMult._lgemm($P, α, transA, A, transB, B, Int(Nc))
        $(Symbol(pfx, "_lgemm!"))(α, transA, A, transB, B, β, C) =
            GenMult._lgemm!($P, α, transA, A, transB, B, β, C)
    end
end

_transp(t::Char, A::AbstractMatrix) =
    t == 'N' ? A :
    t == 'T' ? transpose(A) :
    t == 'C' ? A' : error("invalid transpose character")

# in Julia < 0.7 randn() does not generate complex numbers.
_randn(::Type{T}, ::Tuple{}) where {T} = _randn(T)
_randn(::Type{T}, dims::Integer...) where {T} = _randn(T, map(Int, dims))
_randn(::Type{T}) where {T<:AbstractFloat} = randn(T)
_randn(::Type{<:Complex{T}}) where {T<:AbstractFloat} =
    Complex(randn(T), randn(T))
function _randn(::Type{T}, dims::NTuple{N,Int}) where {T,N}
    A = Array{T,N}(undef, dims)
    @inbounds for i in eachindex(A)
        A[i] =  _randn(T)
    end
    return A
end

function ref_lgemv(α::Number,
                   trans::Char,
                   A::DenseArray{Ta,Na},
                   x::DenseArray{Tx,Nx}) where {Ta,Na,Tx,Nx}
    @assert 1 ≤ Nx < Na
    Ny = Na - Nx
    dims = size(A)
    if trans == 'N'
        rows, cols = dims[1:Ny], dims[Ny+1:end]
        @assert size(x) == cols
    else
        rows, cols = dims[1:Nx], dims[Nx+1:end]
        @assert size(x) == rows
    end
    M = reshape(A, (prod(rows), prod(cols)))
    v = reshape(x, length(x))
    if trans == 'N'
        return α*reshape(M*v, rows)
    elseif trans == 'T'
        return α*reshape(transpose(M)*v, cols)
    elseif trans == 'C'
        return α*reshape(M'*v, cols)
    else
        error("invalid transpose character")
    end
end

function ref_lgemm(α::Number,
                   transA::Char,
                   A::DenseArray{Ta,Na},
                   transB::Char,
                   B::DenseArray{Tb,Nb},
                   Nc::Integer=2) where {Ta,Na,Tb,Nb}
    Ni, Nj, Nk = GenMult._lgemm_ndims(Na, Nb, Int(Nc))
    @assert Na == Ni + Nk
    @assert Nb == Nk + Nj
    @assert Nc == Ni + Nj
    Adims = size(A)
    Bdims = size(B)
    if transA == 'N'
        I, K = Adims[1:Ni], Adims[Ni+1:Ni+Nk]
    else
        K, I = Adims[1:Nk], Adims[Nk+1:Nk+Ni]
    end
    if transB == 'N'
        @assert Bdims[1:Nk] == K
        J = Bdims[Nk+1:Nk+Nj]
    else
        @assert Bdims[Nj+1:Nj+Nk] == K
        J = Bdims[1:Nj]
    end
    m = prod(I)
    n = prod(J)
    p = prod(K)
    return α*reshape(_transp(transA, reshape(A, transA == 'N' ? (m, p) : (p, m)))*
                     _transp(transB, reshape(B, transB == 'N' ? (p, n) : (n, p))),
                     (I..., J...))
end

cnv(::Type{T}, x::T) where {T<:AbstractFloat} = x
cnv(::Type{T}, x::Real) where {T<:AbstractFloat} = convert(T, x)
cnv(::Type{T}, x::Complex{<:Real}) where {T<:AbstractFloat} =
    convert(Complex{T}, x)

function similar_values(::Type{T},
                        A::AbstractArray,
                        B::AbstractArray;
                        atol::Real=zero(T),
                        rtol::Real=sqrt(eps(T))) where {T<:AbstractFloat}
    if axes(A) != axes(B)
        return false
    end
    local anrm2::T = 0
    local bnrm2::T = 0
    local cnrm2::T = 0
    n = 0
    for i in eachindex(A, B)
        a, b = cnv(T, A[i]), cnv(T, B[i])
        anrm2 += abs2(a)
        bnrm2 += abs2(b)
        cnrm2 += abs2(a - b)
        n += 1
    end
    return sqrt(cnrm2) ≤ n*atol + rtol*sqrt(max(anrm2, bnrm2))
end

function similar_values(A::AbstractArray{Ta},
                        B::AbstractArray{Tb};
                        kwds...) where {Ta,Tb}
    similar_values(float(real(promote_type(Ta,Tb))), A, B; kwds...)
end

"""

`worst_type(T1,T2)` yields the smallest floating-point real type of
real/complex types `T1` and `T2`.

"""
function worst_type(::Type{Complex{T1}},
                    ::Type{T2}) where {T1<:AbstractFloat,
                                       T2<:AbstractFloat}
    worst_type(T1, T2)
end

function worst_type(::Type{T1},
                    ::Type{Complex{T2}}) where {T1<:AbstractFloat,
                                                T2<:AbstractFloat}
    worst_type(T1, T2)
end

function worst_type(::Type{Complex{T1}},
                    ::Type{Complex{T2}}) where {T1<:AbstractFloat,
                                                T2<:AbstractFloat}
    worst_type(T1, T2)
end

function worst_type(::Type{T1},
                    ::Type{T2}) where {T1<:AbstractFloat,
                                       T2<:AbstractFloat}
    sizeof(T1) ≤ sizeof(T2) ? T1 : T2
end

function test_lgemv(reduced::Bool=false)
    # Notes:
    #  - π is good for testing special (non floating-point) values.
    #  - Use prime numbers for dimensions so that they cannot be split.
    if reduced
        MULTIPLIERS = (0, 1, π)
        DIMENSIONS = (((3,  ), (5,  )),
                      ((2, 3), (5, 7)))
        TRANS = ('N', 'T', 'C')
        TYPES = (Float32, Float64, ComplexF64)
    else
        MULTIPLIERS = (0, 1, π)
        DIMENSIONS = (((3,  ), (5,  )),
                      ((2, 3), (5,  )),
                      ((3,  ), (2, 5)),
                      ((2, 3), (5, 7)))
        TRANS = ('N', 'T', 'C')
        TYPES = (Float32, Float64, ComplexF32, ComplexF64)
    end
    @testset "LGEMV" begin
        @testset for dims in DIMENSIONS,
            Ta in TYPES, transA in TRANS,
            Tx in TYPES,
            α in MULTIPLIERS, β in MULTIPLIERS
            Ta <: Real && transA == 'C' && continue
            m, n = dims
            Ty = promote_type(Ta, Tx)
            A = _randn(Ta, (m..., n...))
            x = _randn(Tx, transA == 'N' ? n : m)
            y = _randn(Ty, transA == 'N' ? m : n)
            Tw = worst_type(Ta, Tx)
            ref = ref_lgemv(α, transA, A, x) + β*y
            if β == 0
                @test similar_values(Tw, ref,   generic_lgemv(α, transA, A, x))
                @test similar_values(Tw, ref,    linear_lgemv(α, transA, A, x))
                @test similar_values(Tw, ref,           lgemv(α, transA, A, x))
                if ndims(A) == 2 && ndims(x) == 1
                    @test similar_values(Tw, ref, basic_lgemv(α, transA, A, x))
                end
                if Ta == Tx == Ty && Ta <: BlasFloat
                    α′ = convert(Ta, α)
                    @test similar_values(Tw, ref,  blas_lgemv(α′,transA, A, x))
                end
            end
            @test similar_values(Tw, ref,    generic_lgemv!(α, transA, A, x, β, deepcopy(y)))
            @test similar_values(Tw, ref,     linear_lgemv!(α, transA, A, x, β, deepcopy(y)))
            @test similar_values(Tw, ref,            lgemv!(α, transA, A, x, β, deepcopy(y)))
            if ndims(A) == 2 && ndims(x) == 1
                @test similar_values(Tw, ref,  basic_lgemv!(α, transA, A, x, β, deepcopy(y)))
            end
            if Ta == Tx == Ty && Ta <: BlasFloat
                α′ = convert(Ta, α)
                β′ = convert(Ta, β)
                @test similar_values(Tw, ref,   blas_lgemv!(α′,transA, A, x, β′,deepcopy(y)))
            end
        end
    end
end

function test_lgemm(reduced::Bool=false)
    if reduced
        MULTIPLIERS = (0, 1, π)
        DIMENSIONS = (((5,  ), (3,  ), (2,  )),
                      ((5, 2), (4,  ), (3,  )),
                      ((4, 5), (3, 2), (2, 3)))
        TRANS = ('N', 'T', 'C')
        TYPES = (Float32, Float64, ComplexF64)
    else
        MULTIPLIERS = (0, 1, π)
        DIMENSIONS = (((5,  ), (3,     ), (2,  )),
                      ((5, 2), (3,     ), (2,  )),
                      ((5, 2), (4, 3   ), (3, 2)),
                      ((3, 4), (2, 2, 3), (2, 3)))
        TRANS = ('N', 'T', 'C')
        TYPES = (Float32, Float64, ComplexF32, ComplexF64)
    end
    @testset "LGEMM" begin
        @testset for dims in DIMENSIONS,
            α in MULTIPLIERS, β in MULTIPLIERS,
            Ta in TYPES, transA in TRANS,
            Tb in TYPES, transB in TRANS
            Ta <: Real && transA == 'C' && continue
            Tb <: Real && transB == 'C' && continue
            m, n, p = dims
            Tc = promote_type(Ta, Tb)
            C = _randn(Tc, (m..., n...))
            Nc = ndims(C)
            A = _randn(Ta, transA == 'N' ? (m..., p...) : (p..., m...))
            B = _randn(Tb, transB == 'N' ? (p..., n...) : (n..., p...))
            Tw = worst_type(Ta, Tb)
            ref = ref_lgemm(α, transA, A, transB, B, Nc) + β*C
            #ref = generic_lgemm(α, transA, A, transB, B) + β*C
            if β == 0
                @test similar_values(Tw, ref,   generic_lgemm(α, transA, A, transB, B, Nc))
                @test similar_values(Tw, ref,    linear_lgemm(α, transA, A, transB, B, Nc))
                @test similar_values(Tw, ref,           lgemm(α, transA, A, transB, B, Nc))
                if ndims(A) == ndims(B) == 2
                    @test similar_values(Tw, ref, basic_lgemm(α, transA, A, transB, B, Nc))
                end
                if Ta == Tb == Tc && Ta <: BlasFloat
                    α′ = convert(Ta, α)
                    @test similar_values(Tw, ref,  blas_lgemm(α′,transA, A, transB, B, Nc))
                end
            end
            @test similar_values(Tw, ref,   generic_lgemm!(α, transA, A, transB, B, β, deepcopy(C)))
            @test similar_values(Tw, ref,    linear_lgemm!(α, transA, A, transB, B, β, deepcopy(C)))
            @test similar_values(Tw, ref,           lgemm!(α, transA, A, transB, B, β, deepcopy(C)))
            if ndims(A) == ndims(B) == 2
                @test similar_values(Tw, ref, basic_lgemm!(α, transA, A, transB, B, β, deepcopy(C)))
            end
            if Ta == Tb == Tc && Ta <: BlasFloat
                α′ = convert(Ta, α)
                β′ = convert(Ta, β)
                @test similar_values(Tw, ref,  blas_lgemm!(α′,transA, A, transB, B, β′,deepcopy(C)))
            end
        end
    end
end

function test_all(reduced::Bool=false)
    print_test_results(test_lgemv(reduced))
    print_test_results(test_lgemm(reduced))
    nothing
end

end # module

TestingLazyAlgebraGenMult.test_all()
nothing


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.11.3
#   kernelspec:
#     display_name: Julia 1.7.2
#     language: julia
#     name: julia-1.7
# ---

# ## Outline
#
# 1. Start From Step 1 Result
# 1. Uncertainty Profiles
# 1. Seasonal Climatologies
# 1. MITprof File Creation

include("ArgoToMITprof_step1.jl")

using Plots
do_plot=true

# ## Load Gridded Fields
#
# - grid : ECCO global grid (LLC90; incl. land mask).
# - monthly T/S mean climatology (jan,feb, .. dec)
# - annual T/S variance climatology (non-seasonal)

gridded_fields=GriddedFields.load();

# ## Uncertainty Profiles
#
# The result is expressed as a least-squares weight.
#
# 1. spatial interpolation
# 2. vertical interpolation
# 3. combine instrumental and representation error
#
# **Note: σT,σS may still need to account for bounds at this point.**

# +
#1. spatial interpolation

(f,i,j,w)=GriddedFields.InterpolationFactors(gridded_fields.Γ,prof.lon,prof.lat)
📚=(f=f,i=i,j=j,w=w)

# +
prof_σT=[GriddedFields.Interpolate(gridded_fields.σT[:,k],📚.f,📚.i,📚.j,📚.w)[1] for k=1:50]
prof_σS=[GriddedFields.Interpolate(gridded_fields.σS[:,k],📚.f,📚.i,📚.j,📚.w)[1] for k=1:50]

prof_σT_step1=prof_σT #for visual verification

prof_σT[1:5]

# +
#2. vertical interpolation

z_std=meta["z_std"]
prof_σT=ArgoTools.interp_z(-gridded_fields.Γ.RC,prof_σT,z_std)
prof_σT[1:5]
# -

if do_plot
    plot(prof_σT_step1,gridded_fields.Γ.RC,label="native", legend = :bottomright,marker=:+, xlabel="degree K", ylabel="depth")
    plot!(prof_σT[1:55],-z_std[1:55],label="interpolated",marker=:x)
end

# +
#3. combine instrumental and representation error

prof_std.Tweight.=1 ./(prof_σT.^2 .+ prof_std.T_ERR.^2)
prof_std.Tweight[1:5]

# +
#4. salinity

prof_σS=ArgoTools.interp_z(-gridded_fields.Γ.RC,prof_σS,z_std)
prof_std.Sweight.=1 ./(prof_σS.^2 .+ prof_std.S_ERR.^2)
prof_std.Sweight[1:5]
# -

# ## Seasonal Climatology Profiles

# +
#4. spatio-temporal interpolation

fac,rec=ArgoTools.monthly_climatology_factors(prof.date[1])
# -

prof_T1=[GriddedFields.Interpolate(gridded_fields.T[rec[1]][:,k],📚.f,📚.i,📚.j,📚.w)[1] for k=1:50]
prof_T2=[GriddedFields.Interpolate(gridded_fields.T[rec[2]][:,k],📚.f,📚.i,📚.j,📚.w)[1] for k=1:50];

prof_std.Testim.=ArgoTools.interp_z(-gridded_fields.Γ.RC,fac[1]*prof_T1+fac[2]*prof_T2,z_std)
prof_std.Testim[1:5]

if do_plot
    plot(prof_std.Testim,-z_std,marker=:x,label="climatology", legend = :bottomright, xlabel="degree C", ylabel="depth")
    plot!(prof_std.T[1:55],-z_std[1:55],marker=:o,label="data")
end

# +
#4. spatio-temporal interpolation

prof_S1=[GriddedFields.Interpolate(gridded_fields.S[rec[1]][:,k],📚.f,📚.i,📚.j,📚.w)[1] for k=1:50]
prof_S2=[GriddedFields.Interpolate(gridded_fields.S[rec[2]][:,k],📚.f,📚.i,📚.j,📚.w)[1] for k=1:50];
# -

prof_std.Sestim.=ArgoTools.interp_z(-gridded_fields.Γ.RC,fac[1]*prof_S1+fac[2]*prof_S2,z_std)
prof_std.Sestim[1:5]

ArgoTools.prof_test_set2!(prof_std,meta);

# ## MITprof File Creation

# +
meta["fileOut"]=joinpath(tempdir(),"MITprof_example.nc")
isfile(meta["fileOut"]) ? rm(meta["fileOut"], force=true) : nothing

MITprof.write(meta,[prof],[prof_std])
# -

# ### Show Global Maps (via interpolation)

# +
function interpolate_map(Γ)
    lon=[i for i=-179.5:1.0:179.5, j=-89.5:1.0:89.5]
    lat=[j for i=-179.5:1.0:179.5, j=-89.5:1.0:89.5]
    (f,i,j,w)=GriddedFields.InterpolationFactors(Γ,vec(lon),vec(lat))
    (f=f,i=i,j=j,w=w,lon=vec(lon[:,1]),lat=vec(lat[1,:]))
end

🌐=interpolate_map(gridded_fields.Γ);
# -

ni,nj=length(🌐.lon),length(🌐.lat)
σTm=Array{Float64,3}(undef,(ni,nj,50))
σSm=Array{Float64,3}(undef,(ni,nj,50))
for k=1:50
    tmp=GriddedFields.Interpolate(gridded_fields.σT[:,k],🌐.f,🌐.i,🌐.j,🌐.w)
    σTm[:,:,k]=reshape(tmp,(ni,nj))
    tmp=GriddedFields.Interpolate(gridded_fields.σS[:,k],🌐.f,🌐.i,🌐.j,🌐.w)
    σSm[:,:,k]=reshape(tmp,(ni,nj))
end

if do_plot
    contourf(🌐.lon,🌐.lat,log10.(transpose(σTm[:,:,20]).^2),clims=(-2.5,1.0),title="log10(temperature variance) at 300m")
end

if do_plot
    tmp1=GriddedFields.Interpolate(gridded_fields.T[6][:,20],🌐.f,🌐.i,🌐.j,🌐.w)
    contourf(🌐.lon,🌐.lat,tmp1,clims=(-2.0,20.0),title="temperature in June at 300m")
end


text:
module stochastic_schroedinger

export schroedinger, schroedinger_dynamic

using ...bases, ...states, ...operators
using ...operators_dense, ...operators_sparse
using ...timeevolution
import ...timeevolution: integrate_stoch, recast!
import ...timeevolution.timeevolution_schroedinger: dschroedinger, dschroedinger_dynamic, check_schroedinger

"""
    stochastic.schroedinger(tspan, state0, H, Hs[; fout, ...])

Integrate stochastic Schrödinger equation.

# Arguments
* `tspan`: Vector specifying the points of time for which the output should
        be displayed.
* `psi0`: Initial state as Ket.
* `H`: Deterministic part of the Hamiltonian.
* `Hs`: Stochastic part(s) of the Hamiltonian (either an operator or a vector
        of operators).
* `fout=nothing`: If given, this function `fout(t, state)` is called every time
        an output should be displayed. ATTENTION: The given state is neither
        normalized nor permanent!
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function schroedinger(tspan, psi0::Ket, H::Operator, Hs::Vector;
                fout::Union{Function,Void}=nothing,
                kwargs...)
    tspan_ = convert(Vector{Float64}, tspan)

    n = length(Hs)
    dstate = copy(psi0)
    x0 = psi0.data
    state = copy(psi0)

    check_schroedinger(psi0, H)
    dschroedinger_determ(t::Float64, psi::Ket, dpsi::Ket) = dschroedinger(psi, H, dpsi)
    dschroedinger_stoch(t::Float64, psi::Ket, dpsi::Ket, index::Int) = dschroedinger_stochastic(psi, Hs, dpsi, index)

    integrate_stoch(tspan_, dschroedinger_determ, dschroedinger_stoch, x0, state, dstate, fout, n; kwargs...)
end
schroedinger(tspan, psi0::Ket, H::Operator, Hs::Operator; kwargs...) = schroedinger(tspan, psi0, H, [Hs]; kwargs...)

"""
    stochastic.schroedinger_dynamic(tspan, state0, fdeterm, fstoch[; fout, ...])

Integrate stochastic Schrödinger equation with dynamic Hamiltonian.

# Arguments
* `tspan`: Vector specifying the points of time for which the output should
        be displayed.
* `psi0`: Initial state.
* `fdeterm`: Function `f(t, psi, u) -> H` returning the deterministic
        (time- or state-dependent) part of the Hamiltonian.
* `fstoch`: Function `f(t, psi, u, du) -> Hs` returning a vector that
        contains the stochastic terms of the Hamiltonian.
* `fout=nothing`: If given, this function `fout(t, state)` is called every time
        an output should be displayed. ATTENTION: The given state is neither
        normalized nor permanent!
* `noise_processes=0`: Number of distinct white-noise processes in the equation.
        This number has to be equal to the total number of noise operators
        returned by `fstoch`. If unset, the number is calculated automatically
        from the function output.
        NOTE: Set this number if you want to avoid an initial calculation of
        the function output!
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function schroedinger_dynamic(tspan, psi0::Ket, fdeterm::Function, fstoch::Function;
                fout::Union{Function,Void}=nothing, noise_processes::Int=0,
                kwargs...)
    tspan_ = convert(Vector{Float64}, tspan)

    if noise_processes == 0
        fs_out = fstoch(0.0, psi0)
        n = length(fs_out)
    else
        n = noise_processes
    end

    dstate = copy(psi0)
    x0 = psi0.data
    state = copy(psi0)

    dschroedinger_determ(t::Float64, psi::Ket, dpsi::Ket) = dschroedinger_dynamic(t, psi, fdeterm, dpsi)
    dschroedinger_stoch(t::Float64, psi::Ket, dpsi::Ket, index::Int) =
        dschroedinger_stochastic(t, psi, fstoch, dpsi, index)

    integrate_stoch(tspan, dschroedinger_determ, dschroedinger_stoch, x0, state, dstate, fout, n; kwargs...)
end

function dschroedinger_stochastic(psi::Ket, Hs::Vector{T}, dpsi::Ket, index::Int) where T <: Operator
    check_schroedinger(psi, Hs[index])
    dschroedinger(psi, Hs[index], dpsi)
end

function dschroedinger_stochastic(t::Float64, psi::Ket, f::Function, dpsi::Ket, index::Int)
    ops = f(t, psi)
    check_schroedinger(psi, ops[index])
    dschroedinger(psi, ops[index], dpsi)
end

recast!(psi::StateVector, x::SubArray{Complex128, 1}) = (x .= psi.data)
recast!(x::SubArray{Complex128, 1}, psi::StateVector) = (psi.data = x)

end # module


text:
"""
    foo(x,y)

Creates a 2-element static array from the scalars `x` and `y`.
"""
function foo(x::Number, y::Number)
    SA[x, y]
end


text:
# Suppose the lanternfish live forever and have unlimited food and space. Would
# they take over the entire ocean?
#
# After 256 days in the example above, there would be a total of 26984457539
# lanternfish!
#
# How many lanternfish would there be after 256 days?

# read input - we're using a vector where the index represents the number of
# days until a fish replicates, and the value is how many fish have that
# internal state. Since Julia has 1-based indexing, we need to add one to our
# inputs. We also ned to make sure our array is large enough to handle the fact
# that new fish need extra time before they can replicate: one element
# representing today, and 9 more elements to represent the 9 days a new fish
# needs before replication.
fish = zeros(UInt64, 10)
foreach(split(readline("input.txt"), ',')) do d
  fish[d[1] - '0' + 1] += 1
end

# Rather than rotating the array around, just keep an index representing
# "today". On each day, the fish at the current index replicate and the index
# moves one spot.
for dayidx in 0:255
  # 1-based indexing is kind of a pain in the ass
  idx = dayidx % 10 + 1

  # replicate today's fish
  fish[(idx + 8) % 10 + 1] += fish[idx]

  # move today's fish ahead 7 days
  fish[(idx + 6) % 10 + 1] += fish[idx]

  # clear today's fish
  fish[idx] = 0
end

println("Number of fish: ", sum(fish))


text:
module NeuralNetDiffEq

using Reexport, Statistics
@reexport using DiffEqBase

using Flux, Zygote, DiffEqSensitivity, ForwardDiff, Random
using DiffEqFlux, Adapt
import Tracker, Optim

abstract type NeuralNetDiffEqAlgorithm <: DiffEqBase.AbstractODEAlgorithm end

struct TerminalPDEProblem{G,F,Mu,Sigma,X,T,P} <: DiffEqBase.DEProblem
    g::G
    f::F
    μ::Mu
    σ::Sigma
    X0::X
    tspan::Tuple{T,T}
    p::P
    TerminalPDEProblem(g,f,μ,σ,X0,tspan,p=nothing) = new{typeof(g),typeof(f),
                                                         typeof(μ),typeof(σ),
                                                         typeof(X0),eltype(tspan),
                                                         typeof(p)}(
                                                         g,f,μ,σ,X0,tspan,p)
end

Base.summary(prob::TerminalPDEProblem) = string(nameof(typeof(prob)))

function Base.show(io::IO, A::TerminalPDEProblem)
  println(io,summary(A))
  print(io,"timespan: ")
  show(io,A.tspan)
end

struct KolmogorovPDEProblem{ Mu, Sigma, Phi, X , T , D ,P} <: DiffEqBase.DEProblem
    μ::Mu
    sigma::Sigma
    phi::Phi
    xspan::Tuple{X,X}
    tspan::Tuple{T,T}
    d::D
    p::P
    KolmogorovPDEProblem( μ, sigma, phi , xspan , tspan , d, p=nothing) = new{typeof(μ),typeof(sigma),typeof(phi),eltype(tspan),eltype(xspan),typeof(d),typeof(p)}(μ,sigma,phi,xspan,tspan,d,p)
end

Base.summary(prob::KolmogorovPDEProblem) = string(nameof(typeof(prob)))
function Base.show(io::IO, A::KolmogorovPDEProblem)
  println(io,summary(A))
  print(io,"timespan: ")
  show(io,A.tspan)
  print(io,"xspan: ")
  show(io,A.xspan)
  println(io , "μ")
  show(io , A.μ)
  println(io,"Sigma")
  show(io , A.sigma)
end

include("ode_solve.jl")
include("pde_solve.jl")
include("pde_solve_ns.jl")
include("second_ode_solve.jl")
include("kolmogorov_solve.jl")

export NNODE, NNODE2, TerminalPDEProblem, NNPDEHan, NNPDENS, KolmogorovPDEProblem, NNKolmogorov

end # module


text:
function prior(theta, lb, ub)
    a = 0.0
    if(all((theta .>= lb) .& (theta .<= ub)))
        a = 1.0
    end
    return a
end

max_stars_repo_name: mcreel/SV

max_issues_repo_name: mcreel/SV

max_forks_repo_name: mcreel/SV


text:
import Base: eltype

# definition of the `Status` struct.
# TODO: this bakes in a lot of assumptions about ordering and padding
if !@isdefined(Status)
    let
        struct_fields = Any[]
        empty_vals = Cint[]
        offset = 0
        i = 0

        while offset < MPI_Status_Source_offset
            push!(struct_fields, :($(Symbol(:_pad, i+=1))::Cint))
            push!(empty_vals, zero(Cint))
            offset += 4
        end
        @assert offset == MPI_Status_Source_offset
        push!(struct_fields, :(source::Cint))
        push!(empty_vals, MPI_ANY_SOURCE)
        offset += 4

        while offset < MPI_Status_Tag_offset
            push!(struct_fields, :($(Symbol(:_pad, i+=1))::Cint))
            push!(empty_vals, zero(Cint))
            offset += 4
        end
        @assert offset == MPI_Status_Tag_offset
        push!(struct_fields, :(tag::Cint))
        push!(empty_vals, MPI_ANY_TAG)
        offset += 4

        while offset < MPI_Status_Error_offset
            push!(struct_fields, :($(Symbol(:_pad, i+=1))::Cint))
            push!(empty_vals, zero(Cint))
            offset += 4
        end
        @assert offset == MPI_Status_Error_offset
        push!(struct_fields, :(error::Cint))
        push!(empty_vals, MPI_SUCCESS)
        offset += 4

        while offset < MPI_Status_size
            push!(struct_fields, :($(Symbol(:_pad, i+=1))::Cint))
            push!(empty_vals, zero(Cint))
            offset += 4
        end

        @eval begin
            struct Status
                $(struct_fields...)
            end
            const STATUS_EMPTY = Status($(empty_vals...))
        end
    end

    @assert fieldoffset(Status, findfirst(isequal(:source), fieldnames(Status))) == MPI_Status_Source_offset
    @assert fieldoffset(Status, findfirst(isequal(:tag), fieldnames(Status))) == MPI_Status_Tag_offset
    @assert fieldoffset(Status, findfirst(isequal(:error), fieldnames(Status))) == MPI_Status_Error_offset
    @assert sizeof(Status) == MPI_Status_size
end

Base.cconvert(::Type{Ptr{Status}}, x::SentinelPtr) = x
Base.unsafe_convert(::Type{Ptr{Status}}, x::SentinelPtr) = reinterpret(Ptr{Status}, x)

"""
    MPI.Status

The status of an MPI receive communication. It has 3 accessible fields

- `source`: source of the received message
- `tag`: tag of the received message
- `error`: error code. This is only set if a function returns multiple statuses.

Additionally, the accessor function [`MPI.Get_count`](@ref) can be used to determine the
number of entries received.
"""
Status

"""
    StatusRef()

Construct an empty `Ref{Status}`. This can be passed to MPI Test/Wait operations.
"""
StatusRef() = Ref(STATUS_EMPTY)

Get_source(status::Status) = Int(status.source)
Get_tag(status::Status) = Int(status.tag)
Get_error(status::Status) = Int(status.error)

"""
    MPI.Request

An MPI Request object, representing a non-blocking communication. This also contains a
reference to the buffer used in the communication to ensure it isn't garbage-collected
during communication.

The status of a Request can be checked by the [`Wait`](@ref) and [`Test`](@ref) functions
or their multiple-request variants, which will deallocate the request once it is
determined to be complete. Alternatively, it will be deallocated at finalization, meaning
that it is safe to ignore the request objects if the status of the communication can be
checked by other means.

See also [`Cancel!`](@ref).
"""
@mpi_handle Request MPI_Request buffer

const REQUEST_NULL = _Request(MPI_REQUEST_NULL, nothing)
Request() = Request(REQUEST_NULL.val, nothing)
isnull(req::Request) = req.val == REQUEST_NULL.val

function free(req::Request)
    if !isnull(req) && !MPI.Finalized()
        @mpichk ccall((:MPI_Request_free, libmpi), Cint, (Ptr{MPI_Request},), req)
        req.buffer = nothing
    end
    return nothing
end

"""
    status = Probe(src::Integer, tag::Integer, comm::Comm)

Blocks until there is a message that can be received matching `src`, `tag` and
`comm`. Returns the corresponding [`Status`](@ref) object.

# External links
$(_doc_external("MPI_Probe"))
"""
function Probe(src::Integer, tag::Integer, comm::Comm)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    @mpichk ccall((:MPI_Probe, libmpi), Cint,
          (Cint, Cint, MPI_Comm, Ptr{Status}),
          src, tag, comm, stat_ref)
    stat_ref[]
end

"""
    ismessage, (status|nothing) = Iprobe(src::Integer, tag::Integer, comm::Comm)

Checks if there is a message that can be received matching `src`, `tag` and `comm`. If so,
returns a tuple `true` and a [`Status`](@ref) object, otherwise returns a tuple `false, nothing`.

# External links
$(_doc_external("MPI_Iprobe"))
"""
function Iprobe(src::Integer, tag::Integer, comm::Comm)
    flag = Ref{Cint}()
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    @mpichk ccall((:MPI_Iprobe, libmpi), Cint,
          (Cint, Cint, MPI_Comm, Ptr{Cint}, Ptr{Status}),
          src, tag, comm, flag, stat_ref)
    if flag[] == 0
        return false, nothing
    end
    true, stat_ref[]
end

"""
    MPI.Get_count(status::Status, T)

The number of entries received. `T` should match the argument provided by the receive call that set the status variable.

If the number of entries received exceeds the limits of the count parameter, then it returns `MPI_UNDEFINED`.

# External links
$(_doc_external("MPI_Get_count"))
"""
function Get_count(stat::Status, datatype::Datatype)
    count = Ref{Cint}()
    @mpichk ccall((:MPI_Get_count, libmpi), Cint,
                  (Ptr{Status}, MPI_Datatype, Ptr{Cint}),
                  Ref(stat), datatype, count)
    Int(count[])
end
Get_count(stat::Status, ::Type{T}) where {T} = Get_count(stat, Datatype(T))

"""
    Wait(req::Request)
    status = Wait(req::Request, Status)

Block until the request `req` is complete and deallocated.

The `Status` argument returns the `Status` of the completed request.

# External links
$(_doc_external("MPI_Wait"))
"""
function Wait(req::Request, status::Union{Ref{Status},Nothing}=nothing)
    # int MPI_Wait(MPI_Request *request, MPI_Status *status)
    @mpichk ccall((:MPI_Wait, libmpi), Cint,
                  (Ptr{MPI_Request}, MPIPtr),
                  req, something(status, MPI_STATUS_IGNORE))
    if isnull(req) # only clear the buffer for non-persistent requests
        req.buffer = nothing
    end
    return nothing
end
function Wait(req::Request, ::Type{Status})
    status = Ref(STATUS_EMPTY)
    Wait(req, status)
    return status[]
end

"""
    flag = Test(req::Request)
    flag, status = Test(req::Request, Status)

Check if the request `req` is complete. If so, the request is deallocated and `flag = true` is returned. Otherwise `flag = false`.

The `Status` argument additionally returns the `Status` of the completed request.

# External links
$(_doc_external("MPI_Test"))
"""
function Test(req::Request, status::Union{Ref{Status},Nothing}=nothing)
    flag = Ref{Cint}()
    # int MPI_Test(MPI_Request *request, int *flag, MPI_Status *status)
    @mpichk ccall((:MPI_Test, libmpi), Cint,
                  (Ptr{MPI_Request}, Ptr{Cint}, Ptr{Status}),
                  req, flag, something(status, MPI_STATUS_IGNORE))
    if isnull(req)
        req.buffer = nothing
    end
    return flag[] != 0
end
function Test(req::Request, ::Type{Status})
    status = Ref(STATUS_EMPTY)
    flag = Test(req, status)
    return flag, status[]
end

"""
    RequestSet(requests::Vector{Request})

A wrapper for an array of `Request`s that can be used to reduce intermediate memory
allocations in [`Waitall`](@ref), [`Testall`](@ref), [`Waitany`](@ref), [`Testany`](@ref),
[`Waitsome`](@ref) or [`Testsome`](@ref).

"""
mutable struct RequestSet <: AbstractVector{Request}
    requests::Vector{Request}
    vals::Vector{MPI_Request}
end

function RequestSet(requests::Vector{Request})
    n = length(requests)
    vals = Vector{MPI_Request}(undef, n)
    for i = 1:n
        vals[i] = requests[i].val
    end
    return RequestSet(requests, vals)
end

Base.length(reqs::RequestSet) = length(reqs.requests)
Base.getindex(reqs::RequestSet, i::Integer) = reqs.requests[i]
function Base.setindex!(reqs::RequestSet, req::Request, i::Integer)
    reqs.vals[i] = req.val
    reqs.requests[i] = req
end

function Base.push!(reqs::RequestSet, req::Request)
    push!(reqs.vals, req.val)
    push!(reqs.requests, req)
end

function update!(reqs::RequestSet, i::Integer)
    req = reqs[i]
    req.val = reqs.vals[i]
    if isnull(req)
        req.buffer = nothing
    end
end
function update!(reqs::RequestSet)
    n = length(reqs)
    for i = 1:n
        update!(reqs, i)
    end
end

"""
    Waitall(reqs::AbstractVector{Request}[, statuses::Vector{Status}])
    statuses = Waitall(reqs::AbstractVector{Request}, Status)

Block until all active requests in the array `reqs` are complete.

The optional `statuses` or `Status` argument can be used to obtain the return `Status` of
each request.

# External links
$(_doc_external("MPI_Waitall"))
"""
function Waitall(reqs::RequestSet, statuses::Union{AbstractVector{Status},Nothing}=nothing)
    n = length(reqs)
    n == 0 && return nothing
    @assert isnothing(statuses) || length(statuses) >= n
    # int MPI_Waitall(int count, MPI_Request array_of_requests[],
    #                 MPI_Status array_of_statuses[])
    @mpichk ccall((:MPI_Waitall, libmpi), Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Status}),
                  n, reqs.vals, something(statuses, MPI_STATUSES_IGNORE))
    update!(reqs)
    return nothing
end
function Waitall(reqs::RequestSet, ::Type{Status})
    statuses = Array{Status}(undef, length(reqs))
    Waitall(reqs, statuses)
    return statuses
end
Waitall(reqs::AbstractVector{Request}, args...) = Waitall(RequestSet(reqs), args...)

"""
    flag = Testall(reqs::AbstractVector{Request}[, statuses::Vector{Status}])
    flag, statuses = Testall(reqs::AbstractVector{Request}, Status)

Check if all active requests in the array `reqs` are complete. If so, the requests are
deallocated and `true` is returned. Otherwise no requests are modified, and `false` is
returned.

The optional `statuses` or `Status` argument can be used to obtain the return `Status` of
each request.

# External links
$(_doc_external("MPI_Testall"))
"""
function Testall(reqs::RequestSet, statuses::Union{AbstractVector{Status},Nothing}=nothing)
    n = length(reqs)
    flag = Ref{Cint}()
    @assert isnothing(statuses) || length(statuses) >= n
    # int MPI_Testall(int count, MPI_Request array_of_requests[], int *flag,
    #                 MPI_Status array_of_statuses[])
    @mpichk ccall((:MPI_Testall, libmpi), Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Cint}, MPIPtr),
                  n, reqs.vals, flag, something(statuses, MPI_STATUSES_IGNORE))
    update!(reqs)
    return flag[] != 0
end
function Testall(reqs::RequestSet, ::Type{Status})
    statuses = Array{Status}(undef, length(reqs))
    flag = Testall(reqs, statuses)
    return flag, statuses
end
Testall(reqs::Vector{Request}, args...) = Testall(RequestSet(reqs), args...)

"""
    i = Waitany(reqs::AbstractVector{Request}[, status::Ref{Status}])
    i, status = Waitany(reqs::AbstractVector{Request}, Status)

Blocks until one of the requests in the array `reqs` is complete: if more than one is
complete, one is chosen arbitrarily. The request is deallocated and the (1-based) index
`i` of the completed request is returned.

If there are no active requests, then `i = nothing`.

The optional `status` argument can be used to obtain the return `Status` of the request.

# External links
$(_doc_external("MPI_Waitany"))
"""
function Waitany(reqs::RequestSet, status::Union{Ref{Status}, Nothing}=nothing)
    ref_idx = Ref{Cint}()
    n = length(reqs)
    # int MPI_Waitany(int count, MPI_Request array_of_requests[], int *index,
    #                 MPI_Status *status)
    @mpichk ccall((:MPI_Waitany, libmpi), Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Cint}, MPIPtr),
                  n, reqs.vals, ref_idx, something(status, MPI_STATUS_IGNORE))
    idx = ref_idx[]
    if idx == MPI_UNDEFINED
        return nothing
    end
    i = Int(idx) + 1
    update!(reqs, i)
    return i
end
function Waitany(reqs::RequestSet, ::Type{Status})
    status = Ref(STATUS_EMPTY)
    i = Waitany(reqs, status)
    return i, status[]
end
Waitany(reqs::Vector{Request}, args...) = Waitany(RequestSet(reqs), args...)

"""
    flag, idx = Testany(reqs::AbstractVector{Request}[, status::Ref{Status}])
    flag, idx, status = Testany(reqs::AbstractVector{Request}, Status)

Checks if any one of the requests in the array `reqs` is complete.

If one or more requests are complete, then one is chosen arbitrarily, deallocated. `flag =
true` and its (1-based) index `idx` is returned.

If there are no completed requests, then `flag = false` and `idx = nothing` is returned.

If there are no active requests, `flag = true` and `idx = nothing`.

The optional `status` argument can be used to obtain the return `Status` of the request.

# External links
$(_doc_external("MPI_Testany"))
"""
function Testany(reqs::RequestSet, status::Union{Ref{Status}, Nothing}=nothing)
    ref_idx = Ref{Cint}()
    rflag = Ref{Cint}()
    n = length(reqs)
    # int MPI_Testany(int count, MPI_Request array_of_requests[], int *index,
    #                 int *flag, MPI_Status *status)
    @mpichk ccall((:MPI_Testany, libmpi), Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Cint}, Ptr{Cint}, MPIPtr),
                  n, reqs.vals, ref_idx, rflag, something(status, MPI_STATUS_IGNORE))
    idx = ref_idx[]
    flag = rflag[] != 0

    if idx == MPI_UNDEFINED
        return flag, nothing
    end
    i = Int(idx) + 1
    update!(reqs, i)
    return flag, i
end
function Testany(reqs::RequestSet, ::Type{Status})
    status = Ref(STATUS_EMPTY)
    flag, i = Testany(reqs, status)
    return flag, i, status[]
end
Testany(reqs::Vector{Request}, args...) = Testany(RequestSet(reqs), args...)

"""
    inds = Waitsome(reqs::AbstractVector{Request}[, statuses::Vector{Status}])

Block until at least one of the active requests in the array `reqs` is complete. The
completed requests are deallocated, and an array `inds` of their indices in `reqs` is returned.

If there are no active requests, then `inds = nothing`.

The optional `statuses` argument can be used to obtain the return `Status` of each
completed request.

# External links
$(_doc_external("MPI_Waitsome"))
"""
function Waitsome(reqs::RequestSet, statuses::Union{AbstractVector{Status},Nothing}=nothing)
    ref_nout = Ref{Cint}()
    n = length(reqs)
    idxs = Vector{Cint}(undef, n)
    @assert isnothing(statuses) || length(statuses) >= n

    # int MPI_Waitsome(int incount, MPI_Request array_of_requests[],
    #                  int *outcount, int array_of_indices[],
    #                  MPI_Status array_of_statuses[])
    @mpichk ccall((:MPI_Waitsome, libmpi), Cint,
          (Cint, Ptr{MPI_Request}, Ptr{Cint}, Ptr{Cint}, Ptr{Status}),
          n, reqs.vals, ref_nout, idxs, something(statuses, MPI_STATUSES_IGNORE))
    nout = Int(ref_nout[])
    # This can happen if there were no valid requests
    if nout == MPI_UNDEFINED
        return nothing
    end
    update!(reqs)
    return [Int(idxs[i]) + 1 for i = 1:nout]
end
function Waitsome(reqs::RequestSet, ::Type{Status})
    statuses = Array{Status}(undef, length(reqs))
    inds = Waitsome(reqs, statuses)
    resize!(statuses, isnothing(inds) ? 0 : length(inds))
    return inds, statuses
end
Waitsome(reqs::Vector{Request}, args...) = Waitsome(RequestSet(reqs), args...)

"""
    inds = Testsome(reqs::AbstractVector{Request}[, statuses::Vector{Status}])

Similar to [`Waitsome`](@ref) except that if no operations have completed it will return
an empty array.

If there are no active requests, then the function returns `nothing`.

The optional `statuses` argument can be used to obtain the return `Status` of each
completed request.

# External links
$(_doc_external("MPI_Testsome"))
"""
function Testsome(reqs::RequestSet, statuses::Union{AbstractVector{Status},Nothing}=nothing)
    ref_nout = Ref{Cint}()
    n = length(reqs)
    idxs = Vector{Cint}(undef, n)
    @assert isnothing(statuses) || length(statuses) >= n

    # int MPI_Testsome(int incount, MPI_Request array_of_requests[],
    #                  int *outcount, int array_of_indices[],
    #                  MPI_Status array_of_statuses[])
    @mpichk ccall((:MPI_Testsome, libmpi), Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Cint}, Ptr{Cint}, Ptr{Status}),
                  n, reqs.vals, ref_nout, idxs, something(statuses, MPI_STATUSES_IGNORE))
    nout = Int(ref_nout[])
    # This can happen if there were no valid requests
    if nout == MPI_UNDEFINED
        return nothing
    end
    update!(reqs)
    return [Int(idxs[i]) + 1 for i = 1:nout]
end
function Testsome(reqs::RequestSet, ::Type{Status})
    statuses = Array{Status}(undef, length(reqs))
    inds = Testsome(reqs, statuses)
    resize!(statuses, isnothing(inds) ? 0 : length(inds))
    return inds, statuses
end
Testsome(reqs::Vector{Request}, args...) = Testsome(RequestSet(reqs), args...)

"""
    Cancel!(req::Request)

Marks a pending [`Irecv!`](@ref) operation for cancellation (cancelling a [`Isend`](@ref),
while supported in some implementations, is deprecated as of MPI 3.1). Note that the
request is not deallocated, and can still be queried using the test or wait functions.

# External links
$(_doc_external("MPI_Cancel"))
"""
function Cancel!(req::Request)
    # int MPI_Cancel(MPI_Request *request)
    @mpichk ccall((:MPI_Cancel, libmpi), Cint, (Ptr{MPI_Request},), req)
    nothing
end


text:
"Partitioned Ordinary Differential Equation"
immutable PODE{T,FQ,FP} <: Equation{T}
    d::Int
    n::Int
    f::FQ
    g::FP
    t₀::T
    q₀::Array{T, 2}
    p₀::Array{T, 2}

    function PODE(d, n, f, g, t₀, q₀, p₀)
        @assert d == size(q₀,1) == size(p₀,1)
        @assert n == size(q₀,2) == size(p₀,2)
        @assert T == eltype(q₀) == eltype(p₀)

        if ndims(q₀) == 1
            q₀ = reshape(q₀, d, n)
        end

        if ndims(p₀) == 1
            p₀ = reshape(p₀, d, n)
        end

        new(d, n, f, g, t₀, q₀, p₀)
    end
end

function PODE{T,FQ,FP}(f::FQ, g::FP, t₀::Real, q₀::DenseArray{T}, p₀::DenseArray{T})
    @assert size(q₀) == size(p₀)
    PODE{T,FQ,FP}(size(q₀, 1), size(q₀, 2), f, g, t₀, q₀, p₀)
end

function PODE(f, g, q₀, p₀)
    PODE(f, g, 0, q₀, p₀)
end

Base.hash(ode::PODE, h::UInt) = hash(ode.d, hash(ode.n, hash(ode.f, hash(ode.g, hash(ode.t₀, hash(ode.q₀, hash(ode.p₀, h)))))))
Base.:(==){T1, T2, FQ1, FQ2, FP1, FP2}(ode1::PODE{T1,FQ1,FP1}, ode2::PODE{T2,FQ2,FP2}) = (
                                ode1.d == ode2.d
                             && ode1.n == ode2.n
                             && ode1.f == ode2.f
                             && ode1.g == ode2.g
                             && ode1.t₀ == ode2.t₀
                             && ode1.q₀ == ode2.q₀
                             && ode1.p₀ == ode2.p₀)


text:
function connect!(c, j, i, σ = 1e-6)
    W = sparse(c.I, c.J, c.W, length(c.rowptr) - 1, length(c.colptr) - 1)
    W[i, j] = σ * randn(SNNFloat)
    c.rowptr, c.colptr, c.I, c.J, c.index, c.W = dsparse(W)
    c.tpre, c.tpost, c.Apre, c.Apost = zero(c.W), zero(c.W), zero(c.W), zero(c.W)
    return nothing
end

function dsparse(A)
    At = sparse(A')
    colptr = A.colptr
    rowptr = At.colptr
    I = rowvals(A)
    V = nonzeros(A)
    J = zero(I)
    # FIXME: Breaks when A is empty
    for j in 1:(length(colptr) - 1)
        J[colptr[j]:(colptr[j+1] - 1)] .= j
    end
    index = zeros(size(I))
    coldown = zeros(eltype(index), length(colptr) - 1)
    for i in 1:(length(rowptr) - 1)
        for st in rowptr[i]:(rowptr[i+1] - 1)
            j = At.rowval[st]
            index[st] = colptr[j] + coldown[j]
            coldown[j] += 1
        end
    end
    # Test.@test At.nzval == A.nzval[index]
    rowptr, colptr, I, J, index, V
end

function record!(obj)
    for (key, val) in obj.records
        if isa(key, Tuple)
            sym, ind = key
            push!(val, getindex(getfield(obj, sym),ind))
        else
            push!(val, copy(getfield(obj, key)))
        end
    end
end

function monitor(obj, keys)
    for key in keys
        if isa(key, Tuple)
            sym, ind = key
        else
            sym = key
        end
        typ = typeof(getfield(obj, sym))
        obj.records[key] = Vector{typ}()
    end
end

function monitor(objs::Array, keys)
    for obj in objs
        monitor(obj, keys)
    end
end

function getrecord(p, sym)
    key = sym
    for (k,val) in p.records
        isa(k, Tuple) && k[1] == sym && (key = k)
    end
    p.records[key]
end

function clear_records(obj)
    for (key, val) in obj.records
        empty!(val)
    end
end

@inline function exp32(x::SNNFloat)
    x = ifelse(x < -10f0, -32f0, x)
    x = 1f0 + x / 32f0
    x *= x; x *= x; x *= x; x *= x; x *= x
    return x
end

@inline function exp256(x::SNNFloat)
    x = ifelse(x < -10f0, -256f0, x)
    x = 1.0 + x / 256.0
    x *= x; x *= x; x *= x; x *= x
    x *= x; x *= x; x *= x; x *= x
    return x
end

macro symdict(x...)
    ex = Expr(:block)
    push!(ex.args, :(d = Dict{Symbol,Any}()))
    for p in x
        push!(ex.args, :(d[$(QuoteNode(p))] = $(esc(p))))
    end
    push!(ex.args, :(d))
    return ex
end


text:
# model.jl

"""
    abstract type AbstractModel end

Implements a statistical model. Has so far only SDEModel as subtype.
"""
abstract type AbstractModel end

"""
    struct SDEModel{T<:Union{Float64, <:Function}} <: AbstractModel
        σ::T
        beginvalue::Float64
        endtime::Float64
        Δ::Float64
    end

# Constructors

```julia
    SDEModel(σ::Q, beginvalue::R, endtime::S, Δ::T) where {Q<:Real, R<:Real,
        S<:Real, T<:Real}
    SDEModel(σ::Q, beginvalue::R, endtime::S, Δ::T) where {Q<:Function, R<:Real,
        S<:Real, T<:Real}
```

Implements an SDE model with begin value X_0=beginvalue (at time zero) variance
σ, up to end time endtime. The SDE is discretised with precision Δ.

# Warning
It is assumed that
for every b under consideration, the laws of dX_t=b(X_t)dt+σ(X_t)dW_t are
equivalent.

# Example

```julia
SDEModel(1.0, 0.0, 10.0, 0.01)
# Implements the model dX_t=b(X_t)dt+σ(X_t)dW_t, with X_0=0.0, up to time 10,
# with precision 0.01.
#
SDEModel(identity, 0.0, 100.0, 0.1)
# Implements the model dX_t=b(X_t)+X_tdW_t, with X_0=0.0, up to time 100,
# discretised with precision 0.1.
```
"""
struct SDEModel{T} <: AbstractModel where T<:Union{Float64, <:Function}
    σ::T
    beginvalue::Float64
    endtime::Float64
    Δ::Float64

    function SDEModel(σ::Q, beginvalue::Float64, endtime::Float64, Δ::Float64) where {
            Q<:Union{Float64, <:Function}}
        endtime > 0 || throw(ArgumentError("Endtime should be positive."))
        Δ > 0 || throw(ArgumentError("Δ should be positive."))
        new{typeof(σ)}(σ, beginvalue, endtime, Δ)
    end
end


text:
#=
# 1D advection matching the example in the book: Nodal Discontinuous Galerkin Method.
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("advection1d");

useLog("dgnewadvection1dlog", level=3)

domain(1)
solverType(DG)
functionSpace(order=3)
timeStepper(LSRK4, cfl=0.01)

mesh(LINEMESH, elsperdim=10, bids=2, interval=[0,2])

u = variable("u")
testSymbol("v")

T = 1;
timeInterval(T)
initial(u, "sin(x)")

boundary(u, 1, DIRICHLET, "-sin(2*pi*t)")
boundary(u, 2, NO_BC)

coefficient("a", 2*pi)
coefficient("beta", 0.2)

###################################################
# Uncomment one of the following two groups of code(or both)

### To generate the Julia code and export it to a set of files, use these two lines.
weakForm(u, "Dt(u*v) + a*grad(u)*v + surface(-a*normal()*u*v + a*normal()*ave(u)*v - 0.5*(1-beta)*a*jump(u)*v)")
exportCode("advec1d")
###

### To import previously generated or modified code, use this line
#importCode("advec1d")
###

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

solve(u);

# Uncomment to plot
# using Plots
# pyplot();
# display(plot(Finch.grid_data.allnodes[1,:], u.values[1,:], marker=:circle, reuse=false))

# @finalize()

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
module CCPH
#Import packages
import Dates, Optim, DifferentialEquations

export Constants, EnvironmentStruct, PhotoKineticRates, PhotoPar, PhotoPar!,
TreeSize, TreePar, HydraulicsPar, CCPHStruct, CCPHOutput, CCPHTS, WeatherTS, CCPH_run,
CCPHTraitmodel, CCPHStandGrowth!

#include external files
include("Auxiliary.jl")
include("ModelStructs.jl")
include("Photosynthesis.jl")
include("Hydraulic.jl")
include("Model.jl")

end


text:
using HDF5
using JLD
using DataFrames

fname = joinpath(tempdir(), "mydata.jld")

df = DataFrame({[2:6], pi*[1:5]})
df2 = DataFrame(a = [1:5], b = pi * [1:5])

file = jldopen(fname, "w")
write(file, "df", df)
write(file, "df2", df2)
close(file)

file = jldopen(fname, "r")
x = read(file, "df")
y = read(file, "df2")
close(file)

using Base.Test
@test isequal(df, x)
@test isequal(df2, y)


text:
using DataFrames, Query

df1 = DataFrame(a=[1,2,3], b=[1.,2.,3.])
df2 = DataFrame(c=[2.,4.,2.], d=["John", "Jim","Sally"])

q = @from i in df1 begin
    @from j in df2
    @select {i.a,i.b,j.c,j.d}
    @collect DataFrame
end

println(q)

source_dict = Dict(:a=>[1,2,3], :b=>[4,5])

q = @from i in source_dict begin
	@from j in i.second
	@select {Key=i.first,Value=j}
	@collect DataFrame
end

println(q)


text:
module TestDataFrameInterface

using FactCheck
using MultidimensionalTables
import DataFrames: DataFrame
import RDatasets: dataset

facts("DataFrameInterface tests") do
  @fact DictArray(DataFrame(Any[collect(1:10), repmat(["x","y"], 5), fill(:sym,10)], [:A,:B,:C])) -->
    @darr(A=collect(1:10), B=repmat(["x","y"],5), C=fill(:sym,10))
  @fact DataFrame(@darr(A=collect(1:10), B=repmat(["x","y"],5), C=fill(:sym,10))) -->
    DataFrame(Any[collect(1:10), repmat(["x","y"], 5), fill(:sym,10)], [:A,:B,:C])
  @fact LabeledArray(DataFrame(Any[collect(1:10), repmat(["x","y"], 5), fill(:sym,10)], [:A,:B,:C])) -->
    @larr(A=collect(1:10), B=repmat(["x","y"],5), C=fill(:sym,10))
  @fact DataFrame(@larr(A=collect(1:10), B=repmat(["x","y"],5), C=fill(:sym,10))) -->
    DataFrame(Any[collect(1:10), repmat(["x","y"], 5), fill(:sym,10)], [:A,:B,:C])

  @fact DataFrame(@darr(A=reshape(collect(1:10),2,5), B=reshape(repmat(["x","y"],5),2,5), C=reshape(fill(:sym,10),2,5))) -->
    DataFrame(Any[collect(1:10), repmat(["x","y"], 5), fill(:sym,10)], [:A,:B,:C])
  @fact DataFrame(@larr(A=reshape(collect(1:10),2,5), B=reshape(repmat(["x","y"],5),2,5), C=reshape(fill(:sym,10),2,5))) -->
    DataFrame(Any[collect(1:10), repmat(["x","y"], 5), fill(:sym,10)], [:A,:B,:C])
  @fact DataFrame(@larr(A=reshape(collect(1:10),2,5), B=reshape(repmat(["x","y"],5),2,5), C=reshape(fill(:sym,10),2,5), axis1[k1=[:a,:b],k2=[10,11]])) -->
    DataFrame(Any[repmat([:a,:b],5),repmat(10:11,5),collect(1:10), repmat(["x","y"], 5), fill(:sym,10)], [:k1,:k2,:A,:B,:C])
  @fact DataFrame(@larr(A=reshape(collect(1:10),2,5), B=reshape(repmat(["x","y"],5),2,5), C=reshape(fill(:sym,10),2,5), axis1[k1=[:a,:b],k2=[10,11]], axis2[r=['A',3,"5",:x,1.0]])) -->
    DataFrame(Any[repmat([:a,:b],5),repmat(10:11,5),repeat(['A',3,"5",:x,1.0],inner=[2]),collect(1:10), repmat(["x","y"], 5), fill(:sym,10)], [:k1,:k2,:r,:A,:B,:C])

  #using MultidimensionalTables;using DataFrames;using RDatasets
  context("RDatasets readability tests") do
    #for pd in zip(collect(values(peel(@select(@darr(RDatasets.datasets()), :Package, :Dataset))))...)
    #  d = @darr(dataset(map(x->x.value, pd)...))
    #  l = @larr(dataset(map(x->x.value, pd)...))
    #  nothing
    #end
    @fact (dataset("datasets", "iris");nothing) --> nothing
    @fact @darr(dataset("datasets", "iris")) --> convert(DictArray, dataset("datasets", "iris"))
    @fact @larr(dataset("datasets", "iris")) --> convert(LabeledArray, dataset("datasets", "iris"))

    @fact (dataset("boot", "neuro");nothing) --> nothing
    @fact @darr(dataset("boot", "neuro")) --> convert(DictArray, dataset("boot", "neuro"))
    @fact @larr(dataset("boot", "neuro")) --> convert(LabeledArray, dataset("boot", "neuro"))
  end
end

end


text:
struct Variable{Name} <: AbstractVariable{Name}
end

copy(v::Variable) = v

macro polyvar(names...)
    exprs = [
        quote
            $(esc(name)) = Variable{$(esc(Expr(:quote, name)))}()
        end
        for name in names
    ]
    Expr(:block, exprs...)
end

struct Power{V}
    exponent::Int

    function Power{V}(e::Int) where {V}
        new{typeassert(V, Variable)}(e)
    end
end

variable(p::Power{V}) where {V} = V
exponent(p::Power) = p.exponent

struct Monomial{Vars, Exponents} <: AbstractMonomial{Vars}
    function Monomial{V, E}() where {V, E}
        @assert length(V) == length(E)
        new{typeassert(V, Tuple{Vararg{Variable}}),
            typeassert(E, Tuple{Vararg{Int}})}()
    end
end

exponents(m::Monomial{V, E}) where {V, E} = E

@generated function powers(::Type{Monomial{V, E}}) where {V, E}
    Expr(:tuple, [:($(Power{v}(e))) for (v, e) in zip(V, E)]...)
end

struct Term{T, M <: Monomial} <: AbstractTerm{T, M}
    coefficient::T
end

Term{M}(c::T) where {T, M <: Monomial} = Term{T, M}(c)
Term(c::T, m::M) where {T, M <: Monomial} = Term{T, M}(c)

coefficient(t::Term) = t.coefficient
monomial(::Type{Term{T, M}}) where {T, M} = M()
monomial(t::Term) = monomial(typeof(t))

struct Polynomial{Terms <: Tuple{Vararg{Term}}} <: AbstractPolynomial
    terms::Terms
end

terms(p::Polynomial) = p.terms

function literal_pow(^, v::V, ::Type{Val{x}}) where {V <: Variable, x}
    Monomial{(v,), (x,)}()
end

const MonomialLike = Union{<:Variable, <:Monomial}
const TermLike = Union{<:MonomialLike, <:Term}
const PolynomialLike = Union{<:TermLike, <:Polynomial}


text:
module graphops

using PyCall
include("utils.jl")

export
    # graph_transformer
    optimize,
    simplify,
    fuse_einsums,
    # graph_dedup
    dedup,
    # expr_generator
    rewrite_einsum_expr,
    # graph_inv_optimizer
    optimize_inverse,
    prune_inv_node,
    # optimal_tree
    generate_optimal_tree,
    split_einsum,
    generate_optimal_tree_w_constraint,
    # graph_als_optimizer
    generate_sequential_optimal_tree

const eingraph = PyNULL()
const gops = PyNULL()

function __init__()
    copy!(eingraph, pyimport_conda("autohoot.einsum_graph", "autohoot"))
    copy!(gops, pyimport_conda("autohoot.graph_ops", "autohoot"))
end

@func_loader eingraph.expr_generator rewrite_einsum_expr
@func_loader gops.graph_transformer optimize simplify fuse_einsums
@func_loader gops.graph_inv_optimizer optimize_inverse prune_inv_node
@func_loader gops.optimal_tree generate_optimal_tree split_einsum generate_optimal_tree_w_constraint
@func_loader gops.graph_als_optimizer generate_sequential_optimal_tree
@func_loader gops.graph_dedup dedup

end


text:
_hlp_bathymetry_allowed_nonlinear_modes =[:None, :PiecewiseLinear, :Cone, :Nonconvex]
function _hlp_bathymetry_allowed_nonlinear_mode(mode::Symbol)
  if ! in(mode, _hlp_outputs_allowed_nonlinear_modes)
    error("Given nonlinear mode not recognised: " * string(mode) * ". " *
          "Allowed modes are: " * string(_hlp_bathymetry_allowed_nonlinear_modes) * ".")
  end
end

"""
Implements the given bathymetry relationship between the reservoir level and the volume in the reservoir.
"""
function hlp_JuMP_bathymetry(model::JuMP.Model, volume::JuMP.Variable, level::JuMP.Variable, r::Reservoir)
  g = bathymetry(r)
  if isempty(g)
    error("The reservoir has no bathymetry!")
  end
  return hlp_JuMP_bathymetry(model, volume, level, bathymetry(r))
end

function hlp_JuMP_bathymetry(model::JuMP.Model, volume::JuMP.Variable, level::JuMP.Variable, g::ReservoirBathymetry)
  error("[hlp_JuMP_bathymetry] Helper not defined for " * string(typeof(g)) * ".")
end

function hlp_JuMP_bathymetry(model::JuMP.Model, volume::JuMP.Variable, level::JuMP.Variable, g::LinearReservoirBathymetry;
                           nonlinear_mode::Symbol=:None)
  # Ignore nonlinear_mode: not relevant for this bathymetry (purely linear).
  @constraint(model, level == intercept(g) + slope(g) * volume)
end

function hlp_JuMP_bathymetry(model::JuMP.Model, volume::JuMP.Variable, level::JuMP.Variable, g::LinearReservoirBathymetry;
                           nonlinear_mode::Symbol=:None)
  _hlp_bathymetry_allowed_nonlinear_mode(nonlinear_mode)

  if mode == :Nonconvex
    @NLconstraint(m, level == coefficient(g, 0) + coefficient(g, 1) * volume + coefficient(g, 2) * storage ^ 2)
  else
    error("[hlp_JuMP_bathymetry] TODO: mode " * string(nonlinear_mode) * " not implemented for LinearReservoirBathymetry.")
  end
end


text:
function multisolve(m::Model, mdata::MultiData, ::WeightedSum, ::LinearProblem)
    objectives = mdata.objectives
    sensemap = Dict(:Min => 1.0, :Max => -1.0)

    vararr = [JuMP.Variable(m,i) for i in 1:MathProgBase.numvar(m)]

    numobj = length(objectives)
    Phi = zeros(numobj,numobj)

    # Individual minimisations
    for (i, objective) in enumerate(objectives)
        @objective(m, objective.sense, objective.f)
        for (key, value) in objective.initialvalue
            # TODO: What is the correct way to set these values?
            setvalue(m.objDict[key], value)
        end
        status = solve(m, ignore_solve_hook=true);
        if status != :Optimal
            return status
        end

        push!(mdata.utopiavarvalues, [getvalue(var) for var in vararr])
        push!(mdata.paretovarvalues, [getvalue(var) for var in vararr])

        Phi[:,i] = sensevalue.(objectives)

        push!(mdata.paretofront, getvalue(objectives))
    end
    Fmax = maximum(Phi, dims=2)
    Fmin = minimum(Phi, dims=2) # == diag(Phi)?
    if Fmax == Fmin
        error("The Nadir and Utopia points are equal") # I think that's what this means?
    end

    mdata.Phi = Phi

    beta = zeros(numobj); beta[end] = 1.0

    betatree = betas(numobj, mdata.pointsperdim-1)

    for betaval in betatree
        if count(t -> t != 0, betaval) == 1
            # Skip individual optimisations as
            # they are already performed
            continue
        end
        JuMP.setobjective(m, :Min,
            sum(betaval[i]*(sensemap[objectives[i].sense]*objectives[i].f -
            Fmin[i])/(Fmax[i]-Fmin[i]) for i in Base.OneTo(numobj))
        )

        status = solve(m, ignore_solve_hook=true);
        if status != :Optimal
            return status
        end

        push!(mdata.paretofront, getvalue(objectives))
        push!(mdata.paretovarvalues, [getvalue(var) for var in vararr])
    end
    return :Optimal
end

function multisolve(m::Model, mdata::MultiData, ::EpsilonCons, ::LinearProblem)
    objectives = mdata.objectives
    sensemap = Dict(:Min => 1.0, :Max => -1.0)

    vararr = [JuMP.Variable(m,i) for i in Base.OneTo(MathProgBase.numvar(m))]

    numobj = length(objectives)
    Phi = zeros(numobj,numobj)

    if numobj > 2
        # TODO:
        # The logic  here becomes difficult, as the feasible region will require
        # a dependency between the constraints
        Base.error(":EPS is thought through for > 2 objectives yet")
    end

    # Individual minimisations
    for (i, objective) in enumerate(objectives)
        @objective(m, objective.sense, objective.f)
        for (key, value) in objective.initialvalue
            # TODO: What is the correct way to do this? (anonymous variables)
            setvalue(m.objDict[key], value)
        end
        status = solve(m, ignore_solve_hook=true);
        if status != :Optimal
            return status
        end

        push!(mdata.utopiavarvalues, [getvalue(var) for var in vararr])
        push!(mdata.paretovarvalues, [getvalue(var) for var in vararr])

        Phi[:,i] = sensevalue.(objectives)

        push!(mdata.paretofront, getvalue(objectives))
    end
    Fmax = maximum(Phi, dims=2)
    Fmin = minimum(Phi, dims=2) # == diag(Phi)?

    mdata.Phi = Phi
    # maintain last objective, other as constraints
    @objective(m, objectives[end].sense, objectives[end].f)

    @constraint(m, objconstr,
                  sensemap[objectives[1].sense]*objectives[1].f
                  <= Fmax[1])
    betatree = betas(numobj, mdata.pointsperdim-1)

    for betaval in 0.0:0.001:1.0

        JuMP.setRHS(objconstr, betaval*Fmin[1]+(1-betaval)*Fmax[1])

        status = solve(m, ignore_solve_hook=true);
        if status != :Optimal
            return status
        end

        push!(mdata.paretofront, getvalue(objectives))
        push!(mdata.paretovarvalues, [getvalue(var) for var in vararr])
    end

    return :Optimal
end


text:
using SymPy
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

@testset "Logical" begin
    ## More logical expressions
    @vars x y
    (x ≪ 0) ∧ Le(x*y,1) ∨ (x ⩵ y) ∧ (¬(x ≫ 3))
end


text:
# Autogenerated wrapper script for FLAC_jll for x86_64-apple-darwin
export flac, libflac, libflacpp, metaflac

using Ogg_jll
JLLWrappers.@generate_wrapper_header("FLAC")
JLLWrappers.@declare_library_product(libflac, "@rpath/libFLAC.8.dylib")
JLLWrappers.@declare_library_product(libflacpp, "@rpath/libFLAC++.6.dylib")
JLLWrappers.@declare_executable_product(flac)
JLLWrappers.@declare_executable_product(metaflac)
function __init__()
    JLLWrappers.@generate_init_header(Ogg_jll)
    JLLWrappers.@init_library_product(
        libflac,
        "lib/libFLAC.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libflacpp,
        "lib/libFLAC++.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        flac,
        "bin/flac",
    )

    JLLWrappers.@init_executable_product(
        metaflac,
        "bin/metaflac",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
n = 2

mass = 1.0
gravity = 9.81
mass_matrix = Diagonal(mass * ones(n))
bias = [0.0; gravity * mass]
contact_jacobian = Diagonal(ones(n))
ϕ(q) = q[2]

h = 0.01

a = 500.0
softplus(x) = 1 / a * log(1.0 + exp(a * x))
softminus(x) = x - softplus(x)

plot(softplus)
plot(softminus)

function dynamics(x, u, t)
	q = x[1:n]
	v = x[n .+ (1:n)]

	k_impact = 500.0
	b_impact = 10.0
	# λ = [0.0; -k_impact * min(0.0, ϕ(q))]
	λ = [0.0; -k_impact * softminus(ϕ(q)) + max(0.0, - b_impact * v[2])]

	v⁺ = v + h * (mass_matrix \ (-bias + transpose(contact_jacobian) * λ))
	q⁺ = q + h * v⁺

	return [q⁺; v⁺]
end

q1 = [0.0; 1.0]
v1 = [0.0; 0.0]
x1 = [q1; v1]

x_hist = [x1]
T = 1000

for t = 1:T-1
	push!(x_hist, dynamics(x_hist[end], nothing, t))
	println(ϕ(x_hist[end][1:2]))
end

plot(hcat(x_hist...)[1:2, :]', )

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using OrdinaryDiffEq
using Trixi
using Plots

initial_condition = initial_condition_constant

surface_flux = FluxPlusDissipation(flux_chandrashekar, DissipationLocalLaxFriedrichs(max_abs_speed_naive)) # flux_lax_friedrichs #
volume_flux  = flux_chandrashekar
volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
basis = LobattoLegendreBasis(4)
solver = DGSEM(basis, surface_flux, volume_integral)

cells_per_dimension = (16, 16)
CFL = 0.5
tspan = (0.0, 50)
###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

# mapping O-mesh
function mapping(xi_, eta_)

    ξ = xi_
    η = eta_

    x = 100 * (2 + ξ ) * cos(π * (η + 1))
    y = 100 * (2 + ξ ) * sin(π * (η + 1))

    return SVector(x, y)
end

function boundary_condition_constant_farfield(u_inner, orientation, direction, x, t,
    surface_flux_function,
    equations::CompressibleEulerEquations2D)
    # Far Field Conditions
    u_boundary = SVector(1, 0.1, -0.2, 10)

    # Calculate boundary flux
    if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
        flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
    else # direction == 4 # u_boundary is "left" of boundary, u_inner is "right" of boundary
        flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
    end

    return flux
end

function boundary_condition_freeslip(u_inner, orientation, direction, x, t,
    surface_flux_function,
    equations::CompressibleEulerEquations2D)
    # Freeslip wall Conditions
    rho, v1, v2, p = cons2prim(u_inner, equations)
    u_boundary = prim2cons(SVector(rho, -v1, -v2, p), equations)

    # Calculate boundary flux
    if direction == 2 # u_inner is "left" of boundary, u_boundary is "right" of boundary
        flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
    else # direction == 4 # u_boundary is "left" of boundary, u_inner is "right" of boundary
        flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
    end

    return flux
end

mesh = CurvedMesh(cells_per_dimension, mapping,periodicity=(false,true))

#mesh = CurvedMesh(cells_per_dimension, coordinates_min, coordinates_max, periodicity=(false, true))

# Implement boundary conditions
# x neg and x_pos are not periodic
# boundary_conditions = boundary_condition_periodic
boundary_conditions = (x_neg=boundary_condition_constant_farfield,
                       x_pos=boundary_condition_constant_farfield,
                       y_neg=boundary_condition_periodic,
                       y_pos=boundary_condition_periodic)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=CFL)
analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
                                                 extra_analysis_errors=(:conservation_error,),
                                                 extra_analysis_integrals=(entropy, energy_total,
                                                                           energy_kinetic, energy_internal)

                                                )

callbacks = CallbackSet(summary_callback, stepsize_callback, analysis_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

pd=PlotData2D(sol)
b=plot(pd["rho"])
plot!(getmesh(pd))

# plot(sol)
# savefig(plotd, "test.png")

max_stars_repo_path: src/WR2/MS4 Freestream.jl

max_issues_repo_path: src/WR2/MS4 Freestream.jl

max_forks_repo_path: src/WR2/MS4 Freestream.jl


text:
export SpiralTrajectory, spiralNodes, spiralDensity

"""
    SpiralTrajectory(numProfiles, numSamplingPerProfile
                  ; TE::Float64=0.0
                  , AQ::Float64=1.e-3
                  , windings::Real= 6.25
                  , angleOffset= :equispaced
                  , kargs...)

returns a 2d spiral trajectory.

# Arguments
* `numProfiles::Int64`            - number of profiles
* `numSamplingPerProfile::Int64`  - number of sampling points per profile
* (`TE::Float64=0.0`)             - echo time in s
* (`AQ::Float64=1.e-3`)           - readout duration in s (per profile)
* (`windings::Real= 6.25`)        - number of windings of the spiral profiles
* (`angleOffset= :equispaced`)    - spacing of profile angles (`:equispaced` sampling, `:golden` angle sampling or `:random` sampling)
"""
function SpiralTrajectory(::Type{T}, numProfiles, numSamplingPerProfile
                  ; TE=0.0
                  , AQ=1.e-3
                  , windings::Real= 6.25
                  , kmax::Real=0.5
                  , angleOffset::String="equispaced"
                  , kargs...) where T
  nodes = spiralNodes(T, numProfiles, numSamplingPerProfile; windings=windings,
                      angleOffset=angleOffset, kmax=kmax)
  times = readoutTimes(T, numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ)
  return  Trajectory("Spiral", nodes, times, TE, AQ, numProfiles, numSamplingPerProfile, 1, false, true)
end

function spiralNodes(::Type{T}, numProfiles::Int64
                            , numSamplingPerProfile::Int64
                            ; windings::Real= 6.25
                            , kmax::Real=0.5
                            , angleOffset::String="equispaced"
                            , kargs...) where T

  nodes = zeros(T, 2, numSamplingPerProfile, numProfiles)
  if angleOffset == "golden"
      angles = [i*(3-sqrt(5))/2  for i=0:numProfiles-1 ]
  elseif angleOffset == "random"
      angles = sort(rand(numProfiles))
  elseif angleOffset == "equispaced"
      angles = collect((0:numProfiles-1)/numProfiles)
  end

  w = windings # 8/64 * 50 = 6.25 which means we have 6.25 turns aka windings
  for l = 1:numProfiles
    for k = 1:numSamplingPerProfile
      # For the special case that τ(t)=t, the amount of time spent for each
      # winding is constant, regardless of whether the acquired
      # winding is near the center or in the outer part of the spiral. I
      # other words, the readout gradients reach their maximum
      # performance at the end of the acquisition.
      t = sqrt((k-1)/(numSamplingPerProfile-1)) #
      nodes[1,k,l] = kmax*t*cos(2*pi*( w*t + angles[l] ))
      nodes[2,k,l] = kmax*t*sin(2*pi*( w*t + angles[l] ))
    end
  end
  return reshape(nodes, 2, numSamplingPerProfile*numProfiles)
end

"""
Calculating the density using VoronoiCells
"""
function spiralDensity(::Type{T}, numProfiles::Int64, numSamplingPerProfile::Int64) where T
  @error "Not implemented!"
end


text:
using Test

if Sys.iswindows()
    @info "No tests enabled for for your platform by default."
else
    @testset "all thrift tests" begin
        @info("Generating code")
        include("gen.jl")
        @info("Running client server tests")
        include("clntsrvr.jl")

        @info("Running protocol and transport tests")
        include("memtransport_tests.jl")
        include("filetransport_tests.jl")
        include("utils_tests.jl")
    end
end


text:
import Base: getindex, setindex!, firstindex, lastindex, size

mutable struct Box{T}<:AbstractMatrix{T}
    values::T
    xindices::UnitRange{Int}
    yindices::UnitRange{Int}
    middle::Tuple{Int, Int}
    function Box(matrix::AbstractMatrix, startx::V, stopx::V, starty::V, stopy::V) where V <: Integer
        @assert startx > 0 && stopx < size(matrix, 1)
        @assert starty > 0 && stopy < size(matrix, 2)
        mid = (middle(startx, stopx), middle(starty, stopy))
        xindices = startx:stopx
        yindices = starty:stopy
        new(view(matrix, xindices, yindices), xindices, yindices, middle)
    end
end

function Box(matrix::AbstractMatrix, indices)
    Box(matrix, indices[1]..., indices[2]...)
end

size(box::Box) = size(box.values)
getindex(box::Box, i) = getindex(box.values, i)
setindex!(box::Box, i, v) = setindex!(box.values, i, v)
firstindex(box::Box) = firstindex(box.values)
lastindex(box::Box) = lastindex(box.values)

middle(start, stop) = (stop-start)/2 + start |> x -> round(Int64, x)
function boxlim(start, stop, width)
    lower = middle(start, stop) - width/2
    upper = middle(start, stop) + width/2
    (lower, upper) .|> x -> floor(Int, x)
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
using DrMZ
using Documenter
# using DocumenterMarkdown

# push!(LOAD_PATH,"../src/")

DocMeta.setdocmeta!(DrMZ, :DocTestSetup, :(using DrMZ); recursive=true)

makedocs(;
    modules=[DrMZ],
    authors="Brek Meuris",
    repo="https://github.com/brekmeuris/DrMZ.jl/blob/{commit}{path}#{line}",
    sitename="DrMZ.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://brekmeuris.github.io/DrMZ.jl/",
        assets=String[],
    ),
    # format = Markdown(),
    pages=[
        "Home" => "index.md",
    ],
)

# FOR HTML
deploydocs(;
    repo="github.com/brekmeuris/DrMZ.jl",
    devbranch="main",
    target="build",
#     # devurl="docs",
    branch="gh-pages",
)


text:
###################################################
# Exemplo de uso do CPLEX/GLPK no Julia
#
# * CPLEX requer instalação e licença separada.
#   O pacote CPLEX.jl é apenas uma interface.
#
# * Ao instalar GLPK.jl no Julia, GLPK será
#   automaticamente instalado.
#
# Autor: Leonardo D. Secchin
# Data : 16/07/2021
# leonardosecchin.github.io/topicospo/
###################################################

using JuMP

# JuMP é um pacote para escrever modelos de otimização.
# Você pode associar diferentes solvers ao modelo.
# Para usar a interface com o CPLEX, você precisa instalá-lo
# em sua máquina, com uma licença válida. Caso não consiga,
# uma opção gratuita é o GLPK (https://www.gnu.org/software/glpk).
# Descomente as linhas adequadas para o solver que irá utilizar.

using CPLEX
# using GLPK

# INSTÂNCIA UFL DA LISTA DE EXERCÍCIOS, TÓPICO "RELAXAÇÃO LAGRANGIANA"
m = 6
n = 5

f = [2;4;5;3;3]

c = [6  2  1  3  5;
	4  10  2  6  1;
	3  2  4  1  3;
	2  0  4  1  4;
	1  8  6  2  5;
	3  2  4  8  1]

########## UFL (formulação forte) ##########

# estrutura JuMP associada com CPLEX
UFL = Model(CPLEX.Optimizer)

# estrutura JuMP associada com GLPK
# UFL = Model(GLPK.Optimizer)

set_silent(UFL)  # desliga saída solver na tela

# insere variáveis
@variable(UFL, x[1:m,1:n] >= 0)
@variable(UFL, y[1:n], Bin)

# função objetivo
@objective(UFL, Min, sum(c[i,j]*x[i,j] for i in 1:m, j in 1:n) + sum(f[j]*y[j] for j in 1:n))

# restrições
for i in 1:m
    @constraint(UFL, sum(x[i,j] for j in 1:n) == 1)
end

for i in 1:m, j in 1:n
    @constraint(UFL, x[i,j] <= y[j])
end

# Resolve UFL
optimize!(UFL)

# FO ótima UFL
print("UFL Resolvido! FO ótima UFL: ")
println(objective_value(UFL))

max_stars_repo_name: leonardosecchin/temp

max_issues_repo_name: leonardosecchin/temp

max_forks_repo_name: leonardosecchin/temp


text:
# MatrixReading.jl
#
# File for reading a sparse SPAM dataset into julia.

module MatrixReading

export ReadMatrix;

# (sp_matrix, tokenlist, category) = ReadMatrix(filename::ASCIIString)
#
# Reads the file stored at 'filename,' which is of the format of
# MATRIX.TEST, and returns a 3-tuple. The first part is 'sp_matrix',
# an m-by-n sparse matrix, where m is the number of training/testing
# examples and n is the dimension, and each row of sp_matrix consists
# of counts of word appearances. (So sp_matrix[i, j] is the number of
# times word j appears in document i.)
#
# tokenlist is a list of the words, where tokenlist[1] is the first
# word in the dictionary and tokenlist[end] is the last.
#
# category is a {0, 1}-valued vector of positive and negative
# examples. Before using in SVM code, you should transform categories
# to have signs +/-1.
function ReadMatrix(filename::ASCIIString)
  fstream = open(filename);
  # Read header line, discard
  headerline = readline(fstream);
  # Read rows and columns, turn into integers
  row_col_line = readline(fstream);
  row_col_split = split(row_col_line);
  num_rows = parse(Int, row_col_split[1]);
  num_cols = parse(Int, row_col_split[2]);
  # Read the list of tokens - just a long string!
  tokenlist = readline(fstream);
  tokenlist = split(tokenlist);

  # Now to read the matrix into the matrix. Each row represents a
  # document (mail), each column represents a distinct token. As the
  # data isn't actually that big, we just use a full matrix to save
  # time.
  full_mat = zeros(num_rows, num_cols);
  categories = zeros(num_rows);
  for ii = 1:num_rows
    str_inds = readline(fstream);
    str_split = split(str_inds);
    categories[ii] = parse(Int, str_split[1]);
    jj = 2;
    ind = 1;
    while (str_split[jj] != "-1")
      ind = parse(Int, str_split[jj]) + ind;
      count = parse(Int, str_split[jj + 1]);
      full_mat[ii, ind] = count;
      jj = jj + 2;
    end
  end
  close(fstream);
  return (sparse(full_mat), tokenlist, categories);
end

end   # module MatrixReading