text:
using ClimateMachine.Mesh.Elements
using GaussQuadrature
using Test

@testset "GaussQuadrature" begin
    for T in (Float32, Float64, BigFloat)
        let
            x, w = GaussQuadrature.legendre(T, 1)
            @test iszero(x)
            @test w ≈ [2 * one(T)]
        end

        let
            endpt = GaussQuadrature.left
            x, w = GaussQuadrature.legendre(T, 1, endpt)
            @test x ≈ [-one(T)]
            @test w ≈ [2 * one(T)]
        end

        let
            endpt = GaussQuadrature.right
            x, w = GaussQuadrature.legendre(T, 1, endpt)
            @test x ≈ [one(T)]
            @test w ≈ [2 * one(T)]
        end

        let
            endpt = GaussQuadrature.left
            x, w = GaussQuadrature.legendre(T, 2, endpt)
            @test x ≈ [-one(T); T(1 // 3)]
            @test w ≈ [T(1 // 2); T(3 // 2)]
        end

        let
            endpt = GaussQuadrature.right
            x, w = GaussQuadrature.legendre(T, 2, endpt)
            @test x ≈ [T(-1 // 3); one(T)]
            @test w ≈ [T(3 // 2); T(1 // 2)]
        end
    end

    let
        err = ErrorException("Must have at least two points for both ends.")
        endpt = GaussQuadrature.both
        @test_throws err GaussQuadrature.legendre(1, endpt)
    end

    let
        T = Float64
        n = 100
        endpt = GaussQuadrature.both

        a, b = GaussQuadrature.legendre_coefs(T, n)

        err = ErrorException(
            "No convergence after 1 iterations " * "(try increasing maxits)",
        )

        @test_throws err GaussQuadrature.custom_gauss_rule(
            -one(T),
            one(T),
            a,
            b,
            endpt,
            1,
        )
    end
end

@testset "Operators" begin
    P5(r::AbstractVector{T}) where {T} =
        T(1) / T(8) * (T(15) * r - T(70) * r .^ 3 + T(63) * r .^ 5)

    P6(r::AbstractVector{T}) where {T} =
        T(1) / T(16) *
        (-T(5) .+ T(105) * r .^ 2 - T(315) * r .^ 4 + T(231) * r .^ 6)
    DP6(r::AbstractVector{T}) where {T} =
        T(1) / T(16) *
        (T(2 * 105) * r - T(4 * 315) * r .^ 3 + T(6 * 231) * r .^ 5)

    IPN(::Type{T}, N) where {T} = T(2) / T(2 * N + 1)

    N = 6
    for test_type in (Float32, Float64, BigFloat)
        r, w = Elements.lglpoints(test_type, N)
        D = Elements.spectralderivative(r)
        x = LinRange{test_type}(-1, 1, 101)
        I = Elements.interpolationmatrix(r, x)

        @test sum(P5(r) .^ 2 .* w) ≈ IPN(test_type, 5)
        @test D * P6(r) ≈ DP6(r)
        @test I * P6(r) ≈ P6(x)
    end

    for test_type in (Float32, Float64, BigFloat)
        r, w = Elements.lgpoints(test_type, N)
        D = Elements.spectralderivative(r)

        @test sum(P5(r) .^ 2 .* w) ≈ IPN(test_type, 5)
        @test sum(P6(r) .^ 2 .* w) ≈ IPN(test_type, 6)
        @test D * P6(r) ≈ DP6(r)
    end
end


text:
# ----------------------------------------------------------
# Space mappings
# ----------------------------------------------------------

function map_device_to_view(world::Ranger.World,
    dx::Int32, dy::Int32, view_point::Geometry.Point{Float64})
    Math.transform!(world.inv_view_space, Float64(dx), Float64(dy), view_point);
end

# intermediate view-space node.
view_point = Geometry.Point{Float64}()
# Standin signalling no psuedoRoot.
nil_node = Nodes.NodeNil()

function map_device_to_node!(world::Ranger.World,
    dx::Int32, dy::Int32, node::Ranger.AbstractNode, local_point::Geometry.Point{Float64})
    # Mapping from device to node requires to transforms from to "directions"
    # 1st is upwards transform and the 2nd is downwards transform.

    # Upwards from node to world-space (aka view-space)
    wtn = world_to_node_transform(node, nil_node)
    # downwards from device-space to view-space
    map_device_to_view(world, dx, dy, view_point)

    # Now map view-space point to local-space of node
    Math.transform!(wtn, view_point.x, view_point.y, local_point)
    # Optional scaling
    # local_point.x *= node.transform.scale.x
    # local_point.y *= node.transform.scale.y
end

function world_to_node_transform(node::Ranger.AbstractNode, psuedoRoot::Ranger.AbstractNode)
    wtn = node_to_world_transform(node, psuedoRoot)
    Math.invert!(wtn)
    wtn
end

function node_to_world_transform(node::Ranger.AbstractNode, psuedoRoot::Ranger.AbstractNode)
    # println("*****************************")
    aft = Nodes.calc_transform!(node)
    # println("node: ", node, ", aft: ", aft)

    # A transform to accumulate the parent transforms.
    comp = Math.copy(aft)

    out = Math.AffineTransform{Float64}()

    # Iterate "upwards" starting with the child towards the parents
    # starting with this child's parent.
    p = node.base.parent
    while (!Nodes.is_nil(p))
        # println("p: ", p, ", comp: ", comp)
        parentT = Nodes.calc_transform!(p)
        # println("parentT: ", parentT)

        # Because we are iterating upwards we need to pre-multiply each child.
        # Ex: [child] x [parent]
        # ----------------------------------------------------------
        #           [comp] x [parentT]
        #               |
        #               | out
        #               v
        #             [comp] x [parentT]
        #                 |
        #                 | out
        #                 v
        #               [comp] x [parentT...]
        #
        # This is a pre-multiply order
        # [child] x [parent of child] x [parent of parent of child]...
        #
        # In other words the child is mutiplied "into" the parent.
        Math.multiply!(comp, parentT, out)
        Math.set!(comp, out)

        if (p == psuedoRoot)
            println("Hit psuedoRoot")
            break;
        end

        # next parent upwards
        p = p.base.parent
    end

    comp
end

ext: jl

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
struct System
    elems::Vector{Element}
    conss::Vector{AbstractString}
    phass::Vector{Phase}
    nelem::Int # number of elements
    ncons::Int # number of constituents
    nphas::Int # number of phases
    nlatt::Int # maximum number of sublattices in a phase
    nsite::Matrix{Real} # number of sites in each sublattice
    temp::Tuple{Real,Real} # temprature range where all the functions are valid
    model::Model # JuMP model
end

function init_system(db::Database, elems::Vector{Element}, phass::Vector{<:Phase})
    # We do this because we modify Phase structs destructively
    phass = deepcopy(phass)

    # Evaluate parsed CALPHAD functions and define them as Julia functions and
    # determine minimum and maximum temperature simultaneously.
    Tl, Tu = -Inf, Inf
    for func in db.funcs
        func.name == "R" && continue
        julialize!(func, db.funcs)
        eval(Meta.parse(func.funcstr))
        Tl = func.temp[1] > Tl ? func.temp[1] : Tl
        Tu = func.temp[2] < Tu ? func.temp[2] : Tu
    end
    for phas in phass
        for param in phas.params
            julialize!(param, vcat(db.funcs, phas.params))
            eval(Meta.parse(param.funcstr))
            Tl = param.temp[1] > Tl ? param.temp[1] : Tl
            Tu = param.temp[2] < Tu ? param.temp[2] : Tu
        end
    end
    temp = (Tl, Tu)

    # The elements always exists as consitituents
    I = length(elems) # number of components
    conss = [ elems[i].name for i in 1:I ]

    # Register non-element constituents
    for phas in phass
        for param in phas.params
            for latt in param.comb
                for cons in latt
                    !(cons in conss) && push!(conss, cons)
                end
            end
        end
    end

    J = length(conss) # number of constituents
    K = length(phass) # number of phases

    # Maximum number of sublattices in a phase
    S = maximum([ length(phass[k].cons) for k in 1:K ])

    # Indexing constituent
    cons_ids = Dict{AbstractString,Int}()
    for j in 1:J
        push!(cons_ids, conss[j] => j)
    end

    # Put zero or empty elements in each sublattice for later use
    for k in 1:K
        @assert length(phass[k].sites) == length(phass[k].cons)
        while length(phass[k].sites) < S
            push!(phass[k].sites, zero(Int))
            push!(phass[k].cons, Vector{Int}[])
        end
    end

    # Reconstuction of constitution vector into an array of Int
    constitution = [ map.(name -> get(cons_ids, name, 0), phass[k].cons) for k in 1:K ]

    # Chemical stoichiometry of constituents
    comps = [ stoichiometry(conss[j]) for j in 1:J ]

    # Stoichiometry matrix, a[k,s,j,i]
    a = [ j in constitution[k][s] ? get(comps[j], elems[i].name, 0) : 0
          for k in 1:K, s in 1:S, j in 1:J, i in 1:I ]

    # Number of sites on sublattice s in phase k, n[k,s]
    n = [ phass[k].sites[s] for k in 1:K, s in 1:S ]

    #
    # GLPK as the lower solver in EAGO
    #
    glpk = GLPK.Optimizer()
    MOI.set(glpk, MOI.RawParameter("meth"), 2)
    MOI.set(glpk, MOI.RawParameter("msg_lev"), 0)

    @debug begin
        MOI.set(glpk, MOI.RawParameter("msg_lev"), 4) # GLP_MSG_ALL
        MOI.set(glpk, MOI.RawParameter("out_frq"), 1)
        "Enabled GLPK debugging"
    end

    #
    # Ipopt as the upper solver in EAGO
    #
    ipopt = Ipopt.Optimizer()

    eps = 2e-8 # TODO: parameterize this

    MOI.set(ipopt, MOI.RawParameter("tol"), 1e-2)
    MOI.set(ipopt, MOI.RawParameter("dual_inf_tol"), 1e-6)
    MOI.set(ipopt, MOI.RawParameter("constr_viol_tol"), eps)
    MOI.set(ipopt, MOI.RawParameter("max_iter"), 1000)
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)

    @debug begin
        MOI.set(ipopt, MOI.RawParameter("print_level"), 5)
        "Enabled Ipopt debugging"
    end

    #
    # EAGO Optimizer
    #
    model = Model(optimizer_with_attributes(EAGO.Optimizer,
        "relaxed_optimizer" => glpk,
        "upper_optimizer" => ipopt,
        "verbosity" => 0,
        "output_iterations" => 1,
        "iteration_limit" => 3,
        "cp_depth" => 0,
        "obbt_depth" => 0,
        "dbbt_tolerance" => eps,
        "absolute_constraint_feas_tolerance" => eps,
    ))

    @debug begin
        set_optimizer_attributes(model, "verbosity" => 5)
        "Enabled EAGO debugging"
    end

    # We have to use a special function xlogx() provided by EAGO.
    EAGO.register_eago_operators!(model)

    #
    # Definition of variables.
    #
    # We define sparse arrays of variables and add constraints for zero elements,
    # instead of defining as many variables as needed, for the sake of programmability
    # and readability of the code, assuming the optimizer is very efficient.
    #
    # T: temperature
    # X[i]: molar amount of components
    # Y[k]: molar amount of each phase
    # x[k,i]: molar fraction of component i in phase k
    # y[k,s,j]: site fraction of a constituent j in sublattice s of phase k
    #
    @variable(model, Tl <= T <= Tu)
    @variable(model, X[i=1:I] >= 0)
    @variable(model, Y[k=1:K] >= eps)
    # @variable(model, eps <= x[k=1:K,i=1:I] <= 1)
    @variable(model, eps <= y[k=1:K,s=1:S,j=1:J] <= 1)

    # Sum of y[k,j,s] for all constituents in a sublattice equals to unity.
    for k in 1:K, s in 1:S
        if constitution[k][s] ≠ []
            @constraint(model, sum(y[k,s,j] for j in 1:J) == 1)
        end
    end

    # Fix site fractions of zero elements to zero
    for k in 1:K, s in 1:S, j in 1:J
        if !in(j, constitution[k][s])
            fix(y[k,s,j], 0, force=true)
        end
    end

    #
    # x[k,i]: molar fraction of component i in phase k
    #
    # TODO: We Assume that amounts of molecular like constituents are small and
    # there's no vacancy. This is because Using NLexpression in NLconstraint
    # results in not obtaining solutions for upper problems.
    #
    @expression(model, x[k=1:K,i=1:I],
                sum(n[k,s] * sum(a[k,s,j,i] * y[k,s,j]
                                 for j in 1:J
                                 if a[k,s,j,i] ≠ 0)
                    for s in 1:S if constitution[k][s] ≠ []) /
                sum(n[k,s] for s in 1:S))
                # TODO: Use this as the denominator
                # sum(n[k,s] * sum(a[k,s,j,i′] * y[k,s,j]
                #                  for j in 1:J, i′ in 1:I
                #                  if conss[j] ≠ "Va" && a[k,s,j,i′] ≠ 0)
                #     for s in 1:S if constitution[k][s] ≠ []))

    # Relationship between molar amount of components and phases
    for i in 1:I
        @constraint(model, sum(sum(n[k,s] for s in 1:S) * Y[k] * x[k,i] for k in 1:K) == X[i])
    end

    # Register all the functions for parameters
    for phas in phass
        for param in phas.params
            funcname = param.funcname
            funcsym = Symbol(funcname)
            # try
                func = getfield(Calmato, funcsym)
            # catch
            #     error(phas.name * param.name * funcname)
            # end
            # try
                register(model, funcsym, 1, func, autodiff=true)
            # catch
            #     @warn "Duplicated definition of $funcname in the database"
            # end
        end
    end

    function G_phas(k::Int; k_param::Int=k, disordered::Bool=false)
        Gs_param = JuMP.NonlinearExpression[]
        m = 0
        for param in phass[k_param].params
            !in(param.symbol, ['G', 'L']) && continue

            m += 1
            l = param.order

            # Replace element names with indexes
            comb = map.(name -> get(cons_ids, name, 0), param.comb)

            # Value of the parameter
            funcname = param.funcname
            funcsym = Symbol(funcname)
            @eval push!($Gs_param, @NLexpression($model, ($funcsym)($T)))

            # Solo constituents
            s_solo = findall(x -> length(x) == 1, comb)
            solo = [ (s, comb[s][1]) for s in s_solo ]
            for (s, j) in solo
                Gs_param[m] =
                @NLexpression(model, Gs_param[m] * ( disordered ? x[k,j] : y[k,s,j] ))
            end

            # Paired constituents
            s_duo = findall(latt -> length(latt) == 2, comb)
            @assert length(s_duo) ≤ 1
            length(s_duo) == 0 && continue
            s = s_duo[1]
            i = comb[s][1]
            j = comb[s][2]
            Gs_param[m] = disordered ?
            @NLexpression(model, Gs_param[m] * x[k,i] * x[k,j]) :
            @NLexpression(model, Gs_param[m] * y[k,s,i] * y[k,s,j])
            for ν in 1:l
                Gs_param[m] = disordered ?
                @NLexpression(model, Gs_param[m] * ( x[k,i] - x[k,j] )) :
                @NLexpression(model, Gs_param[m] * ( y[k,s,i] - y[k,s,j] ))
            end
        end

        G = m > 1 ? @NLexpression(model, sum(Gs_param[i] for i in 1:m)) : Gs_param[1]

        # Ideal entropy of configuration
        nonzero = [(s,j) for s in 1:S, j in 1:J if j in constitution[k][s] && length(constitution[k][s]) > 1]
        if !disordered && length(nonzero) > 1
            G = @NLexpression(model, G + R*T*sum(n[k,s] * xlogx(y[k,s,j]) for (s,j) in nonzero))
        end

        return G
    end

    # Gibbs energy contribution from each phase
    Gs_phas = JuMP.NonlinearExpression[]
    for k in 1:K
        phas = phass[k]

        k_do = disorder_id(db, phas)
        if k_do ≠ nothing # ordered phase
            push!(Gs_phas, G_phas(k, k_param=k_do, disordered=true)) # disordered part
            # TODO: + 1.0 is an adhock and unphysical parameter of "penalty" on the ordered phase.
            @eval $Gs_phas[$k] = @NLexpression($model, $Gs_phas[$k] + $(G_phas(k)) + 1.0)
            @eval $Gs_phas[$k] = @NLexpression($model, $Gs_phas[$k] - $(G_phas(k, disordered=true)))
        else
            @eval push!($Gs_phas, $(G_phas(k)))
        end
    end

    @NLobjective(model, Min, sum(Y[k] * Gs_phas[k] for k in 1:K))

    @debug model

    return System(elems, conss, phass, I, J, K, S, n, temp, model)
end

function init_system(db::Database)
    elems = Vector{Element}()
    phass = Vector{Phase}()

    for ph in db.phass
        if ph.state == 'G'
            @warn "Gaseous phase is not supported yet. Exclude $(ph.name)."
            continue
        end
        if any(latt -> "Va" in latt, ph.cons)
            @warn "$(ph.name) includes vacancies. Calculation may be inaccurate."
        end
        if any(latt -> any(spec -> length(stoichiometry(spec)) > 1, latt), ph.cons)
            @warn "$(ph.name) includes molecular-like constituents. Calculation may be inaccurate."
        end
        push!(phass, ph)
    end

    for el in db.elems
        el.name in ["/-", "Va"] && continue
        push!(elems, el)
    end

    return init_system(db, elems, phass)
end

function julialize!(arg::AbstractGFunction, funcs::Vector{<:AbstractGFunction})
    if isempty(arg.funcstr) # toml database
        @assert !isempty(arg.temps)
        @assert !isempty(arg.exprs)

        arg.funcstr *= "function $(arg.funcname)(T)::typeof(T)\n"
        N = length(arg.exprs)
        for i in 1:N
            expr = arg.exprs[i]
            arg.funcstr *= i == 1 ? "\tif " : "\telseif "
            arg.funcstr *= "$(arg.temps[i]) ≤ T ≤ $(arg.temps[i+1])\n"
            arg.funcstr *= "\t\t" * expr * "\n"
        end
        arg.funcstr *= "\tend\nend\n"

        for func in funcs
            arg.funcstr = replace(arg.funcstr, func.name => func.funcname)
        end
    end

    for pair in ["T*LN(T)" => "xlogx(T)",
                 "LN(" => "log(",
                 "T*ln(T)" => "xlogx(T)",
                 "ln(" => "log(",
                 "TlnT" => "xlogx(T)",
                 "**" => "^",
                 ".+" => ".0+",
                 ".-" => ".0-",
                 ".*" => ".0*",
                 "./" => ".0/"]
        arg.funcstr = replace(arg.funcstr, pair)
    end
    for func in funcs
        func.name == "R" && continue
        func.name == arg.name && continue
        reg = Regex(func.funcname * "(\\(T\\)){0,}(?=\\W)")
        arg.funcstr = replace(arg.funcstr, reg => "$(func.funcname)(T)")
    end
end

function disorder_id(db::Database, phas::Phase)
    code = phas.type
    for char in code
        i = findfirst(type -> type.code == char, db.types)
        isnothing(i) && continue
        tdef = db.types[i]
        args = split(tdef.args)
        length(args) < 3 && continue
        if args[2] == phas.name && args[3] == "DISORDER_PART"
            @assert args[1] == "AMEND_PHASE_DESCRIPTION"
            phasname = args[4]
            j = findfirst(phas -> phas.name == phasname, db.phass)
            @assert j ≠ nothing
            return j
        end
    end
    return nothing
end

function equiatom(sys::System)
    N = length(sys.elems)
    return fill(1 // N, N)
end

function Base.display(sys::System)
    println("System:")
    nelem = length(sys.elems)
    println("\tElements: $nelem")
    nphas = length(sys.phass)
    println("\tPhases: $nphas")
end

# TODO: Detailed output
Base.print(sys::System) = Base.display(sys)


text:
# Local Tangent Space Alignment (LTSA)
# ---------------------------
# Principal Manifolds and Nonlinear Dimension Reduction via Local Tangent Space Alignment,
# Zhang, Zhenyue; Hongyuan Zha (2004),  SIAM Journal on Scientific Computing 26 (1): 313–338.
# doi:10.1137/s1064827502419154.

"""
    LTSA{NN <: AbstractNearestNeighbors, T <: Real} <: AbstractDimensionalityReduction

The `LTSA` type represents a local tangent space alignment model constructed for `T` type data with a help of the `NN` nearest neighbor algorithm.
"""
struct LTSA{NN <: AbstractNearestNeighbors, T <: Real} <: AbstractDimensionalityReduction
    λ::AbstractVector{T}
    proj::Projection{T}
    nearestneighbors::NN
    component::AbstractVector{Int}
end

## properties
outdim(R::LTSA) = size(R.proj, 1)
eigvals(R::LTSA) = R.λ
neighbors(R::LTSA) = R.nearestneighbors.k
vertices(R::LTSA) = R.component

## show
summary(io::IO, R::LTSA) = print(io, "LTSA(outdim = $(outdim(R)), neighbors = $(neighbors(R)))")

## interface functions
"""
    fit(LTSA, data; k=12, maxoutdim=2, nntype=BruteForce)

Fit a local tangent space alignment model to `data`.

# Arguments
* `data`: a matrix of observations. Each column of `data` is an observation.

# Keyword arguments
* `k`: a number of nearest neighbors for construction of local subspace representation
* `maxoutdim`: a dimension of the reduced space.
* `nntype`: a nearest neighbor construction class (derived from `AbstractNearestNeighbors`)

# Examples
```julia
M = fit(LTSA, rand(3,100)) # construct LTSA model
R = transform(M)           # perform dimensionality reduction
```
"""
function fit(::Type{LTSA}, X::AbstractMatrix{T};
        k::Int=12, maxoutdim::Int=2, ɛ::Real=1.0, nntype=BruteForce) where {T<:Real}
    # Construct NN graph
    NN = fit(nntype, X, k)
    D, E = knn(NN, X)
    G, C = largest_component(SimpleWeightedGraph(adjmat(D,E)))
    XX = @view X[:, C]
    n = length(C)

    S = ones(k)./sqrt(k)
    B = spzeros(T, n,n)
    for i=1:n
        II = @view E[:,i]
        VX = view(XX, :, II)

        # re-center points in neighborhood
        μ = mean(VX, dims=2)
        δ_x = VX .- μ

        # Compute orthogonal basis H of θ'
        θ_t = svd(δ_x).V[:,1:maxoutdim]

        # Construct alignment matrix
        G = hcat(S, θ_t)
        B[II, II] .+= Diagonal(fill(one(T), k)) .- G*transpose(G)
    end

    # Align global coordinates
    λ, V = decompose(B, maxoutdim)
    return LTSA{nntype, T}(λ, transpose(V), NN, C)
end

"""
    transform(R::LTSA)

Transforms the data fitted to the local tangent space alignment model `R` into a reduced space representation.
"""
transform(R::LTSA) = R.proj


text:
struct SOS
	order::Int
	indices::Vector{Int}
	weights::Vector{Float64}
	SOS(order, indices, weights) = new(order, indices, weights)
end

mutable struct MPSModel <: NEOSModel
	solver::NEOSSolver
	xmlmodel::String
	last_results::String
	ncol::Int
	nrow::Int
	A::SparseMatrixCSC{Float64, Int}
	collb::Vector{Float64}
	colub::Vector{Float64}
	f::Vector{Float64}
	rowlb::Vector{Float64}
	rowub::Vector{Float64}
	sense::Symbol
	colcat::Vector{Symbol}
	sos::Vector{SOS}
	objVal::Float64
	reducedcosts::Vector{Float64}
	duals::Vector{Float64}
	solution::Vector{Float64}
	status::Symbol

	MPSModel(solver) = new(solver, "", "", 0, 0, sparse(Int[], Int[], Float64[]), Float64[], Float64[], Float64[], Float64[], Float64[], :Min, Symbol[], SOS[], 0.0, Float64[], Float64[], Float64[], NOTSOLVED)
end

LinearQuadraticModel{S}(s::NEOSSolver{S, :MPS}) = MPSModel(s)

function neos_writexmlmodel!(m::MPSModel)
    io = IOBuffer()
    print(io, "<MPS>")
    mps_writer_sos = [(s.order, s.indices, s.weights) for s in m.sos]
    MPSWriter.write(io, m.A, m.collb, m.colub, m.f, m.rowlb, m.rowub, m.sense, m.colcat, mps_writer_sos, Array{Float64}(0,0))
    print(io, "</MPS>")
    # Convert the model to MPS and add
    m.xmlmodel = replace(m.solver.template, r"<MPS>.*</MPS>"is, String(take!(io)))
end

function loadproblem!(m::MPSModel, A, collb, colub, f, rowlb, rowub, sense)
	# @assert length(collb) == length(colub) == length(f)
	m.ncol = length(f)
	m.nrow = length(rowlb)
	m.A = sparse(A)
	m.collb = collb
	m.colub = colub
	m.colcat = fill(:Cont, m.ncol)
	m.f = f
	m.rowlb = rowlb
	m.rowub = rowub
	m.sense = sense
end

function anyints(m::MPSModel)
	for i in m.colcat
		if i == :Int || i==:Bin
			return true
		end
	end
	return false
end

function parseresults!(m::MPSModel, job)
	m.status = SOVLERERROR
	parse_status!(m.solver, m)
	if m.status == OPTIMAL
		parse_objective!(m.solver, m)
		m.solution = zeros(m.ncol)
		parse_solution!(m.solver, m)
		if m.solver.provides_duals && length(m.sos) == 0 && !anyints(m)
			parse_duals!(m.solver, m)
		end
		if m.sense == :Max
			# Since MPS does not support Maximisation
			m.objVal = -m.objVal
		end
	end
end

function addsos1!(m::MPSModel, indices::Vector, weights::Vector)
	if !m.solver.solves_sos
		throw(NEOSSolverError("Special Ordered Sets of type I are not supported by $(typeof(m.solver)). Try a different solver instead"))
	end
	push!(m.sos, SOS(1, indices, weights))
end
function addsos2!(m::MPSModel, indices::Vector, weights::Vector)
	if !m.solver.solves_sos
		throw(NEOSSolverError("Special Ordered Sets of type II are not supported by $(typeof(m.solver)). Try a different solver instead"))
	end
	push!(m.sos, SOS(2, indices, weights))
end

function status(m::MPSModel)
	m.status
end

function getobjval(m::MPSModel)
	m.objVal
end

function getsolution(m::MPSModel)
	m.solution
end

function getreducedcosts(m::MPSModel)
	if m.solver.provides_duals
		return m.reducedcosts
	else
		warn("Reduced costs are not available from $(typeof(m.solver)). Try a different solver instead")
		return fill(NaN, m.ncol)
	end
end

function getconstrduals(m::MPSModel)
	if m.solver.provides_duals
		return 	m.duals
	else
		warn("Constraint duals are not available from $(typeof(m.solver)). Try a different solver instead")
		return fill(NaN, m.nrow)
	end
end

function getsense(m::MPSModel)
	m.sense
end

function setsense!(m::MPSModel, sense::Symbol)
	m.sense = sense
end

function setvartype!(m::MPSModel, t::Vector{Symbol})
	m.colcat = t
end


text:
module FEOperatorsFromTermsTests

using Test
using Gridap.Helpers
using Gridap.Arrays
using Gridap.Algebra
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Integration
using Gridap.Fields
using Gridap.FESpaces
using Gridap.CellData

import Gridap.Fields: ∇

u(x) = x[1] + x[2]
ufun_grad(x) = VectorValue(1.0,1.0)
∇(::typeof(u)) = ufun_grad
f(x) = -(3.0*x[1]+x[2]+1.0)

@law ν(u,x) = (u+1.0)*x[1]
@law dν(du,x) = x[1]*du

domain = (0,1,0,1)
partition = (4,4)
model = CartesianDiscreteModel(domain,partition)

order = 2
diritag = "boundary"
V = GradConformingFESpace(Float64,model,order,diritag)

U = TrialFESpace(V,u)

trian = get_triangulation(model)
degree = 3*order-1
quad = CellQuadrature(trian,degree)

const x = get_physical_coordinate(trian) # its safer to mark as constant objects that are captured

a(u,v,du) = inner( ∇(v), ν(u,x)*∇(du))
res(u,v) = a(u,v,u) - inner(v,f)
jac(u,v,du) = a(u,v,du) + inner(∇(v),dν(du,x)*∇(u))

t_Ω = FETerm(res,jac,trian,quad)
op = FEOperator(U,V,t_Ω)

uh = solve(op)

e = u - uh

l2(u) = inner(u,u)
sh1(u) = inner(∇(u),∇(u))
h1(u) = sh1(u) + l2(u)

el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))
ul2 = sqrt(sum( integrate(l2(uh),trian,quad) ))
uh1 = sqrt(sum( integrate(h1(uh),trian,quad) ))

@test el2/ul2 < 1.e-8
@test eh1/uh1 < 1.e-7

y = ones(num_free_dofs(U))
uh = FEFunction(U,y)
r = residual(op,uh)
A = jacobian(op,uh)
test_fe_operator(op,y,r,≈,jac=A)

end # module


text:
function michaelis_menten_model(concentration, Vm, K)
    (Vm * concentration) / (K + concentration)
end

# MicMen data is taken from the README of the NLReg.jl package: https://github.com/dmbates/NLreg.jl
# 1st column is concentration and 2nd column is the Rate
MicMenData = [
    0.02 76;
    0.02 47;
    0.06 97;
    0.06 107;
    0.11 123;
    0.11 139;
    0.22 159;
    0.22 152;
    0.56 191;
    0.56 201;
    1.1  207;
    1.1  200
];

MicMenConcentration = MicMenData[:, 1];
MicMenRate = MicMenData[:, 2];

# Fitness function takes a vector of Vm and K and calculates the RSS
function mic_men_fitness(params)
    Vm, K = params
    yhat = Float64[michaelis_menten_model(c, Vm, K) for c in MicMenConcentration]
    sumabs2(MicMenRate .- yhat)
end

using BlackBoxOptim
result = bboptimize(mic_men_fitness;
  SearchRange = (-1000.0, 1000.0), NumDimensions = 2, MaxSteps = 1e4)
Vm, K = best_candidate(result)
RSS = best_fitness(result)

println("NLReg.jl uses specific MicMen object and fit method to find:")
println("Vm = 212.684, K = 0.0641212, RSS = 1195.45")
println("\nUsing BlackBoxOptim.jl we find:")
println("Vm = $(Vm), K = $(K), RSS = $(RSS)")


text:
# Autogenerated wrapper script for MariaDB_Connector_ODBC_jll for i686-w64-mingw32
export libmaodbc

using iODBC_jll
using Libiconv_jll
using unixODBC_jll
using OpenSSL_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("MariaDB_Connector_ODBC")
JLLWrappers.@declare_library_product(libmaodbc, "mariadb_odbc.dll")
function __init__()
    JLLWrappers.@generate_init_header(iODBC_jll, Libiconv_jll, unixODBC_jll, OpenSSL_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libmaodbc,
        "bin\\mariadb\\mariadb_odbc.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    mutable struct Program_

Structure that contains all the pieces of code before they are finally assembled into the final evolve function.
"""
mutable struct Program_

    agent::Agent
    integrator::String
    integratorMedium::String
    neighbors::String
    declareVar::Expr
    declareF::Expr
    args::Array{Symbol,1}
    execInit::Expr
    execInloop::Expr
    execAfter::Expr
    returning::Expr

    update::Dict{String,Dict{Symbol,Int}}
end

function Program_(abm::Agent)
    return Program_(abm,"Euler","FTCS","full",quote end,quote end,Array{Symbol,1}(BASICARGS),
                    quote end,quote end,quote end,quote end,Dict{String,Dict{Symbol,Int}}())
end


text:
# [test/TimeSeriesLogger_tests.jl]
using Test

@testset "Time Series Logger Tests" begin

    using JLD2
    import ODEHybrid.ODE_SET, ODEHybrid.TimeSeriesLogger, ODEHybrid.add!, ODEHybrid.contains, ODEHybrid.get_log, ODEHybrid.initialize!, ODEHybrid.tsl_plot

    @testset "Add - Example" begin

        @load "solutions/tsl_add_noise.jld2" noise

        logger= TimeSeriesLogger();
        xs  = zeros(100, 2);
        xs[1, :] = [0.0; 0.0];
        t   = range(1, 100, step = 1);
        for i in t[2:end]
            xs[i, :] = xs[i-1, :] + noise[i, :]
            add!(logger, "walk", i, xs[i, :])
        end;
        x_jul = xs;
        t_jul = Array(t);

        # mat"""
        #     logger= TimeSeriesLogger();  % Make a new logger
        #     xs = zeros(100, 2);
        #     x   = [0.0; 0.0];
        #     xs(1, :) = x;

        #     for t = 2:100              % Generate a random walk
        #         x = x + $noise(t, :)';

        #         xs(t, :) = x';
        #         logger.add('walk', t, x); % Log sample x at time t
        #     end;

        #     $x_mat = xs;
        #     $t_mat = 1:100;
        #     % logger.plot()
        # """
        # t_mat = t_mat[1, :] # squeeze()

        @load "solutions/tsl_add_example.jld2" t_mat x_mat

        @test t_mat == t_jul
        @test x_mat == x_jul

        # tsl_plot(logger)

    end;

    @testset "Add - Multiple Measurements" begin

        @load "solutions/tsl_add_noise1.jld2" noise1;
        @load "solutions/tsl_add_noise2.jld2" noise2;

        logger= TimeSeriesLogger();
        x1_jul = zeros(100, 3);
        x2_jul = zeros(100, 3);
        x1_jul[1, :] = [0.0; 0.0; 0.0];
        x2_jul[1, :] = [0.0; 0.0; 0.0];
        t = range(1, 100, step = 1);

        for i in t[2:end]
            x1_jul[i, :] .= x1_jul[i - 1, :] .+ noise1[i, :]
            x2_jul[i, :] .= x2_jul[i - 1, :] .+ noise2[i, :]
            add!(logger, "walk", i, x1_jul[i, :])
            add!(logger, "run",  i, x2_jul[i, :])
        end
        t_jul = Array(t);

        # mat"""
        #     logger= TimeSeriesLogger();  % Make a new logger
        #     x1 = zeros(100, 3);
        #     x2 = zeros(100, 3);
        #     x1(1, :) = [0; 0; 0];
        #     x2(1, :) = [0; 0; 0];

        #     for t = 2:100
        #         x1(t, :) = x1(t - 1, :) + $noise1(t, :);
        #         x2(t, :) = x2(t - 1, :) + $noise2(t, :);
        #         logger.add('walk', t, x1(t, :));
        #         logger.add('run',  t, x2(t, :));
        #     end;

        #     $x1_mat = x1;
        #     $x2_mat = x2;
        #     $t_mat = 1:100;
        #     % logger.plot();
        # """
        # t_mat = t_mat[1, :] # squeeze()

        @load "solutions/tsl_add_multiple.jld2" x1_mat x2_mat t_mat

        @test x1_mat == x1_jul
        @test x2_mat == x2_jul
        @test t_mat == t_jul

        # tsl_plot(logger)
    end;

    @testset "Contains" begin

        @load "solutions/tsl_add_noise.jld2" noise

        logger = TimeSeriesLogger()
        xs  = zeros(100, 2)
        xs[1, :] = [0.0; 0.0]
        t   = range(1, 100, step = 1)
        for i in t[2:end]
            xs[i, :] = xs[i-1, :] + noise[i, :]
            add!(logger, "walk", i, xs[i, :])
        end
        x_jul = xs
        t_jul = Array(t)

        @test contains(logger, "walk")
        @test !contains(logger, "run")

        # mat"""
        #     logger= TimeSeriesLogger();  % Make a new logger
        #     xs = zeros(100, 2);
        #     x   = [0.0; 0.0];
        #     xs(1, :) = x;

        #     for t = 2:100              % Generate a random walk
        #         x = x + $noise(t, :)';

        #         xs(t, :) = x';
        #         logger.add('walk', t, x); % Log sample x at time t
        #     end;

        #     $m1 = logger.contains("walk");
        #     $m2 = logger.contains("run");
        # """
        # @test m1
        # @test !m2

    end;

    @testset "Get Log" begin

        @load "solutions/tsl_add_noise.jld2" noise

        logger= TimeSeriesLogger()
        xs  = zeros(100, 2)
        xs[1, :] = [0.0; 0.0]
        t   = range(1, 100, step = 1)
        for i in t[2:end]
            xs[i, :] = xs[i-1, :] + noise[i, :]
            add!(logger, "walk", i, xs[i, :])
        end

        t_jul, x_jul = get_log(logger, "walk")

        # mat"""
        #     logger= TimeSeriesLogger();  % Make a new logger
        #     xs = zeros(100, 2);
        #     x   = [0.0; 0.0];
        #     xs(1, :) = x;

        #     for t = 2:100              % Generate a random walk
        #         x = x + $noise(t, :)';

        #         xs(t, :) = x';
        #         logger.add('walk', t, x); % Log sample x at time t
        #     end;

        #     [$t_mat, $x_mat] = logger.get_log('walk');
        # """

        @load "solutions/tsl_get_log.jld2" t_mat x_mat

        @test t_jul == t_jul
        @test x_jul == x_mat
    end;

    @testset "Add - With Groups" begin

        @load "solutions/tsl_add_noise1.jld2" noise1
        @load "solutions/tsl_add_noise2.jld2" noise2

        logger= TimeSeriesLogger();
        x1_jul = zeros(100, 3);
        x2_jul = zeros(100, 3);
        x1_jul[1, :] = [0.0; 0.0; 0.0];
        x2_jul[1, :] = [0.0; 0.0; 0.0];
        t = range(1, 100, step = 1);

        for i in t[2:end]
            x1_jul[i, :] .= x1_jul[i - 1, :] .+ noise1[i, :]
            x2_jul[i, :] .= x2_jul[i - 1, :] .+ noise2[i, :]
            add!(logger, "walk", i, x1_jul[i, :], true, "1")
            add!(logger, "jog",  i, x1_jul[i, :], true, "2")
            add!(logger, "run",  i, x2_jul[i, :], true, "2")
        end
        t_jul = Array(t);

        # mat"""
        #     logger= TimeSeriesLogger();  % Make a new logger
        #     x1 = zeros(100, 3);
        #     x2 = zeros(100, 3);
        #     x1(1, :) = [0; 0; 0];
        #     x2(1, :) = [0; 0; 0];

        #     for t = 2:100
        #         x1(t, :) = x1(t - 1, :) + $noise1(t, :);
        #         x2(t, :) = x2(t - 1, :) + $noise2(t, :);
        #         logger.add('walk', t, x1(t, :), true, '1');
        #         logger.add('jog',  t, x1(t, :), true, '2');
        #         logger.add('run',  t, x2(t, :), true, '2');
        #     end;

        #     $x1_mat = x1;
        #     $x2_mat = x2;
        #     $t_mat = 1:100;
        #     % logger.plot();
        # """
        # t_mat = t_mat[1, :] # squeeze()

        @load "solutions/tsl_add_groups.jld2" t_mat x1_mat x2_mat

        @test x1_mat == x1_jul
        @test x2_mat == x2_jul
        @test t_mat == t_jul

        # tsl_plot(logger)
    end;

    @testset "Add - With show_it = false" begin

        @load "solutions/tsl_add_noise1.jld2" noise1
        @load "solutions/tsl_add_noise2.jld2" noise2

        logger= TimeSeriesLogger();
        x1_jul = zeros(100, 3);
        x2_jul = zeros(100, 3);
        x1_jul[1, :] = [0.0; 0.0; 0.0];
        x2_jul[1, :] = [0.0; 0.0; 0.0];
        t = range(1, 100, step = 1);

        for i in t[2:end]
            x1_jul[i, :] .= x1_jul[i - 1, :] .+ noise1[i, :]
            x2_jul[i, :] .= x2_jul[i - 1, :] .+ noise2[i, :]
            add!(logger, "walk", i, x1_jul[i, :], true, "1")
            add!(logger, "jog",  i, x1_jul[i, :], false, "2")
            add!(logger, "run",  i, x2_jul[i, :], true, "2")
        end
        t_jul = Array(t);

        # mat"""
        #     logger= TimeSeriesLogger();  % Make a new logger
        #     x1 = zeros(100, 3);
        #     x2 = zeros(100, 3);
        #     x1(1, :) = [0; 0; 0];
        #     x2(1, :) = [0; 0; 0];

        #     for t = 2:100
        #         x1(t, :) = x1(t - 1, :) + $noise1(t, :);
        #         x2(t, :) = x2(t - 1, :) + $noise2(t, :);
        #         logger.add('walk', t, x1(t, :), true, '1');
        #         logger.add('jog',  t, x1(t, :), false, '2');
        #         logger.add('run',  t, x2(t, :), true, '2');
        #     end;

        #     $x1_mat = x1;
        #     $x2_mat = x2;
        #     $t_mat = 1:100;
        #     % logger.plot();
        # """
        # t_mat = t_mat[1, :]; # squeeze()

        @load "solutions/tsl_add_no_show.jld2" x1_mat x2_mat t_mat

        @test x1_mat == x1_jul
        @test x2_mat == x2_jul
        @test t_mat == t_jul

        # tsl_plot(logger)
    end;

end;


text:
struct DensifyInnerMostBlockLevelMap <: Map
end

function return_cache(k::DensifyInnerMostBlockLevelMap,
                      a::VectorBlock{<:Vector{T}}) where {T}
  @check all(a.touched)
  # Precompute the size of each row block.
  # We are assuming that the size of each row
  # block is going to be equivalent for all cells
  s=0
  brs=Vector{Int}(undef,size(a)[1])
  for i=1:size(a)[1]
    s=s+length(a.array[i])
    brs[i]=length(a.array[i])
  end
  CachedArray(Vector{T}(undef,(s,))),brs
end

function return_cache(k   :: DensifyInnerMostBlockLevelMap,
                      brs :: Vector{<:Integer},
                      a   :: VectorBlock{<:Vector{T}}) where {T}
  @check length(brs)==size(a)[1]
  CachedArray(Vector{T}(undef,(sum(brs),)))
end

function return_cache(k::DensifyInnerMostBlockLevelMap,
                      a::MatrixBlock{<:Vector{T}}) where {T}
  @check _check_preconditions(k,a)
  s=Vector{Int}(undef,2)
  s[1]=0
  s[2]=size(a)[2]
  # Precompute the size of each row block.
  # We are assuming that the size of each row
  # block is going to be equivalent for all cells
  brs=Vector{Int}(undef,size(a)[1])
  for i=1:size(a)[1]
    for j=1:size(a)[2]
      if (a.touched[i,j])
        s[1]=s[1]+length(a.array[i,j])
        brs[i]=length(a.array[i,j])
        break
      end
    end
  end
  CachedArray(Array{T,2}(undef,Tuple(s))),brs
end

function return_cache(k   :: DensifyInnerMostBlockLevelMap,
                      brs :: Vector{<:Integer},
                      a   :: MatrixBlock{<:Vector{T}}) where {T}
  @check length(brs)==size(a)[1]
  CachedArray(Array{T,2}(undef,(sum(brs),size(a)[2])))
end

function _check_preconditions(k::DensifyInnerMostBlockLevelMap,
                              a::MatrixBlock{<:Vector{T}}) where {T}
  # Double check that, for each row, there is at least one non-zero block
  found=false
  for i=1:size(a)[1]
    found=false
    for j=1:size(a)[2]
      if (a.touched[i,j])
        found = true
        break
      end
    end
    if !found
      break
    end
  end
  found
end

function return_cache(k::DensifyInnerMostBlockLevelMap,
                      a::VectorBlock{<:Matrix{T}}) where {T}
  @check all(a.touched)
  s=Vector{Int}(undef,2)
  s[1]=0
  s[2]=size(a.array[1])[2]
  for i=1:size(a)[1]
    s[1]=s[1]+size(a.array[i])[1]
  end
  CachedArray(Array{T,2}(undef,Tuple(s)))
end

function return_cache(k  ::DensifyInnerMostBlockLevelMap,
                      brs::Vector{<:Integer},
                      cs ::Integer,
                      a  ::VectorBlock{<:Matrix{T}}) where {T}
  CachedArray(Array{T,2}(undef,(sum(brs),cs)))
end

function return_cache(k::DensifyInnerMostBlockLevelMap,
                      a::MatrixBlock{<:Matrix{T}}) where {T}
  @check _check_preconditions(k,a)
  s=Vector{Int}(undef,2)
  s[1]=0
  s[2]=0
  # Precompute the size of each row/col block.
  # We are assuming that the size of each row/col
  # block is going to be equivalent for all cells
  brs=Vector{Int}(undef,size(a)[1])
  bcs=Vector{Int}(undef,size(a)[2])
  # Row traversal
  for i=1:size(a)[1]
    for j=1:size(a)[2]
      if (a.touched[i,j])
        s[1]=s[1]+size(a.array[i,j])[1]
        brs[i]=size(a.array[i,j])[1]
        break
      end
    end
  end
  #Column traversal
  for j=1:size(a)[2]
     for i=1:size(a)[1]
       if (a.touched[i,j])
        s[2]=s[2]+size(a.array[i,j])[2]
        bcs[j]=size(a.array[i,j])[2]
        break
       end
     end
  end
  CachedArray(Array{T,2}(undef,Tuple(s))),brs,bcs
end

function _check_preconditions(k::DensifyInnerMostBlockLevelMap,
                              a::MatrixBlock{<:Matrix{T}}) where {T}
  # Double check that, for each row,
  # there is at least one non-zero block
  found=false
  for i=1:size(a)[1]
    found=false
    for j=1:size(a)[2]
      if (a.touched[i,j])
        found = true
        break
      end
    end
    if !found
      break
    end
  end
  # Double check that, for each col,
  # there is at least one non-zero block
  if (found)
    for j=1:size(a)[2]
      found=false
      for i=1:size(a)[1]
        if (a.touched[i,j])
          found = true
          break
        end
      end
      if !found
        break
      end
    end
  end
  found
end

function return_cache(k::DensifyInnerMostBlockLevelMap,
                      brs::Vector{<:Integer},
                      bcs::Vector{<:Integer},
                      a::MatrixBlock{<:Matrix{T}}) where {T}
  CachedArray(Array{T,2}(undef,(sum(brs),sum(bcs))))
end

function evaluate!(cache,
                   k::DensifyInnerMostBlockLevelMap,
                   a::VectorBlock{<:Vector{T}}) where {T}
  cache_array, brs = cache
  _evaluate!(cache_array,brs,a)
end

function evaluate!(cache,
                   k   :: DensifyInnerMostBlockLevelMap,
                   brs :: Vector{<:Integer},
                   a   :: VectorBlock{<:Vector{T}}) where {T}
  _evaluate!(cache,brs,a)
end

function _evaluate!(cache,
                    brs::Vector{<:Integer},
                    a::VectorBlock{<:Vector{T}}) where {T}
  output  = cache.array
  output .= zero(eltype(output))
  current_i=1
  for i=1:size(a)[1]
    range = current_i:current_i+brs[i]-1
    if (a.touched[i])
      output[range] = a.array[i]
    end
    current_i = current_i + length(range)
  end
  output
end

function evaluate!(cache,
                   k::DensifyInnerMostBlockLevelMap,
                   a::MatrixBlock{<:Vector{T}}) where {T}
  cache_array, brs = cache
  _evaluate!(cache_array,brs,a)
end

function evaluate!(cache,
                   k   :: DensifyInnerMostBlockLevelMap,
                   brs :: Vector{<:Integer},
                   a   :: MatrixBlock{<:Vector{T}}) where {T}
  _evaluate!(cache,brs,a)
end

function _evaluate!(cache,
                    brs::Vector{<:Integer},
                    a::MatrixBlock{<:Vector{T}}) where {T}
  output  = cache.array
  output .= zero(eltype(output))
  for j=1:size(a)[2]
     current_i=1
     for i=1:size(a)[1]
        range = current_i:current_i+brs[i]-1
        if (a.touched[i,j])
          output[range,j] = a.array[i,j]
        end
        current_i = current_i + length(range)
     end
  end
  output
end

function evaluate!(cache,
                   k::DensifyInnerMostBlockLevelMap,
                   a::MatrixBlock{<:Matrix{T}}) where {T}
  cache_array, brs, bcs = cache
  _evaluate!(cache_array, brs, bcs,a)
end

function evaluate!(cache,
  k::DensifyInnerMostBlockLevelMap,
  brs::Vector{<:Integer},
  bcs::Vector{<:Integer},
  a::MatrixBlock{<:Matrix{T}}) where {T}
  _evaluate!(cache, brs, bcs, a)
end

function _evaluate!(cache_array, brs, bcs, a)
  output  = cache_array.array
  output .= zero(eltype(output))
  current_j=1
  for j=1:size(a)[2]
    current_i=1
    range_j = current_j:current_j+bcs[j]-1
    for i=1:size(a)[1]
      range_i = current_i:current_i+brs[i]-1
      if (a.touched[i,j])
         output[range_i,range_j] = a.array[i,j]
      end
      current_i = current_i + brs[i]
    end
    current_j = current_j + bcs[j]
  end
  output
end

function evaluate!(cache,
                   k::DensifyInnerMostBlockLevelMap,
                   a::VectorBlock{<:Matrix{T}}) where {T}
  @check all(a.touched)
  output = cache.array
  current_i=1
  n=size(a.array[1])[2]
  for i=1:size(a)[1]
    range=current_i:current_i+size(a.array[i])[1]-1
    output[range,1:n] = a.array[i]
    current_i=current_i+length(range)
  end
  output
end

function return_cache(k::DensifyInnerMostBlockLevelMap,
      a::ArrayBlock{<:ArrayBlock{T,M} where {T,M},N}) where {N}
    cache_touched=a.touched
    i=findfirst(isone, cache_touched)
    cache_block=return_cache(k,a.array[i])
    cache_array=Array{typeof(cache_block),N}(undef,size(a))
    output_array=Array{return_type(k,a.array[i]),N}(undef,size(a))
    linds=LinearIndices(size(a))
    cinds=CartesianIndices(size(a))
    while i != nothing
      cache_array[i]=cache_block
      if (linds[i]+1 <= length(linds))
        i=findnext(isone, cache_touched, cinds[linds[i]+1])
        if (i!=nothing)
          cache_block=return_cache(k,a.array[i])
        end
      else
        i=nothing
      end
    end
    ArrayBlock(cache_array,cache_touched),
    ArrayBlock(output_array,cache_touched),
    linds,
    cinds
end

function evaluate!(cache,
    k::DensifyInnerMostBlockLevelMap,
    a::ArrayBlock{<:ArrayBlock{T,M} where {T,M},N}) where {N}
    cache_array, output_array, linds, cinds = cache
    @check cache_array.touched == a.touched
    i=findfirst(isone, cache_array.touched)
    while i != nothing
      output_array.array[i]=evaluate!(cache_array.array[i],k,a.array[i])
      if (linds[i]+1 <= length(linds))
        i=findnext(isone, cache_array.touched, cinds[linds[i]+1])
      else
        i=nothing
      end
    end
    output_array
end


text:
println("Building...")

builddir = mktempdir()
# println("$builddir / $")
# cmd = `cmake `


text:
using MAT

#useproc("CPU")
useproc("GPU")

Ntrain=60000
Ntest=10000
#BatchSize=Ntrain
BatchSize=100
TrainingIts=1000 # number of Nesterov updates
include("loadmnist.jl")
images,label=loadmnist()
label=convert(Array{Int,2},label)
r=randperm(size(images,2))
data=images[:,r]; label=label[r]
class=zeros(10,70000)+0.00001
for i=1:70000
    class[label[i]+1,i]=0.9999
end
traindata=data[:,1:Ntrain]
testdata=data[:,Ntrain+1:Ntrain+Ntest]
trainclass=class[:,1:Ntrain]
testclass=class[:,Ntrain+1:Ntrain+Ntest]

# Construct the DAG function:
H=[784 10] # number of units in each layer
L=length(H) # number of hidden layers
# node indices:
w=Array(ADnode,L) # weight index
h=Array(ADnode,L) # hidden layer index (note that I call the input layer h[1])

StartCode()
x=h[1]=ADnode()
class=ADnode() # node has no parents
for i=2:L-1
    w[i]=ADvariable()
    #h[i]=sigmoid(w[i]*h[i-1])
#   h[i]=abs(w[i]*h[i-1])
#    h[i]=abs(w[i]*h[i-1])+1.5*w[i]*h[i-1]
    h[i]=rectlin(w[i]*h[i-1])
end
w[L]=ADvariable()
h[L]=w[L]*h[L-1]
loss=KLsoftmax(class,h[L]) # the loss we are minimising must be the final node in the graph.
net=EndCode()

# instantiate parameter nodes and inputs:
net.value[x]=traindata[:,1:BatchSize]
net.value[class]=trainclass[:,1:BatchSize]
for i=2:L
    net.value[w[i]]=sign(randn(H[i],H[i-1]))/sqrt(H[i-1])
end

net=compile(net) # compile the DAG and preallocate memory
gradcheck(net)
@gpu CUDArt.init([0]) # let the user do device management
@gpu net=convert(net,"GPU")

println("Training: using $(net.gpu==true? "GPU" : "CPU") ")
starttime=time()
error=Array(Float64,0)
ParsToUpdate=Parameters(net)
velo=NesterovInit(net)
minibatchstart=1 # starting datapoint for the minibatch
for iter=1:TrainingIts
    LearningRate=0.1/(1+iter/500)
    minibatchstart,minibatch=GetBatch(minibatchstart,BatchSize,Ntrain)
    net.value[x]=cArray(PROC,traindata[:,minibatch]) # select batch of data
    net.value[class]=cArray(PROC,trainclass[:,minibatch]) # select batch of data
    ADforward!(net)
    ADbackward!(net)
    push!(error,extract(net.value[net.FunctionNode])[1])
    printover("iteration $iter: training loss = $(error[iter])")
    for par in ParsToUpdate
            #GradientDescentUpdate!(net.value[par],net.gradient[par],LearningRate)
        NesterovGradientDescentUpdate!(net.value[par],net.gradient[par],velo[par],LearningRate,iter)
    end
end
endtime=time()
println("\nTraining took $(endtime-starttime) seconds")

net=convert(net,"CPU")

classpredtrain=argcolmax(softmax(net.value[h[L]]))
classtrain=argcolmax(net.value[class])
println("last batch train accuracy = $(mean(classpredtrain.==classtrain))")

println("Training:")
# evaluate the test predictions:
net.value[x]=traindata
ForwardPassList!(net,ExcludeNodes=[class])
# ADforward! by default computes values inplace. If we change the input value dimensions (which happens here since there are a different number of test points to train points), we need to do an out-of-place forward pass in order to compute the shapes of the values on the graph. This only needs to be done once.
ADforward!(net,AllocateMemory=true)

net=convert(net,"CPU") # to make analysis simpler

classpredtrain=argcolmax(softmax(net.value[h[L]]))
classtrain=argcolmax(trainclass)
println("train accuracy = $(mean(classpredtrain.==classtrain))")

println("Testing:")
# evaluate the test predictions:
net.value[x]=testdata
ForwardPassList!(net,ExcludeNodes=[loss])
ADforward!(net,AllocateMemory=true) # need to reallocate memory since the shape of net.value[x] has changed

net=convert(net,"CPU")
classpredtest=argcolmax(softmax(net.value[h[L]]))
classtest=argcolmax(testclass)
println("test accuracy = $(mean(classpredtest.==classtest))")

#end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


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

using Markdown
using InteractiveUtils

# ╔═╡ 8c8388cf-9891-423c-8db2-d40d870bb38e
using PlutoUI

# ╔═╡ 54778212-3702-467d-8a96-4fa18b3ccd63
using LinearAlgebra

# ╔═╡ 6f4274b5-87e2-420d-83d2-83a8408650cd
# edit the code below to set your name and kerberos ID (i.e. email without @mit.edu)

student = (name = "Jazzy Doe", kerberos_id = "jazz")

# you might need to wait until all other cells in this notebook have completed running.
# scroll around the page to see what's up

# ╔═╡ ea28bf57-ba62-41ce-8be6-d38ca2c5caa3
[j for j=1:1000 if j*(1/j) != 1]

# ╔═╡ f0e59676-fd8e-4609-833e-f7d1a5a6e9e3
smallest_j = 49

# ╔═╡ a5be0775-68de-41ce-95cd-1465723d099b
error = 1 - smallest_j * (1 / smallest_j)

# ╔═╡ 4abec609-09cd-4f86-8d86-a7d02325cc7b
log2(error)

# ╔═╡ d221c61c-a4ab-4a82-b89d-52735d957cae
32 * 23 - 736

# ╔═╡ c3d49b98-a9c1-4aba-becc-7fa84f4bbc75
3.2 * 23 - 73.6

# ╔═╡ 7b33c09c-2ef2-4b97-b5a5-5fdf9268f76b
3.2 * 2.3 - 7.36

# ╔═╡ 615b88b6-8505-11eb-0a7c-294ee5ae7474
begin
	struct FirstRankOneMatrix
		# Your code here
		v::Vector{Float64}
		w::Vector{Float64}
		FirstRankOneMatrix(v, w) = new(v, w)
	end

	# Add the extra constructor here
	FirstRankOneMatrix(v) = FirstRankOneMatrix(v, v)

end

# ╔═╡ b0db7388-850c-11eb-0915-597f1fa5ab93
ten_twelve = FirstRankOneMatrix(1:10, 1:12) # Your code here

# ╔═╡ 2dfc2ef9-3efa-4b82-960f-f2ef0171e9eb
sqrt === Base.sqrt

# ╔═╡ d4a0baea-cc02-4d91-a938-8181dea7b47a
filter === Base.filter

# ╔═╡ fada4734-8505-11eb-3f2b-d1f1ef391ba4
function Base.size(M::FirstRankOneMatrix)
	return length(M.v), length(M.w)
end

# ╔═╡ 590dfe1a-8506-11eb-0069-d7cd91f02a65
function Base.getindex(M::FirstRankOneMatrix, i, j)
	return M.v[i] * M.w[j]
end

# ╔═╡ 941b6b10-8ae9-11eb-3bf8-b732f5f60af3
ten_twelve

# ╔═╡ a243a400-8af3-11eb-0637-cf8f80aae86d
ones(10, 12) # An example matrix (two-dimensional array)

# ╔═╡ aa756549-1b9f-41a1-a2e1-713166e4211a
(ones(12, 12) + ones(12, 12)) * (ones(12, 12) + ones(12, 12))

# ╔═╡ b577420c-8501-11eb-267a-719125315fe1
begin
	struct RankOneMatrix{T} <: AbstractMatrix{T}
		v::AbstractVector{T}
		w::AbstractVector{T}
		RankOneMatrix(v::AbstractVector{T}, w::AbstractVector{T}) where T = new{T}(v, w)
	end

	# Add the two extra constructors
	# (Should we make these missing by default? if so - remove hint below)
	RankOneMatrix(v) = RankOneMatrix(v, v)
end

# ╔═╡ eb612772-8b06-44bb-a36a-827435cbb2ee
RankOneMatrix([1,2], [3,4,5]);

# ╔═╡ d3e7c8d1-4b4a-47b6-9c96-150333078f42
RankOneMatrix(1:10);

# ╔═╡ f2d8b45c-8501-11eb-1c6a-5f819c240d9d
function Base.size(M::RankOneMatrix)
	return length(M.v), length(M.w)
end

# ╔═╡ ed72e880-8afa-11eb-3a4a-175a838188d9
function Base.getindex(M::RankOneMatrix, i, j)
	return M.v[i] * M.w[j]
end

# ╔═╡ 7b3fb0ef-9a9e-401c-8c09-e5615134a4ad
R2 = RankOneMatrix([1,2], [3,4,5])

# ╔═╡ fc962c72-8501-11eb-2821-cbb7a52d5f61
M = RankOneMatrix(1:10)

# ╔═╡ 0887ee78-0e8a-41c7-90e7-44237acc1477
collect(M)

# ╔═╡ 1705dc7b-9f93-440b-acbb-ecb362d08125
typeof(collect(M))

# ╔═╡ ba3cb45a-8502-11eb-2141-6369b0e08807
begin
	struct RankTwoMatrix{T} <: AbstractMatrix{T}
		# Your code here
		A::RankOneMatrix{T}
		B::RankOneMatrix{T}
		RankTwoMatrix(A::RankOneMatrix{T}, B::RankOneMatrix{T}) where T = new{T}(A, B)
	end

	# Add a constructor that uses two vectors/ranges
	RankTwoMatrix(v1, v2) = RankTwoMatrix(RankOneMatrix(v1), RankOneMatrix(v2))
end

# ╔═╡ f5a95dd8-850d-11eb-2aa7-2dcb1868577f
RankTwoMatrix(1.0:10.0, 0.0:0.1:0.9);

# ╔═╡ c784e02c-8502-11eb-3efa-7f4c45f4274c
function Base.getindex(M::RankTwoMatrix, i, j)
	return M.A[i, j] + M.B[i, j]
end

# ╔═╡ 0bab818e-8503-11eb-02b3-178098599847
function Base.size(M::RankTwoMatrix)
	return size(M.A)
end

# ╔═╡ b6717812-8503-11eb-2729-39bfdc1fd2f9
struct LowRankMatrix <: AbstractMatrix{Float64}
	Ms::Vector{RankOneMatrix}
	rank::Int
	LowRankMatrix(Ms, rank) = new(Ms, rank)
end

# ╔═╡ c49e350e-8503-11eb-15de-7308dd03dc08
function Base.getindex(M::LowRankMatrix, i, j)
	return sum(M.Ms[k][i, j] for k=1:M.rank)
end

# ╔═╡ eadb174e-2c1d-48c8-9de2-99cdc2b38d32
md"_homework 5, version 2_"

# ╔═╡ 0560cf7b-9986-402a-9c40-779ea7a7292d
md"""

Submission by: **_$(student.name)_** ($(student.kerberos_id)@mit.edu)
"""

# ╔═╡ 082542fe-f89d-4774-be20-1e1a78f21291
md"""

# **Homework 5**: _Structure_
`18.S191`, Spring 2021

This notebook contains _built-in, live answer checks_! In some exercises you will see a coloured box, which runs a test case on your code, and provides feedback based on the result. Simply edit the code, run it, and the check runs again.

_For MIT students:_ there will also be some additional (secret) test cases that will be run as part of the grading process, and we will look at your notebook and write comments.

Feel free to ask questions!
"""

# ╔═╡ aaa41509-a62d-417b-bca7-a120e3a5e5b2
md"""
#### Intializing packages
_When running this notebook for the first time, this could take up to 15 minutes. Hang in there!_
"""

# ╔═╡ ac35c6b3-a142-4b7a-921b-a2bfe1d84713
md"""
## **Exercise 1:** _When is zero not quite zero?_
Students used to pure math are sometimes surprised to see numbers like `1e-15`, `1e-16`, or `1e-17` appearing in computations when `0` might have been expected. At first glance, this behaviour seems 'random' or 'noisy', but in this exercise, we will try to demonstrate some structure behind floating point artihmetic.

#### Exercise 1.1
👉 Find all integers $j$ with $1 ≤ j ≤ 1000$ for which Julia's result satisfies `j*(1/j) ≠ 1`.

"""

# ╔═╡ a5901f93-007f-4a30-97fc-29b367ec47c6
md"""
Notice that when you re-run the computation, the result does not change. Floating-point arithmetic is not random!

#### Exercise 1.2

👉 Take the smallest number `j` you found above and compute the error, i.e. the distance between `j*(1/j)` and `1`.

Is this an integer power of 2?  Which one? (`log2` might help.)
"""

# ╔═╡ 32e073bb-943f-4fa9-b15f-5ec18feecf15
md"""
The error is indeed a power of 2, it's 1/2^53
"""

# ╔═╡ 7d1e1724-cd1a-483b-af40-f24ae5301849
md"""
#### Exercise 1.3
Caculate all of the following:
* 32 * 23  - 736
* 3.2 * 23 - 73.6
* 3.2 * 2.3 - 7.36
"""

# ╔═╡ 224acb51-c0db-4f04-b865-8349d6a28e98
md"""
We wanted to show you that floating point arithmetic has structure to it.  It is not the fuzz or slop that you may have seen as experimental errors in maybe a chemistry or physics class.  If you are interested, we may show more later in the semester.
"""

# ╔═╡ a7e9ab20-8501-11eb-2a7a-af8fd495b21e
md"""
## **Exercise 2:** _Rank-one matrices_
In this exercise we will go into some more detail about how to define types in Julia in order to make a structured-matrix type, similar to the `OneHotVector` type from class.

To begin, we will make a type to represent a rank-1 matrix. A *rank-1 matrix* is a matrix whose columns are all multiples of each other. This already tells us a lot about the matrix; in fact, we can represent any such matrix as the outer product of two vectors `v` and `w`. *Only* the two vectors will be stored; the matrix elements will be *calculated* on demand, i.e. when we *index* into the object.
"""

# ╔═╡ 8ce8c4bc-8505-11eb-2357-c50a70b8745c
md"""
#### Exercise 2.1
Let's make a `FirstRankOneMatrix` type that contains two vectors of floats, `v` and `w`. Here `v` represents a column and `w` the multpliers for each column.

We include (in the same cell, due to requirements of Pluto) a constructor that takes a single vector `v` and duplicates it.
"""

# ╔═╡ 9f534aa0-8505-11eb-3476-4d326fb19d42
md"""
👉 Create an object of type `FirstRankOneMatrix` representing the multiplication table of the numbers from 1 to 10 and the numbers from 1 to 12. Call it `ten_twelve`.
"""

# ╔═╡ 7aac43b3-49ce-46ea-854b-d292ac0591c7
md"""
#### Exercise 2.2 - _extending Base methods_
Right now, our `FirstRankOneMatrix` is just a container for two arrays. To make it _behave_ like a matrix, we add methods to functions in `Base`, Julia's standard set of functions. Most of the functions you have used so far come from base, for example:
"""

# ╔═╡ ce6e648c-ade5-4882-9f92-1ae3f878a5a3
md"""
These are built-in functions, and each function comes with a set of methods pre-defined (in Julia's source code _(which is mostly written in Julia)_). Uncomment the next cell to see the full list of methods for `size`. We will add a method to this list!
"""

# ╔═╡ 87fce3d8-765c-410c-a059-4dd51bd97254
md"""

Let's extend the `size` and `getindex` methods so that they can also work with our rank-1 matrices. This will allow us to easily access the size and entries of the matrix.
"""

# ╔═╡ e86d8ab6-8505-11eb-1995-137dfbff0191
md"""
👉 Extend `Base.size` for the type `FirstRankOneMatrix` to return the size as a tuple of values corresponding to each direction.
"""

# ╔═╡ 545a5916-8506-11eb-1a31-078957d791f8
md"""
👉 Extend `Base.getindex` for the type `FirstRankOneMatrix` to return the $(i, j)$th entry of the outer product.
"""

# ╔═╡ fe6df9bf-6059-4b76-af39-385d395ece72
Base.getindex(ten_twelve, 5, 11)

# ╔═╡ b9cb6192-8505-11eb-3d2c-790654bbc9a1
md"""
#### Exercise 2.3

If you ask Julia to display the value of an object of type `FirstRankOneMatrix`, you will notice that it doesn't do what we are used to, which is probably to display the whole matrix that it corresponds to. Let's see:
"""

# ╔═╡ 1b3cc404-9da8-4acd-9df5-5815ce167c34
md"""
This is what matrices normally look like in Julia:
"""

# ╔═╡ 4e926bfe-8734-11eb-1355-e1d1d9b36929
md"""
One possible workaround for this is to create a new function specifically to display a matrix of our custom type.

👉 Make a function `print_as_matrix` that prints the entries of a matrix of type `FirstRankOneMatrix` in a grid format. We'll test that it works below with the matrix `R` we already created.
"""

# ╔═╡ 9eb846c0-8737-11eb-101b-0191f715d8c9
md"""
We should now be able to see our `FirstRankOneMatrix`-type matrix displayed in the terminal! In Pluto, we have to use `with_terminal` to capture and display output printed to the _terminal_.
"""

# ╔═╡ a75f3c76-8506-11eb-1d83-cf5781b656a3
md"""
#### Exercise 2.4 - _AbstractMatrix_
In fact, Julia (together with Pluto) can do some of this work for us! Julia provides facilities to make our life easier when we tell Julia that our type *behaves like* an array. We do so by making the type a *subtype of* `AbstractMatrix`. This will let our type **inherit** the methods and attributes that come with the `AbstractMatrix` type by default - including some display functions that Pluto notebooks use.

Let's do all this on a new rank-one matrix type: `RankOneMatrix`. As we do this, we will also remove the restriction to `Float64`-type entries by using a *parametrised type* `T`. Parametrised types allow us flexibility to handle different types of entries without repeating a lot of code; we won't go too in-depth about these for now.
"""

# ╔═╡ 1d16e228-850e-11eb-1dca-41df85b706da
md"""
In the cell above, we added a second ('outer') constructor that takes a single vector as argument.

👉 Make sure that you can use both constructors by trying them out below.
"""

# ╔═╡ 3ed18ba0-850b-11eb-3e90-f188c54d9ce9
md"""
👉 Add `getindex` and `size` methods for this new type. These will allow us to access entries of our custom matrices with the usual index notation `M[i,j]`, as well as quickly retrieving their dimensions.
"""

# ╔═╡ 5a6dd3bc-af85-4e7d-b03c-b8c1bc88d850
Base.getindex(R2, 2, 3)

# ╔═╡ 91b1125a-850b-11eb-1a21-c7af293332fd
md"""
👉 Create an object of the new type for the $10 \times 10$ multiplication matrix, using a single range.
"""

# ╔═╡ 8f5c7c9c-850b-11eb-3d97-bf9c6b5d7d2e
md"""
You should see two things: Firstly, the matrix now contains integers, instead of floats (this is thanks to our parametrised type!). And secondly, Julia *automatically* displays the matrix as we wanted, once you have defined `getindex` and `size`, provided you have told Julia that your type is a subtype of `AbstractArray`!
"""

# ╔═╡ 248ab23e-850c-11eb-396f-5ba365f0f743
md"""
Julia also allows us to automatically convert the matrix to a normal ("dense") matrix, using either `collect` or `Array`. Let's try these out. You may need to re-run the cell below after completing the exercises.
"""

# ╔═╡ 8cab8d7d-32fc-4b43-bd54-fb0c78fa5d00
md"""
> Fun fact: we did not define a method for `Base.collect`, and yet it works! This is because we told Julia that our `RankOneMatrix` is a subtype of `AbstractMatrix`, for which `Base.collect` already has a method defined in Julia's source code. This fallback method uses `Base.getindex` and `Base.size`, which we did define.
"""

# ╔═╡ 6a08cb74-8510-11eb-2212-2f1f676f3d5d
md"""
#### Exercise 2.5

Why do we need a special type to represent special types of structured matrices? One reason is that not only do they give a more efficient representation in space (requiring less memory to store), they can also be more efficient in time, i.e. *faster*.

   For example, let's look at **matrix--vector mutiplication**. This is a *fundamental* part of many, *many* algorithms in scientific computing, and because of this, we usually want it to be as fast as possible.

   For a rank-one matrix given by $M = v w^T$, the matrix--vector product $M \cdot x$ is given by $(w \cdot x) v$. Note that $w \cdot x$ is a number (scalar) which is multiplying the vector element by element. This computation is much faster than the usual matrix-vector multiplication: we are taking advantage of structure!

"""

# ╔═╡ 61cbfc80-86f5-11eb-316f-6307cac7cbbb
md"""
👉 Define a function `matvec` that takes a `RankOneMatrix` `M` and a `Vector` `x` and carries out matrix-vector multiplication. We will be able to compare the result with doing the matrix--vector product using the corresponding dense matrix.
"""

# ╔═╡ 4864e4d6-850c-11eb-210c-0318b8660a9a
md"""
## Exercise 3: Low-rank matrices
"""

# ╔═╡ 58ce8d2a-850c-11eb-29e9-3fc966e129b1
md"""
In this exercise we will combine rank-1 matrices into low-rank matrices, which are sums of rank-1 matrices. Just like in the previous exercise, we will make custom types for these matrices, and we will be able to compute the entries of a matrix on-demand as we request them (instead of storing all of them always).
"""

# ╔═╡ 9d4d943e-8502-11eb-31cb-e11824b2dce0
md"""
#### Exercise 3.1

👉 Make a rank-2 matrix type `RankTwoMatrix` that contains two rank-1 matrices `A` and `B`.  For simplicity, you can use the `FirstRankOneMatrix` type.

   Include in the same cell a constructor for `RankTwoMatrix` that takes two vectors and makes the rank-one matrices from those vectors.

*Note*: In principle we should check when constructing the type that the input matrices have the same dimensions, but we will just assume that they do.
"""

# ╔═╡ 92a91904-850c-11eb-010e-c58ae218f541
md"""
#### Exercise 3.2
👉 Construct a rank two matrix out of the rank-1 matrix representing the multiplication table of $1.0:10.0$, together with the multiplication table of $0.0:0.1:0.9$.

   (Note that both range arguments must contain floats, so that we can add up entries.)
"""

# ╔═╡ dc714540-8afa-11eb-205b-7770074771c8
md"""
#### Exercise 3.3
👉 As with last time, extend the `getindex` and `size` methods for the `RankTwoMatrix` type. Keep in mind they are already defined for `RankOneMatrix`.
"""

# ╔═╡ aca709f0-8503-11eb-1144-1fc01cb85c39
md"""
#### Exercise 3.4

Making a custom type for rank-2 matrices is a step forward from rank-1 matrices: instead of storing two vectors, we store two rank-1 matrices themselves. What if we want to represent a rank-3 matrix? We would need to store *three* rank-1 matrices, instead of just two. What about rank-4, rank-5, and so on?

We can go even further and make a general custom type `LowRankMatrix` for rank-$k$ matrices, for general (ideally low) $k$. In this case, we should store two main things: the list of rank-1 matrices that our low-rank matrix is made up of, and also the *rank* of the matrix (which is how many rank-1 matrices we are storing).

👉 Complete the definition for the type `LowRankMatrix`. Remember to store both the rank-1 matrices and the rank of the matrix itself.
"""

# ╔═╡ 6e280730-8b00-11eb-04b7-615a389e792f
md"""
#### Exercise 3.5
👉 Extend the `getindex` and `size` methods to work with the `LowRankMatrix` type. As before, remember that these are already defined for `RankOneMatrix`.
"""

# ╔═╡ dd27f508-8503-11eb-36b9-33f5f99f78b0
function Base.size(M::LowRankMatrix)
	return size(M.Ms[1])
end

# ╔═╡ ead4c008-0e7e-4414-aced-d4a576423bd3
size === Base.size

# ╔═╡ 7adbcbab-66e7-4b58-b197-5c220ba1e9a2
methods(size);

# ╔═╡ 0b7c6cbe-57de-419d-adcb-8724791f9c89
Base.size(ten_twelve)

# ╔═╡ c7a15c5e-8505-11eb-3af2-2fa84b74b590
function print_as_matrix(M::FirstRankOneMatrix)
	m, n = size(M)
	println(m, "x", n, " FirstRankOneMatrix")

	pad = length(repr(maximum(M.w) * maximum(M.v))) + 1
	for i=1:m
		for j=1:n
			print(lpad(M[i, j], pad))
		end
		println()
	end
end

# ╔═╡ b0a0b9a4-850b-11eb-30f1-11f5270efe02
with_terminal() do
	print_as_matrix(ten_twelve)
end

# ╔═╡ 52451232-2c3d-4425-a94a-da1a955cf784
Base.size(R2)

# ╔═╡ ee58251a-8511-11eb-074c-5b1e27c4ebd4
function matvec(M::RankOneMatrix, x)
	m, _ = size(M)
	return sum(M.w[i] * x[i] for i in 1:m) * M.v
end

# ╔═╡ 93a75ea0-8b0c-11eb-3946-3d5d92487fb5
let
	comp1, comp2 = RankOneMatrix(1.0:1.0:10.0), RankOneMatrix(0.0:0.1:0.9)
	ex3 = LowRankMatrix([comp1, comp2], 2)
end

# ╔═╡ 1083d5ee-8512-11eb-0fdc-ab3987f94bc5
md"""
#### Exercise 3.6
👉 Extend the method `matvec` to work with our custom type `LowRankMatrix` (and a `Vector`). Remember that `matvec` is already defined for `RankOneMatrix`.
"""

# ╔═╡ 3fed837a-8512-11eb-1fdd-c1b72b48d07b
function matvec(M::LowRankMatrix, x)
	return sum(matvec(M.Ms[k], x) for k=1:M.rank)
end

# ╔═╡ 2d65bd1a-8512-11eb-1bd2-0313588dfa0e
md"""
#### Exercise 3.7
One of the big advantages of our rank-1 matrices is its space efficiency: to "store" a $n\times n$ matrix, we only need to *actually* store two vectors of $n$ entries, as opposed of $n^2$ total entries. Rank-2 matrices are a litte less efficient: we store two rank-1 matrices, or four vectors of $n$ vertices. As we increase the rank, we lose space efficiency, at the cost of being able to represent more kinds of matrices.

👉 For what rank will a matrix of rank $k$ need the same amount of storage as the dense version? Explain your answer.
"""

# ╔═╡ f9556098-8504-11eb-08a0-39fbe00892da
answer = md"""

From rank $k = n/2$ onwards.

Because for a rank k matrix, we store k rank-1 matrices, or $2\times k$ vectors of $n$ vertices. This corresponds to $2\times n\times k$ entries. If we solve the equation $n^2$ = 2*n*k, we get $k = n/2$.
"""

# ╔═╡ 295acdac-880a-402e-9f7e-19b0fc801130
md"""
## **Exercise 4:** _The SVD for structured matrices_
In a math class, you may or may not learn about the *singular value decomposition (SVD)*. From a computational thinking point of view, whether you have seen this before or not we hope will not matter.  Here we would like you to _get to "know" the SVD_, through experience, rather than through a math lecture.  This is how we see "computational thinking" after all. (Definitely do not look up some eigenvalue definition.)
"""

# ╔═╡ 3ba3ae50-b9ed-416e-bb47-a14f82eb4e19
md"""
#### Exercise 4.1

The `LinearAlgebra` package defines a function `svd` that computes the decomposition for us. Have a look at the result of calling the `svd` function. Try not to get intimidated by the large output (you may need to scroll the cell output), and look at the docs for `svd`.
"""

# ╔═╡ bdf3954e-2a05-4a0d-899b-1f6f59119550
biggie = rand(100,100)

# ╔═╡ a58c1c54-4f17-4ca4-ac78-22773b526dfb
md"""
👉 What are the _singular values_ of `biggie`?
"""

# ╔═╡ bb649c89-709c-49c8-8111-53044e8e682a
md"""

#### Exercise 4.2

If we try to run the `svd` function on a `RankOneMatrix`, you will see that it does not work. The error is telling us that **no `svd` method has been defined for our type**. Let's extend the `svd` function to work on objects of our type.
"""

# ╔═╡ 6f78b3bb-5440-4f34-849a-fa8e21972fa8
A = RankOneMatrix(rand(3),rand(4))

# ╔═╡ 44b8fbf1-411f-4ce4-b5d7-e12f2d6be51b
md"""
👉 Extend `svd` to work directly on a rank-one matrix, by writing a new method.

Keep things simple. Inside your method, call `LinearAlgebra.svd` on a type that it is _already defined for_.

"""

# ╔═╡ 6c9ae344-084e-459c-841c-8377451507fd
function LinearAlgebra.svd(A::RankOneMatrix)
	return LinearAlgebra.svd(Matrix(A))
end

# ╔═╡ 167b9580-fa18-4248-8643-a0fde723ecc4
svd(biggie)

# ╔═╡ 210392ff-0a22-4e55-be08-9f58804282cf
singular_values_of_biggie = svd(biggie).S

# ╔═╡ bbc4a595-a63e-4a12-9690-f39a23b30ae1
svd(A)

# ╔═╡ 9ae28bc5-d9cb-479d-8c22-7c9248cd5fa0
LinearAlgebra.svd(A)

# ╔═╡ 0cf29d7e-bd83-4c5f-bd41-e8e7c2587188
md"""
#### Exercise 4.3
👉 Look at the singular values, how many of them are approximately non-zero? Does that make sense?
"""

# ╔═╡ ed4b7ea5-3266-466d-a817-7ab2cc82ac9c
count(x -> !isapprox(x, 0, atol=1e-5), LinearAlgebra.svd(A).S)

# ╔═╡ 5c6aee48-017b-49ff-af1a-a31af30a45a9
md"""

This number (the number of singular values that are positive) is something you will (or have) learn(ed) in a linear algebra class: it is known as the *rank* of a matrix, and is usually defined through a complicated elimination procedure.

👉 Write a function `numerical_rank` that calculates the singular values of a matrix `A`, and returns how many of them are approximately zero.

To keep things simple, you can assume that "approximately zero" means: less than `tol=1e5`.
"""

# ╔═╡ a2b7f0c3-488b-4ce9-aed6-b8ccddac6a57
function numerical_rank(A::AbstractMatrix; tol=1e-5)
	return count(x -> !isapprox(x, 0, atol=tol), LinearAlgebra.svd(A).S)
end

# ╔═╡ d3420859-d558-47cb-aaf7-d51a5e2d1f6e
numerical_rank(A)

# ╔═╡ 207d2f4c-cbe8-4828-ab3f-100809bb93e9
numerical_rank(rand(3,3))

# ╔═╡ 7fdc110a-a3e3-44e5-a547-b75e03e0d21e
md"""
#### Exercise 4.4
👉 Write a function that takes an argument $k$ and generates the sum of $k$ random rank-one matrices of size mxn, and counts the number of essentially non-zero singular values.

   Hint: you can do this with a `for` loop, but it can also be done with `sum`, e.g.

   `sum(i for i=1:10)`.

"""

# ╔═╡ 01b12200-2e7e-4f19-96b9-5b6d6cb03233
function k_rank_ones(k, m, n)
	return sum(RankOneMatrix(rand(m),rand(n)) for _=1:k)
end

# ╔═╡ e1e1067b-93ba-40df-bd09-7599538e6181
k_rank_ones(1, 3, 3)

# ╔═╡ 00f5bea8-b808-483a-ad70-332f521481f5
md"""
#### Exercise 4.5
👉 What is the answer when $m$ and $n$ are both $≥ k$?  What if one of $m$ or $n$ is $< k$?
"""

# ╔═╡ a6bb92f9-0cb2-4fdd-8b67-f286edbbdcb6
numerical_rank(k_rank_ones(5, 10, 10))

# ╔═╡ 57c19601-122b-414f-bc99-56f98c794e61
numerical_rank(k_rank_ones(5, 3, 10))

# ╔═╡ 9813e9f9-9970-4a8b-b3ca-e6a699e6fda4
numerical_rank(k_rank_ones(5, 10, 4))

# ╔═╡ c8829a12-917a-4ff4-87c9-fe2b25aaa99c
numerical_rank(k_rank_ones(5, 2, 4))

# ╔═╡ 5aabbec1-a079-4936-9cd1-9c25fe5700e6
md"## Function library

Just some helper functions used in the notebook."

# ╔═╡ 42d6b87d-b4c3-4ccb-aceb-d1b020135f47
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]))

# ╔═╡ a759c4d5-ca00-4dda-b39c-1093ba9bef0e
md"""
Try writing this without a `for` loop. Ideas:

☝ Use `filter`.

☝ Comprehensions can have `if` clauses, as in `[j for j=1:1000 if 5 < j < 17]`.
""" |> hint

# ╔═╡ e9a53ca0-86f2-11eb-2970-5d0996d1603e
hint(md"""
The extra constructor should have the same name as the type; it should take in a single argument `v`, and use it for both the column and the multipliers.
""")

# ╔═╡ 5e85ac82-8734-11eb-2f69-cda569326f68
hint(md"""
Check the docs for functions like `print`, `println`, `lpad` and `rpad`. Also, you don't have to know exactly how the `do` keyword works for now, but feel free to read about it - it can be useful with iterating!
""")

# ╔═╡ 7493be50-8738-11eb-37e7-990f1656998a
hint(md"""
If you want to specify that your argument is a range, use the `AbstractRange` type.
""")

# ╔═╡ f7fb3b32-8b00-11eb-0c0e-8f8d68f849a6
hint(md"""
How would we implement `matvec` for a rank-2 matrix - making use of the two stored rank-1 matrices? Can you generalize this to the rank-$k$ case, when we have many rank-1 matrices?
""")

# ╔═╡ 8fd9e34c-8f20-4979-8f8a-e3e2ae7d6c65
almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]))

# ╔═╡ d375c52d-2126-4594-b819-aba9d34fe77f
still_missing(text=md"Replace `missing` with your answer.") = Markdown.MD(Markdown.Admonition("warning", "Here we go!", [text]))

# ╔═╡ cd02acd0-8af4-11eb-3726-d7a360b98433
let
	example = FirstRankOneMatrix([1,2], [3,4])
	if ismissing(print_as_matrix(example))
		still_missing()
	end
end

# ╔═╡ 2d7565d4-88da-4e41-aad6-958ed6b3b2e0
keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]))

# ╔═╡ d7f45d51-f426-4353-af58-858395295ce0
yays = [md"Fantastic!", md"Splendid!", md"Great!", md"Yay ❤", md"Great! 🎉", md"Well done!", md"Keep it up!", md"Good job!", md"Awesome!", md"You got the right answer!", md"Let's move on to the next question."]

# ╔═╡ 9b2aa95d-4583-40ec-893c-9fc751ea22a1
correct(text=rand(yays)) = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]))

# ╔═╡ 29db89d6-ac40-4121-b996-e083555bc5db
try
	A = RankOneMatrix([1,2,3], [5,6,7])
	result = LinearAlgebra.svd(A)

	if result isa Missing
		still_missing()
	elseif !(result isa LinearAlgebra.SVD)
		almost(md"Like all other `svd` methods (implemented in `LinearAlgebra`), your method should return an SVD factorization (of type `LinearAlgebra.SVD`).")
	else
		reconstructed = result.U[:,1:1] * result.S[1:1] * result.Vt[1:1,:]

		if sum(abs.(A .- reconstructed)) < 1e-5
			correct()
		else
			keep_working()
		end
	end
catch e
	if (e isa MethodError) && e.f === LinearAlgebra.svd
		still_missing(md"Define a new method for `LinearAlgebra.svd` that takes a `RankOneMatrix` as argument.")
	else
		rethrow(e)
	end
end

# ╔═╡ 5dca2bb3-24e6-49ae-a6fc-d3912da4f82a
not_defined(variable_name) = Markdown.MD(Markdown.Admonition("danger", "Oopsie!", [md"Make sure that you define a variable called **$(Markdown.Code(string(variable_name)))**"]))

# ╔═╡ 30e5a813-9ea1-48d9-9039-012b8b267256
let
	if !@isdefined(ten_twelve)
		not_defined(:ten_twelve)
	else
		if ten_twelve isa Missing
			still_missing()
		elseif !(ten_twelve isa FirstRankOneMatrix)
			keep_working(md"`ten_twelve` should be a `FirstRankOneMatrix`.")
		elseif (ten_twelve.v == 1:10 && ten_twelve.w == 1:12) ||
				(ten_twelve.w == 1:10 && ten_twelve.v == 1:12)

			correct()
		else
			keep_working()
		end
	end
end

# ╔═╡ 06040694-b324-4e6b-8314-3bffc7b8010f
let
	if !@isdefined(RankOneMatrix)
		not_defined(:RankOneMatrix)
	else
		R2 = FirstRankOneMatrix([1,2], [3,4,5])

		if ismissing(Base.size(R2))
			still_missing()
		elseif size(R2) != (2, 3)
			keep_working(md"`Base.size` is not implemented correctly.")
		else
			correct()
		end
	end
end

# ╔═╡ 4ba8a0f2-850b-11eb-1c09-b7bba1763ed7
let
	if !@isdefined(FirstRankOneMatrix)
		not_defined(:FirstRankOneMatrix)
	else
		R2 = FirstRankOneMatrix([1,2], [3,4,5])

		if ismissing(Base.getindex(R2, 1, 1))
			still_missing(md"Implement both `Base.size` and `Base.getindex`.")
		elseif ismissing(Base.size(R2))
			still_missing(md"Implement both `Base.size` and `Base.getindex`.")
		elseif size(R2) != (2, 3)
			keep_working(md"`Base.size` is not implemented correctly.")
		else
			shouldbe = [
					3 4 5
					6 8 10
				]
			r = all(CartesianIndices(shouldbe)) do I
				isapprox(R2[Tuple(I)...], shouldbe[I], rtol=1-3)
			end
			if r
				correct()
			else
				keep_working(md"`Base.getindex` is not implemented correctly.")
			end
		end
	end
end

# ╔═╡ 111a8db4-a808-4f37-a75e-08d418bdfc98
let
	if !@isdefined(RankOneMatrix)
		not_defined(:RankOneMatrix)
	else
		R2 = RankOneMatrix([1,2], [3,4,5])

		if ismissing(Base.size(R2))
			still_missing()
		elseif size(R2) != (2, 3)
			keep_working(md"`Base.size` is not implemented correctly.")
		else
			correct()
		end
	end
end

# ╔═╡ 47f17b30-d87b-4d71-abf1-7a8b7f938095
let
	if !@isdefined(RankOneMatrix)
		not_defined(:RankOneMatrix)
	else
		R2 = RankOneMatrix([1,2], [3,4,5])

		if ismissing(Base.getindex(R2, 1, 1))
			still_missing(md"Implement both `Base.size` and `Base.getindex`.")
		elseif ismissing(Base.size(R2))
			still_missing(md"Implement both `Base.size` and `Base.getindex`.")
		elseif size(R2) != (2, 3)
			keep_working(md"`Base.size` is not implemented correctly.")
		else
			shouldbe = [
					3 4 5
					6 8 10
				]
			r = all(CartesianIndices(shouldbe)) do I
				isapprox(R2[Tuple(I)...], shouldbe[I], rtol=1-3)
			end
			if r
				correct()
			else
				keep_working(md"`Base.getindex` is not implemented correctly.")
			end
		end
	end
end

# ╔═╡ 3e72f920-8739-11eb-1354-57b18fa53103
if !@isdefined(matvec)
	not_defined(:matvec)
else
	let
		M = RankOneMatrix(1:10)

		if ismissing(matvec(M, ones(10)))
			still_missing()
		end

		incorrect = false
		for i in 1:5
			r = ones(10).*i
			for j in 1:10
				if matvec(M, r)[j] != (collect(M) * r)[j]
					incorrect = true
				end
			end
		end

		if incorrect
			keep_working()
		else
			correct()
		end
	end
end

# ╔═╡ 0d81b3d3-cfb4-4c27-bb00-18a55e00ab7a
let
	if !@isdefined(RankTwoMatrix)
		not_defined(:RankTwoMatrix)
	else
		R1 = RankOneMatrix([1,2], [3,4,5])
		R2 = RankTwoMatrix(R1, R1)

		if ismissing(Base.size(R2))
			still_missing()
		elseif size(R2) != (2, 3)
			keep_working(md"`Base.size` is not implemented correctly.")
		else
			correct()
		end
	end
end

# ╔═╡ 2e2f503f-c1ad-4d7f-8c0d-cd41fad1f334
let
	if !@isdefined(RankTwoMatrix)
		not_defined(:RankTwoMatrix)
	else
		R1a = RankOneMatrix([1,2], [3,4,5])
		R1b = RankOneMatrix([2,1], [3,4,5])
		R2 = RankTwoMatrix(R1a, R1b)

		if ismissing(Base.getindex(R2, 1, 1))
			still_missing(md"Implement both `Base.size` and `Base.getindex`.")
		elseif ismissing(Base.size(R2))
			still_missing(md"Implement both `Base.size` and `Base.getindex`.")
		elseif size(R2) != (2, 3)
			keep_working(md"`Base.size` is not implemented correctly.")
		else
			shouldbe = [
					3 4 5
					6 8 10
				]
			s2 = shouldbe + shouldbe[[2,1],:]
			r = all(CartesianIndices(s2)) do I
				isapprox(R2[Tuple(I)...], s2[I], rtol=1-3)
			end
			if r
				correct()
			else
				keep_working(md"`Base.getindex` is not implemented correctly.")
			end
		end
	end
end

# ╔═╡ 84d7bc07-e211-4296-98e6-f03723481e37
let
	if !@isdefined(LowRankMatrix)
		not_defined(:LowRankMatrix)
	else
		R1 = RankOneMatrix([1,2], [3,4,5])
		R2 = LowRankMatrix([R1, R1],2)

		if ismissing(Base.size(R2))
			still_missing()
		elseif size(R2) != (2, 3)
			keep_working(md"`Base.size` is not implemented correctly.")
		else
			correct()
		end
	end
end

# ╔═╡ 4ed7c189-3130-4f6a-b57b-5322fdd4cef6
let
	if !@isdefined(LowRankMatrix)
		not_defined(:LowRankMatrix)
	else
		R1a = RankOneMatrix([1,2], [3,4,5])
		R1b = RankOneMatrix([2,1], [3,4,5])
		R2 = LowRankMatrix([R1a, R1b],2)

		if ismissing(Base.getindex(R2, 1, 1))
			still_missing(md"Implement both `Base.size` and `Base.getindex`.")
		elseif ismissing(Base.size(R2))
			still_missing(md"Implement both `Base.size` and `Base.getindex`.")
		elseif size(R2) != (2, 3)
			keep_working(md"`Base.size` is not implemented correctly.")
		else
			shouldbe = [
					3 4 5
					6 8 10
				]
			s2 = shouldbe + shouldbe[[2,1],:]
			r = all(CartesianIndices(s2)) do I
				isapprox(R2[Tuple(I)...], s2[I], rtol=1-3)
			end
			if r
				correct()
			else
				keep_working(md"`Base.getindex` is not implemented correctly.")
			end
		end
	end
end

# ╔═╡ 863f4490-8b06-11eb-352b-61cec188ae57
if !@isdefined(matvec)
	not_defined(:matvec)
else

	let

		comp1, comp2 = RankOneMatrix(1.0:1.0:10.0), RankOneMatrix(0.0:0.1:0.9)
		if ismissing(comp1) | ismissing(comp2)
			still_missing()
		end

		ex4 = LowRankMatrix([comp1, comp2], 2)
		if ismissing(ex4)
			still_missing
		end

		if ismissing(matvec(ex4, ones(10)))
			still_missing()
		end

		incorrect = false
		for i in 1:10
			r = ones(10).*i
			v, c = matvec(ex4, r), collect(ex4) * r
			for j in 1:10
				if abs(v[j] - c[j]) > 1e-10
					incorrect = true
					break
				end
			end
		end

		if incorrect
			keep_working()
		else
			correct()
		end
	end
end

# ╔═╡ bec548dd-e2b9-4057-b0ac-1ab4e5b441d4
if !@isdefined(singular_values_of_biggie)
	not_defined(:singular_values_of_biggie)
else
	if singular_values_of_biggie isa Missing
		still_missing()
	elseif !(singular_values_of_biggie isa AbstractVector)
		keep_working(md"Return the singular values as a _vector_.")
	elseif isapprox(singular_values_of_biggie, svdvals(biggie), rtol=1e-5)
		correct()
	else
		keep_working()
	end
end

# ╔═╡ 2f5eadb0-4fdb-45b6-83fd-116dd4c1f9be
let
	if !@isdefined(numerical_rank)
		not_defined(:numerical_rank)
	else
		result1 = numerical_rank([
				1 2 0
				0 0 1e-10
				-5 -10 0
				])

		if result1 isa Missing
			still_missing()
		elseif !(result1 isa Integer)
			keep_working(md"You should return a number.")
		elseif result1 != 1
			keep_working()
		else
			result2 = numerical_rank([
				1 3 0
				0 0 1e-10
				-5 -10 0
				])
			if result2 == 2
				correct()
			else
				keep_working()
			end
		end
	end
end

# ╔═╡ 1276103b-8c58-4757-ae80-5ffb7f870d09
let
	if !@isdefined(k_rank_ones)
		not_defined(:k_rank_ones)
	else
		result = k_rank_ones(1, 3, 4)

		if result isa Missing
			still_missing()
		elseif !(result isa AbstractMatrix)
			keep_working(md"Make sure that you return a matrix.")
		elseif size(result) != (3, 4)
			keep_working(md"Make sure that you return a matrix of the correct size.")
		else

			numerical_rank(A) = count(>(1e-5), svdvals(A))

			if numerical_rank(result) == 1
				correct()
			else
				result2 = k_rank_ones(2, 10, 4)

				if !(result2 isa AbstractMatrix)
					keep_working(md"Make sure that you return a matrix.")
				elseif size(result2) != (10, 4)
					keep_working(md"Make sure that you return a matrix of the correct size.")
				elseif numerical_rank(result2) == 2
					correct()
				else
					keep_working()
				end
			end

		end
	end
end

# ╔═╡ ddf2a828-7823-4fc0-b944-72b60b391247
todo(text) = HTML("""<div
	style="background: rgb(220, 200, 255); padding: 2em; border-radius: 1em;"
	><h1>TODO</h1>$(repr(MIME"text/html"(), text))</div>""")

# ╔═╡ a7feaaa4-618b-4afe-8050-4bf7cc609c17
bigbreak = html"<br><br><br><br><br>";

# ╔═╡ 4ce0e43d-63d6-4cb2-88b8-8ba80e17012a
bigbreak

# ╔═╡ 36a63540-8739-11eb-21b2-a9d5ee60a0e5
bigbreak

# ╔═╡ 35d22fd9-9cfc-47c5-8adb-a01f14586be3
bigbreak

# ╔═╡ 0d1dac23-bdbd-4a90-ab13-e4c33d83777e
bigbreak

# ╔═╡ a5234680-8b02-11eb-2574-15489d0d49ea
bigbreak

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
"""

# ╔═╡ Cell order:
# ╟─eadb174e-2c1d-48c8-9de2-99cdc2b38d32
# ╟─0560cf7b-9986-402a-9c40-779ea7a7292d
# ╟─082542fe-f89d-4774-be20-1e1a78f21291
# ╠═6f4274b5-87e2-420d-83d2-83a8408650cd
# ╟─aaa41509-a62d-417b-bca7-a120e3a5e5b2
# ╠═8c8388cf-9891-423c-8db2-d40d870bb38e
# ╟─4ce0e43d-63d6-4cb2-88b8-8ba80e17012a
# ╟─ac35c6b3-a142-4b7a-921b-a2bfe1d84713
# ╠═ea28bf57-ba62-41ce-8be6-d38ca2c5caa3
# ╟─a759c4d5-ca00-4dda-b39c-1093ba9bef0e
# ╟─a5901f93-007f-4a30-97fc-29b367ec47c6
# ╠═f0e59676-fd8e-4609-833e-f7d1a5a6e9e3
# ╠═a5be0775-68de-41ce-95cd-1465723d099b
# ╠═4abec609-09cd-4f86-8d86-a7d02325cc7b
# ╟─32e073bb-943f-4fa9-b15f-5ec18feecf15
# ╟─7d1e1724-cd1a-483b-af40-f24ae5301849
# ╠═d221c61c-a4ab-4a82-b89d-52735d957cae
# ╠═c3d49b98-a9c1-4aba-becc-7fa84f4bbc75
# ╠═7b33c09c-2ef2-4b97-b5a5-5fdf9268f76b
# ╟─224acb51-c0db-4f04-b865-8349d6a28e98
# ╟─a7e9ab20-8501-11eb-2a7a-af8fd495b21e
# ╟─8ce8c4bc-8505-11eb-2357-c50a70b8745c
# ╠═615b88b6-8505-11eb-0a7c-294ee5ae7474
# ╟─e9a53ca0-86f2-11eb-2970-5d0996d1603e
# ╟─9f534aa0-8505-11eb-3476-4d326fb19d42
# ╠═b0db7388-850c-11eb-0915-597f1fa5ab93
# ╟─30e5a813-9ea1-48d9-9039-012b8b267256
# ╟─7aac43b3-49ce-46ea-854b-d292ac0591c7
# ╠═2dfc2ef9-3efa-4b82-960f-f2ef0171e9eb
# ╠═ead4c008-0e7e-4414-aced-d4a576423bd3
# ╠═d4a0baea-cc02-4d91-a938-8181dea7b47a
# ╟─ce6e648c-ade5-4882-9f92-1ae3f878a5a3
# ╠═7adbcbab-66e7-4b58-b197-5c220ba1e9a2
# ╟─87fce3d8-765c-410c-a059-4dd51bd97254
# ╟─e86d8ab6-8505-11eb-1995-137dfbff0191
# ╠═fada4734-8505-11eb-3f2b-d1f1ef391ba4
# ╠═0b7c6cbe-57de-419d-adcb-8724791f9c89
# ╟─06040694-b324-4e6b-8314-3bffc7b8010f
# ╟─545a5916-8506-11eb-1a31-078957d791f8
# ╠═590dfe1a-8506-11eb-0069-d7cd91f02a65
# ╠═fe6df9bf-6059-4b76-af39-385d395ece72
# ╟─4ba8a0f2-850b-11eb-1c09-b7bba1763ed7
# ╟─b9cb6192-8505-11eb-3d2c-790654bbc9a1
# ╠═941b6b10-8ae9-11eb-3bf8-b732f5f60af3
# ╟─1b3cc404-9da8-4acd-9df5-5815ce167c34
# ╠═a243a400-8af3-11eb-0637-cf8f80aae86d
# ╠═aa756549-1b9f-41a1-a2e1-713166e4211a
# ╟─4e926bfe-8734-11eb-1355-e1d1d9b36929
# ╠═c7a15c5e-8505-11eb-3af2-2fa84b74b590
# ╟─9eb846c0-8737-11eb-101b-0191f715d8c9
# ╠═b0a0b9a4-850b-11eb-30f1-11f5270efe02
# ╟─cd02acd0-8af4-11eb-3726-d7a360b98433
# ╟─5e85ac82-8734-11eb-2f69-cda569326f68
# ╟─a75f3c76-8506-11eb-1d83-cf5781b656a3
# ╠═b577420c-8501-11eb-267a-719125315fe1
# ╟─1d16e228-850e-11eb-1dca-41df85b706da
# ╠═eb612772-8b06-44bb-a36a-827435cbb2ee
# ╠═d3e7c8d1-4b4a-47b6-9c96-150333078f42
# ╟─7493be50-8738-11eb-37e7-990f1656998a
# ╟─3ed18ba0-850b-11eb-3e90-f188c54d9ce9
# ╠═f2d8b45c-8501-11eb-1c6a-5f819c240d9d
# ╠═ed72e880-8afa-11eb-3a4a-175a838188d9
# ╠═7b3fb0ef-9a9e-401c-8c09-e5615134a4ad
# ╠═5a6dd3bc-af85-4e7d-b03c-b8c1bc88d850
# ╠═52451232-2c3d-4425-a94a-da1a955cf784
# ╟─111a8db4-a808-4f37-a75e-08d418bdfc98
# ╟─47f17b30-d87b-4d71-abf1-7a8b7f938095
# ╟─91b1125a-850b-11eb-1a21-c7af293332fd
# ╠═fc962c72-8501-11eb-2821-cbb7a52d5f61
# ╟─8f5c7c9c-850b-11eb-3d97-bf9c6b5d7d2e
# ╟─248ab23e-850c-11eb-396f-5ba365f0f743
# ╠═0887ee78-0e8a-41c7-90e7-44237acc1477
# ╠═1705dc7b-9f93-440b-acbb-ecb362d08125
# ╟─8cab8d7d-32fc-4b43-bd54-fb0c78fa5d00
# ╟─6a08cb74-8510-11eb-2212-2f1f676f3d5d
# ╟─61cbfc80-86f5-11eb-316f-6307cac7cbbb
# ╠═ee58251a-8511-11eb-074c-5b1e27c4ebd4
# ╟─3e72f920-8739-11eb-1354-57b18fa53103
# ╟─36a63540-8739-11eb-21b2-a9d5ee60a0e5
# ╟─4864e4d6-850c-11eb-210c-0318b8660a9a
# ╟─58ce8d2a-850c-11eb-29e9-3fc966e129b1
# ╟─9d4d943e-8502-11eb-31cb-e11824b2dce0
# ╠═ba3cb45a-8502-11eb-2141-6369b0e08807
# ╟─92a91904-850c-11eb-010e-c58ae218f541
# ╠═f5a95dd8-850d-11eb-2aa7-2dcb1868577f
# ╟─dc714540-8afa-11eb-205b-7770074771c8
# ╠═c784e02c-8502-11eb-3efa-7f4c45f4274c
# ╠═0bab818e-8503-11eb-02b3-178098599847
# ╟─0d81b3d3-cfb4-4c27-bb00-18a55e00ab7a
# ╟─2e2f503f-c1ad-4d7f-8c0d-cd41fad1f334
# ╟─aca709f0-8503-11eb-1144-1fc01cb85c39
# ╠═b6717812-8503-11eb-2729-39bfdc1fd2f9
# ╟─6e280730-8b00-11eb-04b7-615a389e792f
# ╠═c49e350e-8503-11eb-15de-7308dd03dc08
# ╠═dd27f508-8503-11eb-36b9-33f5f99f78b0
# ╠═93a75ea0-8b0c-11eb-3946-3d5d92487fb5
# ╟─84d7bc07-e211-4296-98e6-f03723481e37
# ╟─4ed7c189-3130-4f6a-b57b-5322fdd4cef6
# ╟─1083d5ee-8512-11eb-0fdc-ab3987f94bc5
# ╠═3fed837a-8512-11eb-1fdd-c1b72b48d07b
# ╟─863f4490-8b06-11eb-352b-61cec188ae57
# ╟─f7fb3b32-8b00-11eb-0c0e-8f8d68f849a6
# ╟─2d65bd1a-8512-11eb-1bd2-0313588dfa0e
# ╟─f9556098-8504-11eb-08a0-39fbe00892da
# ╟─35d22fd9-9cfc-47c5-8adb-a01f14586be3
# ╟─295acdac-880a-402e-9f7e-19b0fc801130
# ╠═54778212-3702-467d-8a96-4fa18b3ccd63
# ╟─3ba3ae50-b9ed-416e-bb47-a14f82eb4e19
# ╠═bdf3954e-2a05-4a0d-899b-1f6f59119550
# ╠═167b9580-fa18-4248-8643-a0fde723ecc4
# ╟─a58c1c54-4f17-4ca4-ac78-22773b526dfb
# ╠═210392ff-0a22-4e55-be08-9f58804282cf
# ╟─bec548dd-e2b9-4057-b0ac-1ab4e5b441d4
# ╟─0d1dac23-bdbd-4a90-ab13-e4c33d83777e
# ╟─bb649c89-709c-49c8-8111-53044e8e682a
# ╠═6f78b3bb-5440-4f34-849a-fa8e21972fa8
# ╠═bbc4a595-a63e-4a12-9690-f39a23b30ae1
# ╟─44b8fbf1-411f-4ce4-b5d7-e12f2d6be51b
# ╠═6c9ae344-084e-459c-841c-8377451507fd
# ╠═9ae28bc5-d9cb-479d-8c22-7c9248cd5fa0
# ╟─29db89d6-ac40-4121-b996-e083555bc5db
# ╟─0cf29d7e-bd83-4c5f-bd41-e8e7c2587188
# ╠═ed4b7ea5-3266-466d-a817-7ab2cc82ac9c
# ╟─5c6aee48-017b-49ff-af1a-a31af30a45a9
# ╠═a2b7f0c3-488b-4ce9-aed6-b8ccddac6a57
# ╠═d3420859-d558-47cb-aaf7-d51a5e2d1f6e
# ╠═207d2f4c-cbe8-4828-ab3f-100809bb93e9
# ╟─2f5eadb0-4fdb-45b6-83fd-116dd4c1f9be
# ╟─7fdc110a-a3e3-44e5-a547-b75e03e0d21e
# ╠═01b12200-2e7e-4f19-96b9-5b6d6cb03233
# ╠═e1e1067b-93ba-40df-bd09-7599538e6181
# ╟─1276103b-8c58-4757-ae80-5ffb7f870d09
# ╟─00f5bea8-b808-483a-ad70-332f521481f5
# ╠═a6bb92f9-0cb2-4fdd-8b67-f286edbbdcb6
# ╠═57c19601-122b-414f-bc99-56f98c794e61
# ╠═9813e9f9-9970-4a8b-b3ca-e6a699e6fda4
# ╠═c8829a12-917a-4ff4-87c9-fe2b25aaa99c
# ╟─a5234680-8b02-11eb-2574-15489d0d49ea
# ╟─5aabbec1-a079-4936-9cd1-9c25fe5700e6
# ╟─42d6b87d-b4c3-4ccb-aceb-d1b020135f47
# ╟─8fd9e34c-8f20-4979-8f8a-e3e2ae7d6c65
# ╟─d375c52d-2126-4594-b819-aba9d34fe77f
# ╟─2d7565d4-88da-4e41-aad6-958ed6b3b2e0
# ╟─d7f45d51-f426-4353-af58-858395295ce0
# ╟─9b2aa95d-4583-40ec-893c-9fc751ea22a1
# ╟─5dca2bb3-24e6-49ae-a6fc-d3912da4f82a
# ╟─ddf2a828-7823-4fc0-b944-72b60b391247
# ╟─a7feaaa4-618b-4afe-8050-4bf7cc609c17
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
#########################################################################
#    Helper functions for derivation tests
#########################################################################

# FIXME : remove calls to getSymbols and SubstSymbols

## error thresholds
DIFF_DELTA = 1e-9
ERROR_THRESHOLD = 2e-2

good_enough(x,y) = isfinite(x) ? (abs(x-y) / max(ERROR_THRESHOLD, abs(x))) < ERROR_THRESHOLD : isequal(x,y)
good_enough(t::Tuple) = good_enough(t[1], t[2])

##  gradient check by comparing numerical gradient to automated gradient
function deriv1(ex::Expr, x0::Union(Float64, Vector{Float64}, Matrix{Float64}))

	nx = length(x0)

	ex2 = :( $(MCMC.ACC_SYM) += $ex)
	fex2 = MCMC.generateModelFunction(ex2, gradient=true, debug=true, x=x0)
	myf, dummy = MCMC.generateModelFunction(ex2, gradient=true, x=x0)

	pars0 = vec([x0])
	l0, grad0 = myf(pars0)
	gradn = Array(Float64,nx)
	for i in 1:nx
		l, grad = myf( Float64[ pars0[j] + (j==i)*DIFF_DELTA for j in 1:nx] )
		gradn[i] = (l-l0)/DIFF_DELTA
	end

	if ! all(good_enough, zip([grad0], [gradn]))
		println("Gradient false for $ex at x=$x0, expected $(round(gradn,5)), got $(round(grad0,5))")
		println()
		println( fex2 )
		error()
	end
end

##  tests derivation on all parameters, for all combinations of arguments dimension
macro test_combin(func::Expr, constraints...)
	constraints = collect(constraints)
	parnames = collect(MCMC.ReverseDiffSource.getSymbols(func))
	# println("parnames : $parnames")
	# dump(constraints)
	#  args to derive against
	dargs = filter(ex->isa(ex,Symbol), constraints)
	length(dargs)==0 && (dargs = parnames)
	@assert all(map(x->in(x,parnames), dargs)) "some of specified derivation args $(dargs) not found in tested expression $func"
	# println("dargs : $dargs")

	#  transformations on args to have valid calls
	trans = filter(ex->isa(ex, Expr) && ex.head == :(->), constraints)
	# println("trans : $trans")

	#  dimension validity rules
	rules = filter(ex->isa(ex, Expr) && ex.head != :(->), constraints)
	# println("rules : $rules")

	arity = length(parnames)
	par = [ symbol("arg$i") for i in 1:arity]

	# try each arg dim combination
	combin = [ [:v0ref, :v1ref, :v2ref][1+ ifloor((i-1) / 3^(j-1)) % 3] for i in 1:3^arity, j in 1:arity]
	for ic in 1:size(combin,1)  # try each arg dim in combin
		# create variables
		for i in 1:arity
			eval(:( $(par[i]) = copy($(combin[ic,i]))) )
		end

		# reject combination if one of rules fails
		if all( [ eval( MCMC.ReverseDiffSource.substSymbols(r, Dict(parnames, par))) for r in rules] )
			# println("## combin = $(combin[ic,:])")
			# apply transformations on args
			for t in trans
				pos = find(parnames .== t.args[1]) # find the arg the rules applies to
				@assert length(pos)>0 "arg of transfo ($t.args[1]) not found among $parnames"
				vn = symbol("arg$(pos[1])")
				eval(:( $vn = map($t, $vn)))
			end

			# try derivation for each allowed argument
			for p in dargs
				tpar = copy(par)
				tpar[p .== parnames] = :x  # replace tested args with parameter symbol :x for deriv1 testing func
				fex = MCMC.ReverseDiffSource.substSymbols(func, Dict( parnames, tpar))
				# println("##+## $fex  $(Dict( parnames, tpar))")
				x0 = eval(par[p .== parnames][1])  # set x0 for deriv 1
				# println("##-## $fex  ##-##   $x0")
				deriv1(fex, x0)
			end
		end
	end

end