text:
type ModVar{T,S}
    name::Symbol
    def::Expr
    val::T
    init::Float64
    bounds::Vector{Float64}
end

abstract Agg
abstract Indv
abstract ModelVariable
    abstract State <: ModelVariable
        type Endogenous <: State
            x::Vector{Float64}
            Endogenous(x::Vector{Float64})=new(x)
        end
        abstract Stochastic <: State
            type AR <: Stochastic
                μ::Float64
                ρ::Float64
                σ::Float64
                x::Vector{Float64}
                T::Array{Float64,2}
            end
            type Markov <: Stochastic
                x::Vector{Float64}
                T::Array{Float64,2}
                Markov(x::Vector{Float64},T::Array{Float64,2}) = new(x,T)
            end
    type Policy <: ModelVariable end
    abstract Static <: ModelVariable
        type Dependant  <: Static end
        type Exogenous  <: Static end
    type Aggregate <: ModelVariable
        target::Tuple{Symbol,Int}
        Aggregate(s::Symbol,t::Int)=new((s,t))
    end
    type Future <: ModelVariable
        targname::Symbol
        targXid::Int
        Future(targname::Symbol,targXid::Int) = new(targname,targXid)
    end

function (::Type{Endogenous})(ex::Expr)
    return  ModVar{Endogenous,Indv}(
                ex.args[1],
                ex.args[2],
                Endogenous(CC.g(ex.args[2].args[3]+1)*(ex.args[2].args[2]-ex.args[2].args[1])+ex.args[2].args[1]),
                0.0,
                ex.args[2].args[1:2])
end

function (::Type{AR})(ex::Expr)
    x,T = rouwenhorst(ex.args[2].args[1:3]...,CC.M(1+ex.args[2].args[4]))

    ModVar{AR,isag(ex)}(
        ex.args[1],
        ex.args[2],
        AR(ex.args[2].args[1:3]...,x,T),
        0.0,
        collect(extrema(x)))
end

function (::Type{Markov})(ex::Expr)
    x = ex.args[2].args[1]
    x = collect(linspace(x.args[1],x.args[2],CC.M(1+ex.args[2].args[3])))
    T = eval(current_module(),ex.args[2].args[2])
    @assert length(x) == size(T,1) == size(T,2)

    ModVar{Markov,isag(ex)}(
        ex.args[1],
        ex.args[2],
        Markov(x,T),
        0.0,
        collect(extrema(x)))
end

function (::Type{Policy})(ex::Expr)
    ModVar{Policy,Indv}(
        ex.args[1],
        ex.args[2],
        Policy(),
        ex.args[2].args[end],
        map(Float64,ex.args[2].args[1:2]))
end

function (::Type{Dependant})(ex::Expr)
    ModVar{Dependant,Indv}(
        ex.args[1],
        ex.args[2],
        Dependant(),
        0.0,
        [-Inf,Inf])
end

function (::Type{Exogenous})(ex::Expr)
    ModVar{Exogenous,Indv}(
        ex.args[1],
        :([$(ex.args[2])]),
        Exogenous(),
        ex.args[2],
        [-Inf,Inf])
end

function (::Type{Aggregate})(ex::Expr,slist::Vector{ModVar},dlist::Vector{ModVar})
    if isa(ex.args[2].args[2],Symbol)
        if in(ex.args[2].args[2],names(Endogenous(slist)))
            return [ModVar{Aggregate,Agg}(
                        ex.args[1],
                        ex.args[2],
                        Aggregate(ex.args[2].args[2],-1),
                        ex.args[2].args[3],
                        [-Inf,Inf])
                    ModVar{Aggregate,Agg}(
                                ex.args[1],
                                ex.args[2],
                                Aggregate(ex.args[2].args[2],0),
                                ex.args[2].args[3],
                                [-Inf,Inf])]
        else
            return ModVar{Aggregate,Agg}(
                        ex.args[1],
                        ex.args[2],
                        Aggregate(ex.args[2].args[2],0),
                        ex.args[2].args[3],
                        [-Inf,Inf])
        end

    elseif isa(ex.args[2].args[2],Expr)
        dex = :($(Symbol(string(ex.args[1])*"targ")) = $(ex.args[2].args[2]))
        t = maximum([x.args[2] for x in unique(getv(addindex(dex.args[2])))])
        push!(dlist,Dependant(dex))

        return ModVar{Aggregate,Agg}(
                    ex.args[1],
                    dex.args[2],
                    Aggregate(Symbol(string(ex.args[1])*"targ"),t),
                    ex.args[2].args[3],
                    [-Inf,Inf])
    else
        error("Integrand in $(ex.args[1]) is neither a symbol nor an expression")
    end
end

function (::Type{Future})(f2::Expr,slist::Vector{ModVar},vlist::Vector{ModVar})
    futurevar = Symbol[x.args[1] for x in filter(x->x.args[2]>0,getv(f2))]
    filter!(x->!in(x,names(Stochastic(slist))),futurevar)
    list = ModVar[]

    return [ModVar{Future,Indv}(v,Expr(:ref,v,1),Future(v,find(State(vlist,false),v)),0.0,State(vlist,false)[v].bounds) for v in futurevar]
end

length{T<:State}(s::T) = length(s.x)

ARSim(x::Vector{Float64},e::AR,s=randn(size(x))) = clamp((1-e.ρ)*e.μ + e.ρ*x+s*e.σ,e.x[1],e.x[end])
function MarkovSim!(ID::Vector{Int},s::Stochastic,r::Vector{Float64}=rand(length(ID)))
    csT=  cumsum(s.T,2)
    for i = 1:length(ID)
        j = 1
        while csT[ID[i],j]<r[i]
            j+=1
        end
        ID[i] = j
    end
    return ID
end

(s::AR)(ID::Vector{Int},r::Vector{Float64}=rand(length(ID))) = MarkovSim!(ID,s,r)
(s::AR)(ID::Int,r::Vector{Float64}=rand(length(ID))) = MarkovSim!([ID],s,r)[1]
(s::Markov)(ID::Vector{Int},r::Vector{Float64}=rand(length(ID))) = MarkovSim!(ID,s,r)

function getindex(list::Vector{ModVar},v::Symbol)
    for i = 1:length(list)
        list[i].name == v && (return list[i])
    end
end

function find(list::Vector{ModVar},v::Symbol)
    for i = 1:length(list)
        list[i].name == v && (return i)
    end
end

display(V::Vector{ModVar}) = [println(x.name,"\t",split(string(vtype(x)),'.')[end]) for x in V]
names(V::Vector{ModVar}) = map(x->x.name,V)
length{T<:Stochastic,S}(s::ModVar{T,S}) = length(s.val.x)
vtype{T,S}(x::ModVar{T,S})=(T)
isag{T,S}(x::ModVar{T,S})=(S==Agg)

(::Type{T}){T<:ModelVariable}(V::Array{ModVar},B::Bool=true) = filter(x->B ? (vtype(x)==T || vtype(x)<:T) : !(vtype(x)==T || vtype(x)<:T),V)::Vector{ModVar}

isag(v::Expr) = v.head==:(=) ? Indv : Agg

sort(list::Vector{ModVar}) = vcat(State(list),Policy(list),Static(list),Aggregate(list),Future(list))

function parsestate(states)
    slist = ModVar[]
    for v ∈ states.args
        if length(v.args[2].args)==3 && reduce(&,map(x->isa(x,Number),v.args[2].args))
            push!(slist,Endogenous(v))
        elseif length(v.args[2].args)==3
            push!(slist,Markov(v))
        elseif length(v.args[2].args)==4 && reduce(&,map(x->isa(x,Number),v.args[2].args))
            push!(slist,AR(v))
        end
    end
    slist = vcat(Endogenous(slist),Stochastic(slist))

    G = NGrid(CC,Int[x.def.args[end] for x in slist],hcat(Vector{Float64}[x.bounds for x in slist]...),B=QuadraticBF)
	for s in Stochastic(slist)
        s.val.x = sort(unique(values(G)[:,s.name.==names(slist)]))
	end
    S = values(G)

    exogenousfuture=map(x->repmat(x[:]',length(G),1),ndgrid([s.val.x for s in Stochastic(slist)]...))
    nP=prod(map(length,Stochastic(slist)))
    ProbWeights =  Float64[prod([Stochastic(slist)[N].val.T[findfirst(Stochastic(slist)[N].val.x,S[i,findfirst(names(State(slist)),names(Stochastic(slist))[N])]),findfirst(Stochastic(slist)[N].val.x,exogenousfuture[N][i,j])] for N = 1:length(Stochastic(slist))]) for i = 1:length(G),j=1:nP]
    SP = zeros(length(G)*nP,length(slist))
    for i = 1:length(Stochastic(slist))
        SP[:,findfirst(names(slist),Stochastic(slist)[i].name)] = exogenousfuture[i][:]
    end

    slist,G,S,SP,ProbWeights
end

function parsevars(slist,policy,vars,params)
    vlist = ModVar[]
    for v ∈ policy.args
		push!(vlist,Policy(v))
	end

    dlist = ModVar[]
    for v ∈ vars.args
        @assert v.head==:(=)
        if isa(v.args[2],Expr) && v.args[2].args[1]!=:∫
            push!(dlist,Dependant(v))
        elseif isa(v.args[2],Float64)
            push!(vlist,Exogenous(v))
        elseif isa(v.args[2],Expr) && v.args[2].args[1]==:∫
            vlist = vcat(vlist,Aggregate(v,slist,dlist))
        else
            warn("Failed to parse $(v.args[1])")
        end
    end

    plist   = Dict{Symbol,Float64}(zip([x.args[1] for x in params.args],[x.args[2] for x in params.args]))
    dlist1 = Dict()
	for i = 1:length(dlist)
		p = addindex!(dlist[i].name)=>subs(addindex!(subs(dlist[i].def,plist)),dlist1)
		push!(dlist1,p)
		push!(dlist1,tchange(p[1],1)=>tchange(p[2],1))
	end

    return vlist,dlist,dlist1,plist
end

function (::Type{Dependant})(vlist::Vector{ModVar},plist::Dict{Symbol,Float64})
    list = Dependant(vlist)
	dlist1 = Dict()
	for i = 1:length(list)
		p = addindex!(list[i].name)=>subs(addindex!(subs(list[i].def,plist)),dlist)
		push!(dlist,p)
		push!(dlist,tchange(p[1],1)=>tchange(p[2],1))
	end
	return dlist
end

module ProblemParse
using Calculus
import Base:Expr,parse
import EconModel:subs!,operators,addindex!,tchange

function parse(O_::Expr,C_::Expr,controls_::Expr,params::Dict)

    O = deepcopy(O_)
    C = deepcopy(C_)
    controls = addindex!(deepcopy(controls_))
    for v in (O,C)
        subs!(v,params)
        addindex!(v)
    end
    CP = tchange(C,1)

    focs = :[]
    for i = 1:length(C.args)
        c=C.args[i]
        cp=CP.args[i]
        if c.head==:comparison && c.args[2]==:≤
            C.args[i] = c.args[1]==0 ? c.args[3]+0 : c.args[3]-c.args[1]
            CP.args[i] = cp.args[1]==0 ? cp.args[3]+0 : cp.args[3]-cp.args[1]
            push!(focs.args,:(($(symbol("λ$(i)"))^2+($(c.args[3]-c.args[1]))^2)^(1/2)-$(symbol("λ$(i)"))-($(c.args[3]-c.args[1]))))
        else
            push!(focs.args,c)
        end
    end

    for x in controls.args
        focO = differentiate(O,x)
        focC = +([Expr(:ref,symbol("λ$i"),0)*differentiate(C.args[i],x) for i = 1:length(C.args)]...)
        focCP = simplify(+([Expr(:ref,symbol("λ$(i)"),1)*differentiate(CP.args[i],x) for i = 1:length(C.args)]...))
        (focC!=0 && focCP!=0) ? (focC=focC*x;focCP=focCP*x) : nothing
        focCP = focCP==0 ? focCP : params[:β]*Expr(:call,:Expect,focCP)
        push!(focs.args,simplify(focO+focC+focCP))
    end

    return focs
end

end


text:
using UniversalLogger
using Test
using Logging

@testset "UniversalLogger.jl" begin

    let logger = ULogger()
        with_logger(logger) do
            @info "a" x = zeros(5, 5) y = "a1"
            @info :b x = ones(5, 5) y = "b1"

            @info "a" x = zeros(10, 10) y = "a2"
            @info :b x = ones(10, 10) y = "b2"

        end

        lg = logger.log

        @test lg["a"][:x] == [zeros(5, 5), zeros(10, 10)]
        @test lg["a"][:y] == ["a1", "a2"]
        @test lg[:b][:x] == [ones(5, 5), ones(10, 10)]
        @test lg[:b][:y] == ["b1", "b2"]
    end

    let logger = ULogger()
        with_logger(logger) do
            for i = 1:10
                @info "foo" i val1 = i * 2
                if mod(i, 2) == 0
                    @info :bar val2 = i * 3
                end
            end
        end

        log = get(logger)
        @test log["foo"][:val1] == (1:10) .* 2
        @test log[:bar][:val2] == [6, 12, 18, 24, 30]
    end

    let logger = ULogger()
        with_logger(logger) do
          for i=1:10
            @info "foo" i fooval=i*2
            if mod(i, 2) == 0
              @info :bar barval=i*3
            end
          end
        end

        push!(logger, "baz", :baz1=>1)
        push!(logger, "baz", "baz2"=>10)
        push!(logger, "baz", Dict(:baz1=>2, "baz2"=>20))

        @test logger["foo"][:i] == 1:10
        @test logger["foo"][:fooval] == [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]

        @test logger[:bar][:barval] == [6, 12, 18, 24, 30]

        @test logger["baz", :baz1] == [1, 2]
        @test logger["baz", "baz2"] == [10, 20]
    end

end

ext: jl


text:
function remove_dashes(input_flag)
    first_non_dash = findfirst(!isequal('-'), input_flag)
    return first_non_dash - 1, input_flag[first_non_dash:end]
end

function longest(strs::Vector{String})
    idx = argmax(length.(strs))
    return strs[idx]
end

function is_numeric(flag::String)
    return tryparse(Float64, flag) !== nothing
end


text:
# Scientific Pitch Notation
module SPN # Musical

export NoteOctve
export note, octave, freq
export A4

import ..Notations: Notation, A
import MIDI # MIDI.Note

struct NoteOctave
    notation::Notation
    octave::Int
    midipitch::Int
end

const A4 = NoteOctave(A, 4, 69)

note(no::NoteOctave) = no.notation
octave(no::NoteOctave) = no.octave
freq(no::NoteOctave) = 440 * 2^((no.midipitch-69)/12)
(::Type{MIDI.Note})(no::NoteOctave) = MIDI.Note(no.midipitch, 0, 0, 0)

end # module Musical.SPN


text:
const GCC = "gcc"
const build_sub = ("automotive","consumer","network","office","security","telecomm")
const parallel_build = true

max_stars_repo_name: emcores/MiBench

max_issues_repo_name: emcores/MiBench

max_forks_repo_name: emcores/MiBench


text:
using Compat

"""
    RomanNumeral(::Integer)
    RomanNumeral(::AbstractString)

Roman number: it contains both the standard number and its roman representation.
It behaves like standard numbers.
"""
@compat struct RomanNumeral <: Integer
    val::Int
    str::String
    RomanNumeral(int::Integer) = new(int, toroman(int))
    function RomanNumeral(str::AbstractString)
        num = parse(RomanNumeral, str)
        new(num, toroman(num))
    end
end

"""rn"MMXV" syntactic sugar"""
macro rn_str(str)
    RomanNumeral(str)
end

const RN = RomanNumeral

# Standard functions
# Conversion + promotion
Base.convert(::Type{Bool}, ::RN) = true
Base.convert(::Type{T}, x::RN) where {T<:Real} = T(x.val)

Base.promote_rule(::Type{RN}, ::Type{T}) where {T <: Integer} = T

# IO
Base.show(io::IO, num::RN) = write(io, num.str)
Base.length(num::RN) = length(num.str)
Base.hash(num::RN) = @compat xor(hash(num.str), hash(num.val))


text:
using SparseGrids
using Test

@testset "Unique nodes" begin
    @testset "All nodes are different" begin
        nodes = [[0;], [1;]]
        weights = [1; 1]

        un, uw = SparseGrids.uniquenodes(nodes, weights)

        @test un == nodes
        @test uw == weights
    end

    @testset "Two identical nodes" begin
        nodes = [[0;], [0;]]
        weights = [1; 1]

        un, uw = SparseGrids.uniquenodes(nodes, weights)

        @test un == unique(nodes)
        @test uw == [sum(weights)]
    end

    @testset "Two of three nodes identical" begin
        nodes = [[0;], [1;], [0;]]
        weights = [1, 2, 3]

        un, uw = SparseGrids.uniquenodes(nodes, weights)

        @test un == [[0;], [1;]]
        @test uw == [4; 2]
    end
end


text:
#' ---
#' title: Nonlinear Modelling
#' ---

#' **Originally Contributed by**: Arpit Bhatia

#' This tutorial provides a breif introduction to nonlinear modelling in JuMP.
#' For more details and specifics, visit the [JuMP docs](http://www.juliaopt.org/JuMP.jl/stable/nlp/).

#' ## Nonlinear Programs
#' While we have already seen examples of linear, quadratic and conic programs,
#' JuMP also supports other general smooth nonlinear (convex and nonconvex) optimization problems.

#' A JuMP model object can contain a mix of linear, quadratic, and nonlinear contraints or objective functions.
#' Thus, a model object for a nonlinear program is constructed in the same way as before.

using JuMP, Ipopt
model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0));

#' ### Variables
#' Variables are modelled using the `@variable` macro as usual and
#' a starting point may be provided by using the `start` keyword argument

@variable(model, x, start = 4)
@variable(model, y, start = -9.66);

#' ### Parameters
#' Only in the case of nonlinear models, JuMP offers a syntax for "parameter" objects
#' which can refer to a numerical value.

@NLparameter(model, p == 0.003); # Providing a starting value is necessary for parameters
@NLparameter(model, l[i = 1:10] == 4 - i); # A collection of parameters

#' The `value` and `set_value` functions are used to query and update the value of a parameter respectively.

value(l[1])

#+

set_value(l[1], -4)
value(l[1])

#' Parameters are useful since it's faster to modify a model in-place by changing the value of the parameter
#' compared to creating an entirely new model object.

#' ### Expressions
#' JuMP also supports the creation of arithmetic expressions which can then be inserted into
#' constraints, the objective and other expressions.

@NLexpression(model, expr_1, sin(x))
@NLexpression(model, expr_2, asin(expr_1)); # Inserting one expression into another

#' There are some [syntax rules](https://pkg.julialang.org/docs/JuMP/DmXqY/0.19.2/nlp/#Syntax-notes-1)
#' which must be followed while writing a nonlinear expression.

#' Note that JuMP also supports linear and quadratic expression.
#' You can find out more about this functionality in the [docs](https://pkg.julialang.org/docs/JuMP/DmXqY/0.19.2/expressions/).

#' ### Nonlinear Objectives and Constraints
#' Nonlinear objectives and constraints are specified by using the `@NLobjective` and `@NLconstraint` macros.

@NLconstraint(model, exp(x) + y^4 <= 0)
@NLobjective(model, Min, tan(x) + log(y))

#' ### User-defined Functions
#' In addition to supporting a library of built-in functions,
#' JuMP supports the creation of user-defined nonlinear functions to use within nonlinear expressions.
#' The `register` function is used to enable this functionality.

my_function(a,b) = (a * b)^-6 + (b / a)^3
register(model, :my_function, 2, my_function, autodiff = true)

#' The arguements to the function are:
#' - model for which the function is being registered
#' - Julia symbol object corresponding to the name of the function
#' - Number of arguments the function takes
#' - name of the Julia method
#' - instruction for JuMP to compute exact gradients automatically

#' ## MLE using JuMP

#' Since we already have a bit of JuMP experience at this point,
#' let's try a modelling example and apply what we have learnt.
#' In this example, we compute the maximum likelihood estimate (MLE) of
#' the parameters of a Gaussian distribution i.e. the sample mean and variance.

#' If $X_{1}, \ldots, X_{n}$ are an id sample from a population with pdf or pmf
#' $f\left(x | \theta_{1}, \ldots, \theta_{k}\right),$ the likelihood function is defined by

#' $$
#' L(\theta | \mathbf{x})=L\left(\theta_{1}, \ldots, \theta_{k} | x_{1}, \ldots, x_{n}\right)=\prod_{i=1}^{n} f\left(x_{i} | \theta_{1}, \ldots, \theta_{k}\right)
#' $$

#' For each sample point $\mathbf{x}$, let $\hat{\theta}(\mathbf{x})$ be a parameter value
#' at which $L(\theta | \mathbf{x})$ attains its maximum as a function of $\theta,$ with $\mathbf{x}$ held fixed.
#' A maximum likelihood estimator (MLE) of the parameter $\theta$ based on a sample $\mathbf{X}$ is
#' $\hat{\theta}(\mathbf{X})$.

#'The Gaussian likelihood is -

#' $$
#' L(\theta | \mathbf{x})=\prod_{i=1}^{n} \frac{1}{(2 \pi)^{1 / 2}} e^{-(1 / 2)\left(x_{i}-\theta\right)^{2}}=\frac{1}{(2 \pi)^{n / 2}} e^{(-1 / 2) \Sigma_{i=1}^{n}\left(x_{i}-\theta\right)^{2}}
#' $$

#' In most cases, the natural logarithm of
#' $L(\theta | \mathbf{x}), \log L(\theta | \mathbf{x})$ (known as the log likelihood),
#' is used rather than $L(\theta | \mathbf{x})$ directly.
#' The reason is that the log likelihood is easier to differentiate.
#' This substituion is possible because the log function is strictly increasing on $(0, \infty)$,
#' which implies that the extrema of $L(\theta | \mathbf{x})$ and $\log L(\theta | \mathbf{x})$ coincide.

using Random, Statistics

Random.seed!(1234)

n = 1_000
data = randn(n)

mle = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
@NLparameter(mle, problem_data[i = 1:n] == data[i])
μ0 = randn()
σ0 = rand() + 1
@info "Starting guess, mean: $μ0, std: $σ0"
@variable(mle, μ, start = μ0)
@variable(mle, σ >= 0.0, start = σ0)
@NLexpression(mle, loglikelihood,
    -(n / 2) * (log(2π) + 2 * log(σ)) - inv(2 * σ^2) * sum((xi - μ)^2 for xi in problem_data)
)

@NLobjective(mle, Max, loglikelihood)

optimize!(mle)

println("μ = ", value(μ))
println("mean(data) = ", mean(data))
println("σ^2 = ", value(σ)^2)
println("var(data) = ", var(data))
println("MLE value: ", exp(objective_value(mle)))

#+

# Changing the data

data = randn(n)
optimize!(mle)

println("μ = ", value(μ))
println("mean(data) = ", mean(data))
println("σ^2 = ", value(σ)^2)
println("var(data) = ", var(data))
println("MLE objective: ", objective_value(mle))

#' ## Writing Convex Models

#' Nonlinear solvers like Ipopt are usually local solvers.
#' For convex problems, the local optima is also the global optima,
#' and thus Ipopt is able to provide us with the correct solution.
#' However, in case a problem is not written in the convex form, then we may be unable to solve it.

#' A tool that helps us in dealing with this issue is [Disciplined Convex Programming](https://dcp.stanford.edu) (DCP).
#' DCP is a system for constructing mathematical expressions with known curvature from a given library of base functions.
#' Specifically, it helps us to construct convex optimization models when possible,
#' i.e. minimize convex function or maximize concave function and use constraints that are convex $f$ <= concave $g$


text:
abstract type AlignmentMethod end
struct ΛTAME_M <: AlignmentMethod end
struct ΛTAME_MultiMotif_M <: AlignmentMethod end
struct LowRankTAME_M <: AlignmentMethod end
struct LowRankTAME_MultiMotif_M <: AlignmentMethod end
struct TAME_M <: AlignmentMethod end
struct TAME_MultiMotif_M <: AlignmentMethod end
struct EigenAlign_M <: AlignmentMethod end
struct LowRankEigenAlign_M <: AlignmentMethod end
struct LowRankEigenAlignOnlyEdges_M <: AlignmentMethod end
struct Degree_M <: AlignmentMethod end
struct Random_M <: AlignmentMethod end

"""------------------------------------------------------------------------------
  This function takes in a pair of sparse matrices and aligns them with with the
  specified method. If the method requested is a tensor method, the adjacency
  tensor is built from the triangles enumerated. Additional kwargs are passed

  Inputs
  ------
  * A, B - (SparseMatrixCSC{T,Int}):
    The sparse matrices to be aligned.
  * profile - (Bool):
    Whether or not to profile the methods run. The type of profiling done is
    dependent on the method run.
  * method - (String):
    The method to run. Options include 'Degree', 'Random', 'LowRankEigenAlign',
    'LambdaTAME', 'LowRankTAME', "TAME". Currently supporting 'EigenAlign', but
    may remove in the future as LowRankEigenAlign produces the same result but
    with much better scalability.

  Outputs
  -------
  output depends on whether a tensor method is used or not. Matrix based methods
  return the trianlges in A, triangles in B, triangles matched, the matching
  used, runtime.

  If a tensor method is used, the number of triangles in A and B are returned, in
  addition to whatever is returned by 'align_tensors(_profiled)'.
------------------------------------------------------------------------------"""
function align_matrices(A::SparseMatrixCSC{T,Int},B::SparseMatrixCSC{S,Int};
                        postProcessing::PostProcessingMethod=noPostProcessing(),
                        profile=false,motif=Clique(),kwargs...) where {T,S}

    method = typeof(kwargs[:method])

    A_ten = graph_to_ThirdOrderTensor(A)
    B_ten = graph_to_ThirdOrderTensor(B)

    if method === ΛTAME_M || method === LowRankTAME_M
        if profile
            alignment_output = align_tensors_profiled(A_ten,B_ten;kwargs...)
        else
            alignment_output = align_tensors(A_ten,B_ten;kwargs...)
        end

        if typeof(postProcessing) === noPostProcessing
            return alignment_output
        else
            return alignment_output, post_process_alignment(A,B,A_ten,B_ten,alignment_output,postProcessing;profile,kwargs...)
        end
    elseif method === TAME_M

        if profile
            alignment_output = align_tensors_profiled(A_ten,B_ten;kwargs...)
        else
            alignment_output = align_tensors(A_ten,B_ten;kwargs...)
        end

        return alignment_output

    elseif method === ΛTAME_MultiMotif_M || method == LowRankTAME_MultiMotif_M || method == TAME_MultiMotif_M

        A_tensors = tensors_from_graph(A,kwargs[:orders],kwargs[:samples],motif)
        B_tensors = tensors_from_graph(B,kwargs[:orders],kwargs[:samples],motif)

        A_motifDistribution = [contraction(tensor,ones(tensor.n))./factorial(tensor.order-1) for tensor in A_tensors]
        B_motifDistribution = [contraction(tensor,ones(tensor.n))./factorial(tensor.order-1) for tensor in B_tensors]

        #TODO: standardize kwarg consumption
        subkwargs = Dict([(k,v) for (k,v) in kwargs if k != :orders && k != :samples])

        if profile
            alignment_output = align_tensors_profiled(A_tensors,B_tensors;subkwargs...)
        else
            alignment_output = align_tensors(A_tensors,B_tensors;subkwargs...)
        end

        if typeof(postProcessing) === noPostProcessing
            return A_motifDistribution, B_motifDistribution, alignment_output
        else
            return A_motifDistribution, B_motifDistribution, alignment_output, post_process_alignment(A,B,A_tensors,B_tensors,alignment_output,postProcessing;profile,subkwargs...)
        end
    elseif method === LowRankEigenAlign_M
            iters = 10
            (U,V, ma,mb,_,_),t = @timed align_networks_eigenalign(A,B,iters,"lowrank_svd_union",3)
            matching = Dict{Int,Int}([i=>j for (i,j) in zip(ma,mb)])
            triangle_count, _  = TAME_score(A_ten,B_ten,matching)
            alignment_output = LowRankEigenAlign_Return(triangle_count,
                                                        (size(A_ten.indices,2), size(B_ten.indices,2)),
                                                        matching,(U,V),t)

            if typeof(postProcessing) === noPostProcessing
                return alignment_output
            else
                return alignment_output, post_process_alignment(A,B,A_ten,B_ten,alignment_output,postProcessing;profile,kwargs...)
            end

    elseif method === EigenAlign_M || method === Degree_M || method === Random_M || method === LowRankEigenAlignOnlyEdges_M

        if method === LowRankEigenAlignOnlyEdges_M
            matching,t = @timed lowRankEigenAlignEdgesOnly(A,B)
        elseif method === EigenAlign_M
            (ma,mb),t = @timed NetworkAlignment.EigenAlign(A,B)
            matching = Dict{Int,Int}(zip(ma,mb))
        elseif method === Degree_M
            (ma,mb),t = @timed degree_based_matching(A,B)
            matching = Dict{Int,Int}(zip(ma,mb))
        elseif kwargs[:method] === Random_M
            n,n = size(B)

            matching,t = @timed Dict{Int,Int}(enumerate(shuffle(1:n)))
        else
            error("Invalid input, must be ")
        end

        return size(A_ten.indices,1), size(B_ten.indices,1), triangle_count, matching, t
    else
        throw(ArgumentError("method must be of type ΛTAME_M, LowRankTAME_M, TAME_M, EigenAlign_M, LowRankEigenAlign_M, LowRankEigenAlignEdgesOnly_M, Degree_M, or Random_M."))
    end
end

function align_matrices(fileA,fileB;kwargs...)
    A = MatrixNetworks.readSMAT(fileA)
    B = MatrixNetworks.readSMAT(fileB)
    return align_matrices(A,B;kwargs...)
end

struct KlauPostProcessReturn <: returnType
    original_edges_matched::Int
    matching::Union{Dict{Int,Int},Vector{Int}}
    klau_edges_matched::Int
    klau_tris_matched::Int
    setup_rt::Union{Nothing,Float64}
    klau_rt::Union{Nothing,Float64}
    L_sparsity::Float64
    f_status::Tuple{Float64,Float64}
end

function post_process_alignment(A::SparseMatrixCSC{T,Int},B::SparseMatrixCSC{S,Int},
                                A_ten,B_ten,
                                alignment_output::returnType,postProcessing::KlauAlgo;
                                profile=false,kwargs...) where {T,S}

    method = typeof(alignment_output)
    if method <: ΛTAME_Return || method <: LowRankTAME_Return || method <: ΛTAME_MultiMotif_Return || method <: LowRankEigenAlign_Return

        best_U, best_V = alignment_output.embedding
        best_matching = alignment_output.matching

        if postProcessing.k == -1
            k = size(alignment_output.embedding[1],2)*2
        else
            k = postProcessing.k
        end

        if profile
            (pp_matching,Klau_rt,L_sparsity,fvals),full_rt = @timed netalignmr(A,B,best_U, best_V, best_matching,k,postProcessing)
            setup_rt = full_rt - Klau_rt
        else
            pp_matching, Klau_rt, L_sparsity,fvals = netalignmr(A, B,best_U, best_V, best_matching,k,postProcessing)
        end

        pp_matched_motif,_= TAME_score(A_ten,B_ten,pp_matching)
        pp_matched_edges = Int(edges_matched(A,B,pp_matching)[1]/2) # code doesn't handle symmetries
        original_matched_edges = edges_matched(A,B,alignment_output.matching)[1]/2

        if profile
            return KlauPostProcessReturn(original_matched_edges,pp_matching,pp_matched_edges,pp_matched_motif, setup_rt, Klau_rt, L_sparsity,Tuple(fvals))
        else
            return KlauPostProcessReturn(original_matched_edges,pp_matching,pp_matched_edges,pp_matched_motif, nothing, nothing, L_sparsity,Tuple(fvals))
        end

    else
        throw(ArgumentError("Post processing is only supported for ΛTAME_M got $(method)"))
    end
end

struct SuccessiveKlauPostProcessReturn <: returnType
    original_edges_matched::Int
    matching::Union{Dict{Int,Int},Vector{Int}}
    klau_edges_matched::Int
    klau_tris_matched::Int
    profiling::Union{Nothing,Dict{String,Union{Vector{Float64},Vector{Int64},Vector{Vector{Float64}}}}}
end

function post_process_alignment(A::SparseMatrixCSC{T,Int},B::SparseMatrixCSC{S,Int},
                            A_ten,B_ten,
                            alignment_output::returnType,postProcessing::SuccessiveKlauAlgo;
                            profile=false,kwargs...) where {T,S}

    method = typeof(alignment_output)
    if method <: ΛTAME_Return || method <: LowRankTAME_Return || method <: ΛTAME_MultiMotif_Return || method <: LowRankEigenAlign_Return

        best_U, best_V = alignment_output.embedding
        best_matching = alignment_output.matching

        k = Int(ceil(.02*min(size(A,1),size(B,1))))
        k = 5

        if postProcessing.k == -1
            k = size(alignment_output.embedding[1],2)*2
        else
            k = postProcessing.k
        end

        if profile
            pp_matching, profiling_results = successive_netalignmr_profiled(A, B, best_U, best_V, best_matching, postProcessing)
        else
            pp_matching = successive_netalignmr(A, B, best_U, best_V, best_matching, postProcessing)
        end

        pp_matched_motif,_= TAME_score(A_ten,B_ten,pp_matching)
        pp_matched_edges = Int(edges_matched(A,B,pp_matching)[1]/2) # code doesn't handle symmetries
        original_matched_edges = Int(edges_matched(A,B,alignment_output.matching)[1]/2)

        if profile
            return SuccessiveKlauPostProcessReturn(original_matched_edges,pp_matching,pp_matched_edges,pp_matched_motif, profiling_results)
        else
            return SuccessiveKlauPostProcessReturn(original_matched_edges,pp_matching,pp_matched_edges,pp_matched_motif,nothing)
        end

    else
        throw(ArgumentError("Post processing is only supported for ΛTAME_M, ΛTAME_MultiMotif_M, LowRankTAME_M got $(method)"))
    end
end

struct LocalSearchPostProcessReturn <: returnType
    original_edges_matched::Int
    matching::Union{Dict{Int,Int},Vector{Int}}
    localsearch_edges_matched::Int
    localsearch_tris_matched::Int
    full_runtime::Union{Nothing,Float64}
    profiling::Union{Nothing,Dict{String,Union{Vector{Tuple{Int64,Int64}},Vector{Float64}}}}
end

function post_process_alignment(A::SparseMatrixCSC{T,Int},B::SparseMatrixCSC{S,Int},
                                A_ten, B_ten, alignment_output::returnType,
                                postProcessing::LocalSearch; profile=false,
                                kwargs...) where {T,S}

    if typeof(A_ten) <: Vector
        A_ten =  A_ten[end]
    end

    if typeof(B_ten) <: Vector
        B_ten =  B_ten[end]
    end

    method = typeof(kwargs[:method])
   # A_ten = graph_to_ThirdOrderTensor(A)
    # B_ten = graph_to_ThirdOrderTensor(B)
    method = typeof(alignment_output)

    if method <: ΛTAME_Return || method <: LowRankTAME_Return || method <: ΛTAME_MultiMotif_Return || method <: LowRankEigenAlign_Return

        best_U, best_V = alignment_output.embedding
        best_matching = alignment_output.matching

        if profile
            #(pp_matching,Klau_rt,L_sparsity,fvals),full_rt = @timed netalignmr(A,B,best_U, best_V, best_matching,k,postProcessing)
            #setup_rt = full_rt - Klau_rt
            (pp_matching,profile_results),full_rt = @timed local_search_profiled(A,B,A_ten, B_ten, best_U, best_V, best_matching,  postProcessing)
        else
            pp_matching = local_search(A, B, A_ten, B_ten, best_U, best_V, best_matching,  postProcessing)
        end

        pp_matched_motif,_= TAME_score(A_ten,B_ten,pp_matching)
        pp_matched_edges = Int(edges_matched(A,B,pp_matching)[1]/2) # code doesn't handle symmetries
        original_matched_edges = edges_matched(A,B,alignment_output.matching)[1]/2

        if profile
            return LocalSearchPostProcessReturn(original_matched_edges,pp_matching,pp_matched_edges,pp_matched_motif,full_rt, profile_results)
        else
            return LocalSearchPostProcessReturn(original_matched_edges,pp_matching,pp_matched_edges,pp_matched_motif,nothing, nothing,)
        end

    else
        throw(ArgumentError("Post processing is only supported for ΛTAME_M got $(method)"))
    end
end

"""-----------------------------------------------------------------------------
  This function takes in a sparse matrix and builds a ThirdOrderSymTensor from
  the triangles present in the networks. If the network isn't symmetric, then
  it's symmetrized. The bool 'use_lcc' is used to indicate whether or not to find
  the largest strongly connected component before symmetrizing.
-----------------------------------------------------------------------------"""
function graph_to_ThirdOrderTensor(A;use_lcc=false)

    n,n = size(A)

    if use_lcc
        A,_ = largest_component(A)
    end

    if !issymmetric(A)
        println("Symmetrising matrix")
		A = max.(A,A')  #symmetrize for Triangles routine
	end

    #build COOTens from graphs
    tris = collect(MatrixNetworks.triangles(A))
    nnz = length(tris)

    indices = Array{Int,2}(undef,3,nnz)
                                # julia uses col-major
                                # formatting.
    for i =1:nnz
        indices[1,i]= tris[i][1]
        indices[2,i]= tris[i][2]
        indices[3,i]= tris[i][3]
    end

    vals = ones(nnz)

    return ThirdOrderSymTensor(n,indices,vals)

end

"""------------------------------------------------------------------------------
  This function aligns graphs using their tensor representations and returns the
  profiled version of the algorithms. These routines call the 'param_search'
  functions for the associated method used.

  Inputs
  ------
  * A, B - (ThirdOrderSymTensor):
	Two third order tensors representing the presence of triangles within the
	network. A must be larger than B, else the routines will be called with the
	parameters swapped.
  * method - (String):
	The choice of method used to align the methods. Options include 'LambdaTAME',
	'LowRankTAME', and 'TAME'.
  * 'no_matching' - (Bool):
	Will not run the matching routines when True. This is useful when studying
	the ranks of the iterates. Any counts which may rely on the matchings are
	replaced by -1.
  Outputs
  -------
  * 'best_TAME_PP_tris' - (Int):
    The largest number of triangles matched over all iterations.
  * 'max_triangle_match' - (Int):
	The maximum number of triangles matchable. This is the minimum between the
	number of triangles in graphs A and B.
  * 'best_TAME_PP_x' - (Array{Float,2}):
	The best iterate found over all the alphas and betas specified by the user.
	when 'LambdaTAME' and 'LowRankTAME' are called, this is replaced by the U
	and V components of the best iterate.
  * 'best_matching' - (Dict{Int,Int}):
	The best matching found, maps from A to B. When the method flips A and B,
	the dictionary is also flipped when returned.
  * profile - (Dict):
	A dictionary storing the profiling results from each of the methods. Please
	see the '_profiled' versions of the code to see what is returned by each
	function.
------------------------------------------------------------------------------"""
function align_tensors_profiled(A::Union{ThirdOrderSymTensor,SymTensorUnweighted{S},Vector{SymTensorUnweighted{S}}},
                                B::Union{ThirdOrderSymTensor,SymTensorUnweighted{S},Vector{SymTensorUnweighted{S}}};
					            method::AlignmentMethod=ΛTAME_M(),no_matching=false,kwargs...) where {S <: Motif}

	if typeof(method) == ΛTAME_M || typeof(method) === ΛTAME_MultiMotif_M
		return ΛTAME_param_search_profiled(A,B;kwargs...)

    elseif typeof(method) == LowRankTAME_MultiMotif_M
        if length(A) > 1
            @warn "multi motifs are not actually supported, only one motif at a time currently"
        end
        return LowRankTAME_param_search_profiled(A[end],B[end];no_matching,kwargs...)
	elseif typeof(method) === LowRankTAME_M
		return LowRankTAME_param_search_profiled(A,B;no_matching,kwargs...)
	elseif typeof(method) === TAME_M
		return TAME_param_search_profiled(A,B;no_matching,kwargs...)
    elseif typeof(method) == TAME_MultiMotif_M
        if length(A) > 1
            @warn "multi motifs are not actually supported, only one motif at a time currently"
        end
        return TAME_param_search_profiled(A[end],B[end];no_matching,kwargs...)
	else
		throw(ArgumentError("method must be one of 'LambdaTAME', ΛTAME_MultiMotif_M, 'LowRankTAME','LowRankTAME_MultiMotif_M','TAME_MultiMotif_M', or 'TAME'."))
	end
end

"""------------------------------------------------------------------------------
  This function aligns graphs using their tensor representations. These routines
  call the param_search functions for the associated method used.

  Inputs
  ------
  * A, B - (ThirdOrderSymTensor,SymTensorUnweighted,Array{SymTensorUnweighted,1}):
	Two third order tensors representing the presence of triangles within the
	network. A must be larger than B, else the routines will be called with the
	parameters swapped.
  * method - (AlignmentMethod):
	The choice of method used to align the methods. Options include LambdaTAME_M,
	LowRankTAME_M, and TAME_M.

  Outputs
  -------
	* 'best_TAME_PP_tris' - (Int):
		The largest number of triangles matched over all iterations.
	* 'max_triangle_match' - (Int):
	  The maximum number of triangles matchable. This is the minimum between the
	  number of triangles in graphs A and B.
	* 'best_TAME_PP_x' - (Array{Float,2}):
	  The best iterate found over all the alphas and betas specified by the user.
	  when 'LambdaTAME' and 'LowRankTAME' are called, this is replaced by the U
	  and V components of the best iterate.

-----------------------------------------------------------------------------"""
function align_tensors(A::Union{ThirdOrderSymTensor,SymTensorUnweighted{S},Vector{SymTensorUnweighted{S}}},
	                   B::Union{ThirdOrderSymTensor,SymTensorUnweighted{S},Vector{SymTensorUnweighted{S}}};
					   method::AlignmentMethod=ΛTAME_M(),no_matching=false,kwargs...) where {S <: Motif}

	if typeof(method) === ΛTAME_M || typeof(method) === ΛTAME_MultiMotif_M
		return ΛTAME_param_search(A,B;kwargs...)
    elseif typeof(method) == LowRankTAME_MultiMotif_M
        if length(A) > 1
            @warn "multi motifs are not actually supported, only one motif at a time currently"
        end
        return LowRankTAME_param_search(A[end],B[end];no_matching,kwargs...)
	elseif typeof(method) === LowRankTAME_M
		return LowRankTAME_param_search(A,B;no_matching,kwargs...)
	elseif typeof(method) === TAME_M
		return TAME_param_search(A,B;no_matching,kwargs...)
    elseif typeof(method) == TAME_MultiMotif_M
        if length(A) > 1
            @warn "multi motifs are not actually supported, only one motif at a time currently"
        end
        return TAME_param_search(A[end],B[end];no_matching,kwargs...)
	else
		throw(ArgumentError("method must be one of LambdaTAME_M, ΛTAME_MultiMotif_M, LowRankTAME_M, or TAME_M."))
	end

end

#=
function align_tensors(A::Array{SymTensorUnweighted{S},1}, B::Array{SymTensorUnweighted{S},1};
			           method::AlignmentMethod=ΛTAME_M(),no_matching=false,kwargs...) where {S <: Motif}

	if typeof(method) === ΛTAME_M || typeof(method) === ΛTAME_MultiMotif_M
		return ΛTAME_param_search(A,B;kwargs...)
	elseif typeof(method) === LowRankTAME_M
		return LowRankTAME_param_search(A,B;no_matching = no_matching,kwargs...)
	elseif typeof(method) === TAME_M
		return TAME_param_search(A,B;no_matching = no_matching,kwargs...)
	else
		throw(ArgumentError("method must be one of LambdaTAME_M,ΛTAME_MultiMotif_M, LowRankTAME_M, or TAME_M."))
	end

end
=#

function align_tensors(graph_A_file::String,graph_B_file::String;
                       ThirdOrderSparse=true,profile=false,
                       kwargs...)
    if haskey(kwargs,:postProcessing)
        kwargs = Dict([(k,v) for (k,v) in kwargs if k != :postProcessing])
    end #TODO: this is a quick fix

    if ThirdOrderSparse
        A = load_ThirdOrderSymTensor(graph_A_file)
        B = load_ThirdOrderSymTensor(graph_B_file)
    else
        A = load(graph_A_file,false,"COOTen")
        B = load(graph_B_file,false,"COOTen")
    end

    if profile
        return align_tensors_profiled(A,B;kwargs...)
    else
        return align_tensors(A,B;kwargs...)
    end
end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
# export readuwevt, uwpf, uwpf!
#
"""
    writesac(W::SeisEvent[, v=0])

Write all data in SeisEvent structure `W` to auto-generated SAC files. Event
header information is written from W.hdr; W.source is not used as there is no
standard header position for event source information.
"""
function writesac(S::SeisEvent; v::Integer=KW.v)
  tdata = Array{Float32,1}(undef, 0)
  reset_sacbuf()

  evid  = codeunits(S.hdr.id == "" ? "-12345  " : S.hdr.id)
  evid  = evid[1:min(length(evid),16)]
  BUF.sac_cv[9:length(evid)+8] .= evid

  # Values from event header
  S.hdr.loc.lat == 0.0 || setindex!(BUF.sac_fv, Float32(S.hdr.loc.lat), 40)
  S.hdr.loc.lon == 0.0 || setindex!(BUF.sac_fv, Float32(S.hdr.loc.lon), 41)
  S.hdr.loc.dep == 0.0 || setindex!(BUF.sac_fv, Float32(S.hdr.loc.dep), 42)
  S.hdr.mag.val == -5.0f0 || setindex!(BUF.sac_fv, Float32(S.hdr.mag.val), 44)
  BUF.sac_cv[9:length(evid)+8] .= evid
  t_evt = d2u(S.hdr.ot)

  # Ints, always
  BUF.sac_iv[7] = Int32(6)

  for i = 1:S.data.n
    BUF.sac_fv[8] = t_evt - S.data.t[i][1,2]*μs
    write_sac_channel(S.data, i, "", v)
  end
  return nothing
end


text:
@testset "Elastic net ($T)" for T in [Float32, Float64, ComplexF32, ComplexF64]
    using ProximalOperators
    using ProximalAlgorithms
    using LinearAlgebra

    A = T[
        1.0 -2.0 3.0 -4.0 5.0
        2.0 -1.0 0.0 -1.0 3.0
        -1.0 0.0 4.0 -3.0 2.0
        -1.0 -1.0 -1.0 1.0 3.0
    ]
    b = T[1.0, 2.0, 3.0, 4.0]

    m, n = size(A)

    R = real(T)

    reg = ElasticNet(R(1), R(1))
    reg1 = NormL1(R(1))
    reg2 = SqrNormL2(R(1))
    loss = Translate(SqrNormL2(R(1)), -b)
    cost = LeastSquares(A, b)

    L = opnorm(A)^2

    x_star = T[-0.6004983388704322, 0.0, 0.0, 0.195182724252491, 0.764119601328903]

    @testset "DavisYin" begin

        # with known initial iterate

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        solver = ProximalAlgorithms.DavisYin(tol = R(1e-6))
        xf_dys, xg_dys, it_dys = solver(x0, f = reg1, g = reg2, h = cost, Lh = L)
        @test eltype(xf_dys) == T
        @test eltype(xg_dys) == T
        @test norm(xf_dys - x_star, Inf) <= 1e-3
        @test norm(xg_dys - x_star, Inf) <= 1e-3
        @test it_dys <= 140
        @test x0 == x0_backup

        # with random initial iterate

        x0 = randn(T, n)
        x0_backup = copy(x0)
        solver = ProximalAlgorithms.DavisYin(tol = R(1e-6))
        xf_dys, xg_dys, it_dys = solver(x0, f = reg1, g = reg2, h = cost, Lh = L)
        @test eltype(xf_dys) == T
        @test eltype(xg_dys) == T
        @test norm(xf_dys - x_star, Inf) <= 1e-3
        @test norm(xg_dys - x_star, Inf) <= 1e-3
        @test x0 == x0_backup

    end

    afba_test_params = [
        (R(2), R(0), 130),
        (R(1), R(1), 1890),
        (R(0), R(1), 320),
        (R(0), R(0), 194),
        (R(1), R(0), 130),
    ]

    @testset "AFBA" for (theta, mu, maxit) in afba_test_params

        # with known initial iterates

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        y0 = zeros(T, m)
        y0_backup = copy(y0)

        solver = ProximalAlgorithms.AFBA(theta = theta, mu = mu, tol = R(1e-6))
        x_afba, y_afba, it_afba =
            solver(x0, y0, f = reg2, g = reg1, h = loss, L = A, beta_f = R(1))
        @test eltype(x_afba) == T
        @test eltype(y_afba) == T
        @test norm(x_afba - x_star, Inf) <= 1e-4
        @test it_afba <= maxit
        @test x0 == x0_backup
        @test y0 == y0_backup

        # with random initial iterates

        x0 = randn(T, n)
        x0_backup = copy(x0)
        y0 = randn(T, m)
        y0_backup = copy(y0)

        solver = ProximalAlgorithms.AFBA(theta = theta, mu = mu, tol = R(1e-6))
        x_afba, y_afba, it_afba =
            solver(x0, y0, f = reg2, g = reg1, h = loss, L = A, beta_f = R(1))
        @test eltype(x_afba) == T
        @test eltype(y_afba) == T
        @test norm(x_afba - x_star, Inf) <= 1e-4
        @test x0 == x0_backup
        @test y0 == y0_backup

    end
end


text:
struct DESolutionIterator{T,S}
    sol::S
end

IteratorInterfaceExtensions.isiterable(sol::DESolution) = true
TableTraits.isiterabletable(sol::DESolution) = true

function IteratorInterfaceExtensions.getiterator(sol::DESolution)
    timestamp_type = eltype(sol.t)
    value_type = eltype(sol.u)

    value_types = Type[]
    value_names = Symbol[]

    push!(value_types, timestamp_type)
    push!(value_names, :timestamp)

    if value_type<:AbstractArray
        for i in 1:length(sol.u[1])
            push!(value_types, eltype(sol.u[1]))
            if has_syms(sol.prob.f)
                push!(value_names, sol.prob.f.syms[i])
            else
                push!(value_names, Symbol("value$i"))
            end
          end
    else
        push!(value_types, value_type)
        if has_syms(sol.prob.f)
            push!(value_names, sol.prob.f.syms[1])
        else
            push!(value_names, :value)
        end
    end

    si = DESolutionIterator{NamedTuple{(value_names...,),Tuple{value_types...}},typeof(sol)}(sol)

    return si
end

function Base.length(iter::DESolutionIterator)
    return length(iter.sol)
end

Base.eltype(::Type{DESolutionIterator{T,TS}}) where {T,TS} = T

@generated function Base.iterate(iter::DESolutionIterator{T,S}, state=1) where {T,S}
    columns = []
    push!(columns, :(iter.sol.t[state]))
    nfield = length(fieldnames(T))
    for i in 1:nfield-1
        if nfield > 2
            push!(columns, :(iter.sol.u[state][$i]))
        else
            push!(columns, :(iter.sol.u[state]))
        end
    end

    quote
        if state > length(iter)
            return nothing
        else
            return $(T)(($(columns...),)), state+1
        end
    end
end


text:
module RunTests
import Example

Example.domath(0)

end


text:
# v0.12 deprecations

function ones(dims...)
  Base.depwarn("Flux.ones(size...) is deprecated, please use Flux.ones32(size...) or Base.ones(Float32, size...)", :ones, force=true)
  Base.ones(Float32, dims...)
end
ones(T::Type, dims...) = Base.ones(T, dims...)

function zeros(dims...)
  Base.depwarn("Flux.zeros(size...) is deprecated, please use Flux.zeros32(size...) or Base.zeros(Float32, size...)", :zeros, force=true)
  Base.zeros(Float32, dims...)
end
zeros(T::Type, dims...) = Base.zeros(T, dims...)

ones32(::Type, dims...) = throw(ArgumentError("Flux.ones32 is always Float32, use Base.ones to specify the element type"))
zeros32(::Type, dims...) = throw(ArgumentError("Flux.zeros32 is always Float32, use Base.zeros to specify the element type"))

# v0.13 deprecations

function Broadcast.broadcasted(f::Recur, args...)
  # This had an explicit @adjoint rule, calling Zygote.∇map(__context__, f, args...), until v0.12
  Base.depwarn("""Broadcasting is not safe to use with RNNs, as it does not guarantee an iteration order.
    Re-writing this as a comprehension would be better.""", :broadcasted)
  map(f, args...)  # map isn't really safe either, but
end

@deprecate frequencies(xs) group_counts(xs)

struct Zeros
  function Zeros()
    Base.depwarn("Flux.Zeros is no more, has ceased to be, is bereft of life, is an ex-boondoggle... please use bias=false instead", :Zeros)
    false
  end
end
Zeros(args...) = Zeros()  # was used both Dense(10, 2, initb = Zeros) and Dense(rand(2,10), Zeros())

function Optimise.update!(x::AbstractArray, x̄)
  depwarn("`Flux.Optimise.update!(x, x̄)` was not used internally and has been removed. Please write `x .-= x̄` instead.", :update!)
  x .-= x̄
end

function Diagonal(size::Integer...; kw...)
  Base.depwarn("Flux.Diagonal is now Flux.Scale, and also allows an activation function.", :Diagonal)
  Scale(size...; kw...)
end
function Diagonal(size::Tuple; kw...)
  Base.depwarn("Flux.Diagonal is now Flux.Scale, and also allows an activation function.", :Diagonal)
  Scale(size...; kw...)
end

# Deprecate this eventually once saving models w/o structure is no more
function loadparams!(m, xs)
  Base.depwarn("loadparams! will be deprecated eventually. Use loadmodel! instead.", :loadparams!)
  for (p, x) in zip(params(m), xs)
    size(p) == size(x) ||
      error("Expected param size $(size(p)), got $(size(x))")
    copyto!(p, x)
  end
end

# Channel notation: Changed to match Conv, but very softly deprecated!
# Perhaps change to @deprecate for v0.14, but there is no plan to remove these.
Dense(in::Integer, out::Integer, σ = identity; kw...) =
  Dense(in => out, σ; kw...)
Bilinear(in1::Integer, in2::Integer, out::Integer, σ = identity; kw...) =
  Bilinear((in1, in2) => out, σ; kw...)
Embedding(in::Integer, out::Integer; kw...) = Embedding(in => out; kw...)

RNNCell(in::Integer, out::Integer, σ = tanh; kw...) = RNNCell(in => out, σ; kw...)
LSTMCell(in::Integer, out::Integer; kw...) = LSTMCell(in => out; kw...)

GRUCell(in::Integer, out::Integer; kw...) = GRUCell(in => out; kw...)
GRUv3Cell(in::Integer, out::Integer; kw...) = GRUv3Cell(in => out; kw...)


text: const Ctime_t = Int


text:
#This module handles importing local data from the eve client
# as well as parsing logs

module LocalData
using eve.SDE
using eve.MarketDatas
using DataFrames
# Log import
export latestLog
export latestOrderLog
export latestMarketLog
export importOrders
export importMarketLog

function logPath()
  joinpath(homedir(), "Documents", "EVE", "logs", "Marketlogs")
end
function latestOrderLog(typePrefix = "Corporation")
  fullPrefix = typePrefix * " Orders"
  reg = Regex("$fullPrefix-(.+)\\.txt")
  local files :: Array{Tuple{ASCIIString, DateTime}} = []
  df = Dates.DateFormat("yyyy.mm.dd HHMM")
  for file in readdir(logPath())
    m = match(reg, file)
    if m != nothing
      push!(files, (joinpath(logPath(), file), DateTime(m.captures[1], df)))
    end
  end
  sort!(files, by=x->x[2], rev=true)
  if isempty(files)
    return Nullable{Tuple{ASCIIString,DateTime}}()
  else
    return Nullable(files[1])
  end
end
immutable MarketLog
  path :: ASCIIString
  timestamp :: DateTime
  region :: ASCIIString
  item :: ASCIIString
end
function latestMarketLog(item = "", region = "")
  exclude = Regex("Corporation Orders|My Orders")
  groups = Regex("(.+)-(.+)-(.+)\\.txt" )
  local files :: Array{MarketLog} = []
  df = Dates.DateFormat("yyyy.mm.dd HHMMSS")
  for file in readdir(logPath())
    m = match(exclude, file)
    info = match(groups, file)
    if m == nothing && contains(info.captures[1], region) && contains(info.captures[2], item)
      push!(files, MarketLog(joinpath(logPath(), file), DateTime(info.captures[3], df), info.captures[1], info.captures[2]))
    end
  end
  sort!(files, by=x->x.timestamp, rev=true)
  if isempty(files)
    return Nullable{MarketLog}()
  else
    return Nullable(files[1])
  end
end
function importOrders(typePrefix = "Corporation")
  orders = latestOrderLog()
  if isnull(orders)
    return DataFrame()
  end
  path = get(latestOrderLog())[1]

  readtable(path, truestrings=["True"], falsestrings=["False"])[1:23]
end
function importMarketLog(item = "", region = "")
  o = get(latestMarketLog(item, region))
  res = readtable(o.path, truestrings=["True"], falsestrings=["False"])[1:14]
  names!(res, [:Price,
               :Quantity,
               :TypeID,
               :Range,
               :OrderID,
               :QuantityEntered,
               :MinQuantity,
               :bid,
               :Issued,
               :Duration,
               :StationID,
               :RegionID,
               :solarSystemID,
               :Jumps])
  (sells, buys) = groupby(res, :bid)
  MarketData(itemID(o.item),
             Logs,
             o.timestamp,
             sells,
             buys)
end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    @static

Partially evaluate an expression at parse time.

For example, `@static Sys.iswindows() ? foo : bar` will evaluate `Sys.iswindows()` and insert
either `foo` or `bar` into the expression.
This is useful in cases where a construct would be invalid on other platforms,
such as a `ccall` to a non-existent function.
`@static if Sys.isapple() foo end` and `@static foo <&&,||> bar` are also valid syntax.
"""
macro static(ex)
    if isa(ex, Expr)
        if ex.head === :if || ex.head === :&& || ex.head === :||
            cond = eval(__module__, ex.args[1])
            if xor(cond, ex.head === :||)
                return esc(ex.args[2])
            elseif length(ex.args) == 3
                return esc(ex.args[3])
            elseif ex.head === :if
                return nothing
            else
                return cond
            end
        end
    end
    throw(ArgumentError("invalid @static macro"))
end

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
################################################################################
########## Plotting functions for Network using Gadfly                ##########
########## Maxime Legrand - Berkeley 2015                             ##########
################################################################################

# Pkg.add("Gadfly")
# using Gadfly
# include("../MCL_Julia/mcl.jl")

"Plots n vertices given through a n*2 array."
function plot_vertices(Coordinates)
  x_coor = Coordinates[:,1]
  y_coor = Coordinates[:,2]
  plot(x=x_coor, y=y_coor, Geom.point)
end

"""Graph Plot
Input : V, n*2 matrix of vextices coordinates
    A, n*n adjacency matrix
    C, a vector of clusters to differentiate by vertices' color, Vector[] by default
Output : Corresponding undirected n-network plot"""
function graphplot(V, A; C=Vector[], title="", labels=false, key_position=:none)
  n = size(V,1)
  xn1 = Float64[]
  yn1 = Float64[]
  xn2 = Float64[]
  yn2 = Float64[]
  for i in 1:n
    for j in 1:i
      if A[i,j] != 0
        push!(xn1, V[i,1])
        push!(yn1, V[i,2])
        push!(xn2, V[j,1])
        push!(yn2, V[j,2])
      end
    end
  end
  layers = Layer[]
  if length(C) == 0
        push!(layers, layer(x = V[:,1], y = V[:,2], Geom.point)[1])
  else
    part = map(k -> "$(vvfind(k,C))", 1:n)
        push!(layers, layer(x = V[:,1], y = V[:,2], color = part, Geom.point)[1])
  end
  for i=1:length(xn1)
    push!(layers,
          layer(x = [xn1[i], xn2[i]], y = [yn1[i], yn2[i]], Geom.line, Theme(default_color=color("lightgrey")))[1])
  end
  if labels
    return plot(x = V[:,1], y = V[:,2], label = map(x->"$x",[1:n]), Geom.label,
        layers, Guide.xlabel(""), Guide.ylabel(""), Guide.title(title),
        Theme(grid_line_width=0mm, minor_label_font_size=0mm, key_position=key_position))
  else
    return plot(layers, Guide.xlabel(""), Guide.ylabel(""), Guide.title(title),
        Theme(grid_line_width=0mm, minor_label_font_size=0mm, key_position=key_position))
  end
end

"Regular polygon positionning"
function rpolygon(n, center, radius)
  nodes = zeros(n,2)
  for i in 1:n
    nodes[i,1] = center[1] + radius * cos(2*i*pi/n)
    nodes[i,2] = center[2] + radius * sin(2*i*pi/n)
  end
  return nodes
end

"Superior part"
function spart(x::Float64)
  if x >= 0.
    return x
  else
    return 0.
  end
end

"""Force-directed layout
Inspired from physics, but not physically correct
Notations of "Simple Algorithms for Network Visualization" by M. J. McGuffin
Input : A, adjacency matrix of an undirected network
    K, number of steps taken by the algorithm
    C, clusters to color
    init, string setting the original position of the vertices
    nodes, if init is set to "custom", initial position of the nodes
    key_position, among {:right, :left, :top, :bottom, :none}
Output : corresponding graph"""
function netplot(A; K=100, C=Vector[], init="rand", nodes=[], title="", labels=false, key_position=:none)
  n = size(A,1)
  R = 10
  L = 2*R/(n/2)
  Kr = 1
  Ks = 1
  delta = 2*R/K
  d2max = 2*R
  if init == "polygon"
    nodes = convert(Array{Float64,2}, rpolygon(n, [0 0], R))
  elseif init == "rand"
    nodes = -R + 2*R*rand(n,2)
  end
  for k in 1:K
    forces = zeros(n,2)
    for i1 in 1:(n-1)
      for i2 in (i1+1):n
        dx = nodes[i2,1] - nodes[i1,1]
        dy = nodes[i2,2] - nodes[i1,2]
        if (dx != 0) || (dy != 0)
          d2 = dx*dx + dy*dy
          d = sqrt(d2)
          # Repulsive force Fr = Kr/d^2
          Fr = Kr / ((L/10) + d2)
          fx = Fr * dx / d
          fy = Fr * dy / d
          forces[i1,1] = forces[i1,1] - fx
          forces[i1,2] = forces[i1,2] - fy
          forces[i2,1] = forces[i2,1] + fx
          forces[i2,2] = forces[i2,2] + fy
          # Spring force Fs = Ks(d-L)
          if A[i1,i2] != 0
            Fs = Ks * spart(d - L)
            fx = Fs * dx / d
            fy = Fs * dy / d
            forces[i1,1] = forces[i1,1] + fx
            forces[i1,2] = forces[i1,2] + fy
            forces[i2,1] = forces[i2,1] - fx
            forces[i2,2] = forces[i2,2] - fy
          end
        end
      end
    end
    # Updating positions
    for i in 1:n
      dx = delta * forces[i,1]
      dy = delta * forces[i,2]
      nodes[i,1] = nodes[i,1] + dx
      nodes[i,2] = nodes[i,2] + dy
    end
  end
    return graphplot(nodes, A, C=C, title=title, labels=labels, key_position=key_position)
end


text:
const PTS_VERTICAL_LINE =
[
    20 20;
    20 25;
    20 30;
    20 35;
    20 40;
    20 45.0
]

const PTS_ELLIPTIC =
[       3 3;
    2.759 3.363;
    2.372 3.663;
    1.890 3.862;
    1.355 3.943
]

const PTS_PAIR_1 =
[
    1.768 2.3311;
    1.947 2.6271;
    1.595 2.7951;
    1.019 2.7241;
    0.479 2.4281;
    0.126 2.0521;
   -0.001 1.720;
    0.103 1.514;
    0.442 1.549;
    1.055 1.905
]

const PTS_PAIR_2 =
[
   1.9592 2.44973;
    2.168 2.675;
    1.821 2.804;
    1.244 2.720;
    0.705 2.437;
    0.346 2.104;
    0.195 1.833;
    0.356 1.680;
    0.558 1.742;
    1.186 2.088
]


text:
abstract type AbstractExaminee end
include("latent/latent.jl")

"""
    Examinee <: AbstractExaminee

# Description

An immutable containing the information about the examinee (the test taker).

# Fields

  - **`idx::Int64`**: An integer that identifies the Examinee in this session.
  - **`id::String`**: A string that identifies the Examinee.
  - **`latent::Latent`**: A mutable latent struct associated with the Examinee.

# Factories
    Examinee(idx::Int64, id::String, latent::AbstractLatent) = new(idx, id, latent)
    Examinee(idx::Int64, latent::AbstractLatent) = new(idx, string(idx), latent)

Creates a new Examinee with custom index, (id) and a latent variable.

# Random default initializers (1-D latent)
    Examinee(idx::Int64) = new(idx, string(idx), Latent1D())
    Examinee(idx::Int64, id::String) = new(idx, id, Latent1D())

# Description

Randomly generates an Examinee with custom index and id and with a default 1-dimensional latent variable
(Look at (`Latent1D`)[#Psychometrics.Latent1D] for the defaults).
"""
struct Examinee <: AbstractExaminee
    idx::Int64
    id::String
    latent::AbstractLatent

    # Factories
    Examinee(idx::Int64, id::String, latent::AbstractLatent) = new(idx, id, latent)
    Examinee(idx::Int64, latent::AbstractLatent) = new(idx, string(idx), latent)

    # Random default initializers
    Examinee(idx::Int64) = new(idx, string(idx), Latent1D())
    Examinee(idx::Int64, id::String) = new(idx, id, Latent1D())
end

function Base.copy(examinee::Examinee)
    e = Examinee(examinee.idx, examinee.id, examinee.latent)
    return e::Examinee
end

"""
    get_examinee_by_id(examinee_id::String, examinees::Vector{<:AbstractExaminee})

It returns the Examinee with index `examinee_id` from a Vector of AbstractExaminee.
"""
function get_examinee_by_id(examinee_id::String, examinees::Vector{<:AbstractExaminee})
    filter(e -> e.id == examinee_id, examinees)[1]
end

"""
    empty_chain!(examinee::AbstractExaminee)
"""
function empty_chain!(examinee::AbstractExaminee)
    _empty_chain!(examinee.latent)
end

"""
    set_val!(examinee::AbstractExaminee, val::Float64)
"""
function set_val!(examinee::AbstractExaminee, val::Float64)
    _set_val!(examinee.latent, val)
end

"""
    set_val_from_chain!(examinee::AbstractExaminee)
"""
function set_val_from_chain!(examinee::AbstractExaminee)
    _set_val_from_chain!(examinee.latent)
end

"""
    update_estimate!(examinee::AbstractExaminee; sampling = true)
"""
function update_estimate!(examinee::AbstractExaminee; sampling = true)
    _update_estimate!(examinee.latent; sampling = sampling)
end

"""
    chain_append!(examinee::AbstractExaminee; sampling = false)
"""
function chain_append!(examinee::AbstractExaminee; sampling = false)
    _chain_append!(examinee.latent; sampling = sampling)
end

"""
    set_prior!(
        examinee::AbstractExaminee,
        prior::Union{Distributions.DiscreteUnivariateDistribution, Distributions.ContinuousUnivariateDistribution}
    )
"""
function  set_prior!(
    examinee::AbstractExaminee,
    prior::Union{Distributions.DiscreteUnivariateDistribution, Distributions.ContinuousUnivariateDistribution}
)
    _set_prior!(examinee.latent, prior)
    return nothing
end

function get_latents(examinee::AbstractExaminee)
    return examinee.latent::AbstractLatent
end

function get_latents_vals(examinee::AbstractExaminee)
    return examinee.latent.val
end


text:
# ---
# title: JuliaRL\_VPG\_CartPole
# cover: assets/JuliaRL_VPG_CartPole.png
# description: VPG applied to CartPole
# date: 2021-05-22
# author: "[norci](https://github.com/norci)"
# ---

#+ tangle=true
using ReinforcementLearning
using StableRNGs
using Flux
using Flux.Losses
using Distributions

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:VPG},
    ::Val{:CartPole},
    ::Nothing;
    seed = 123,
)
    rng = StableRNG(seed)
    env = CartPoleEnv(; T = Float32, rng = rng)
    ns, na = length(state(env)), length(action_space(env))

    agent = Agent(
        policy = VPGPolicy(
            approximator = NeuralNetworkApproximator(
                model = Chain(
                    Dense(ns, 128, relu; init = glorot_uniform(rng)),
                    Dense(128, 128, relu; init = glorot_uniform(rng)),
                    Dense(128, na; init = glorot_uniform(rng)),
                ),
                optimizer = ADAM(),
            ) |> cpu,
            baseline = NeuralNetworkApproximator(
                model = Chain(
                    Dense(ns, 128, relu; init = glorot_uniform(rng)),
                    Dense(128, 128, relu; init = glorot_uniform(rng)),
                    Dense(128, 1; init = glorot_uniform(rng)),
                ),
                optimizer = ADAM(),
            ) |> cpu,
            action_space = action_space(env),
            dist = Categorical,
            γ = 0.99f0,
            rng = rng,
        ),
        trajectory = ElasticSARTTrajectory(state = Vector{Float32} => (ns,)),
    )
    stop_condition = StopAfterEpisode(500, is_show_progress=false)

    hook = TotalRewardPerEpisode()
    description = "# Play CartPole with VPG"
    Experiment(agent, env, stop_condition, hook, description)
end

#+ tangle=false
using Plots
pyplot() #hide
ex = E`JuliaRL_VPG_CartPole`
run(ex)
plot(ex.hook.rewards)
savefig("assets/JuliaRL_VPG_CartPole.png") #hide

# ![](assets/JuliaRL_VPG_CartPole.png)


text:
# perfctr_gpu_convenience.jl
using LIKWID
using LinearAlgebra
using CUDA

@assert CUDA.functional()

GPUMarker.init()

# Note: CUDA defaults to Float32
Agpu = CUDA.rand(128, 64)
Bgpu = CUDA.rand(64, 128)
Cgpu = CUDA.zeros(128, 128)

@gpuregion "matmul" for _ in 1:100
    mul!(Cgpu, Agpu, Bgpu)
end

GPUMarker.close()


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

using Markdown
using InteractiveUtils

# ╔═╡ d16725a0-6cd3-11ec-0998-fd4af26f10dd
begin
	using Pkg; Pkg.activate(@__DIR__); Pkg.instantiate()
	Pkg.precompile()

	using PlutoUI
	using Plots
	using Statistics
	using StatsBase
	using PyCall
	using DataFrames
	using GLM
	using Tables
	using XLSX
	using DataStructures
	using MLBase
	using RDatasets
	using LsqFit
end;

# ╔═╡ 91471950-703f-4927-8bce-f4478f4a53ba
PlutoUI.TableOfContents(aside=true, indent=true, depth=3)

# ╔═╡ a53d44c5-75c4-4508-aa19-33e838c9e0bb
md"""
# 📈 Regression

Regression Analysis is a set of statistical processes for estimating the relationships between a _dependent variable_ (often called the _outcome_ or _response_ variable) and one or more _independent variables_ (often called _predictors_, _covariates_, _explanatory variables_ or _features_). The most common form of regression analysis is linear regression, in which one finds the line (or a more complex linear combination) that most closely fits the data according to a specific mathematical criterion.

[YouTube Link](https://www.youtube.com/watch?v=5TCbIK_cpZE)
"""

# ╔═╡ e9458f99-4ec0-4167-b771-3090dfec6320
md"""
## Simple linear regression
"""

# ╔═╡ 784c5b0e-59dc-4b2d-bf40-112261ac8e85
md"""
### Building from scratch

Firstly, let's generate some random data and then plot it...
"""

# ╔═╡ 65e22edc-ef5c-412f-82f4-1c84da209d6d
begin
	# x values
	xvals = repeat(1:0.5:10, inner=2)
	# y values (with some noise)
	yvals = 3 .+ xvals .+ 2 .* rand(length(xvals)) .-1

	# plotting
	scatter(xvals, yvals, color=:black, leg=false,
			xlabel="X values", ylabel="Y values")
end

# ╔═╡ db2c3459-2127-4b46-8e73-546850767dc0
md"""
Assuming an ordinary least squares regression, the slope `a` of the fitted line is equal to the correlation between `y` and `x` corrected by the ratio of standard deviations of these variables. The intercept `b` of the fitted line is such that the line passes through the center of mass (`x`, `y`) of the data points.

Thus, we can easily build the `find_best_fit` function, which returns `a` and `b` coefficients.
"""

# ╔═╡ 3bb8c382-f97d-48dc-8cbe-99bf77de0784
function find_best_fit(x, y)
	X̅, Y̅ = mean(x), mean(y)
	Sˣ, Sʸ = std(x), std(y)

	a = cor(x,y) * (Sˣ/Sʸ)
	b = Y̅ - a*X̅

	return a, b
end

# ╔═╡ 75059ff4-6b9d-4e01-9aed-33ea8b0e9a71
md"""
We can now use the function above to calculate `a` and `b` and then compute the estimated values `ŷvals`...
"""

# ╔═╡ e83331a1-8655-4de0-a36f-b74b5e43674b
begin
	a, b = find_best_fit(xvals, yvals)
	ŷvals₁ = a .* xvals .+ b
end

# ╔═╡ d81e448b-b4e4-4d6f-8ca7-42c41569977d
md"""
### numpy

We can also import the `numpy` python package to perform a simple linear regression. We are able to use any `numpy` function by adding the prefix `np`, just like in Python.
"""

# ╔═╡ a047abfb-178c-40ac-8f8c-c3a8a73c1849
# importing numpy
np = pyimport("numpy");

# ╔═╡ 208af138-5131-4034-885b-c41f6ccf36f5
# numpy model
np_model = np.polyfit(xvals, yvals, 1)

# ╔═╡ 39884fc4-c41c-46e5-a92a-44d766243715
# estimated ŷ values
ŷvals₂ = np_model[1] .* collect(xvals) .+ np_model[2]

# ╔═╡ 13719feb-60ae-433b-bf1e-df6ac045215d
md"""
### GLM

Finally, we can use the [GLM.jl](https://github.com/JuliaStats/GLM.jl) package to perform simple linear regression using the `lm` function.
"""

# ╔═╡ 10c072d4-b793-458c-b947-d5dc3aae5e25
# converting data into a DataFrame
data = DataFrame(X=xvals, Y=yvals);

# ╔═╡ 68a0574f-5b17-48ae-b8c4-a03c7e385561
# GLM model
glm_model = lm(@formula(Y ~ X), data)

# ╔═╡ ecb1cb29-2b57-4197-9327-e762fc876281
# estimated ŷ values
ŷvals₃ = predict(glm_model)

# ╔═╡ 056b28f3-8da4-4a94-a33c-72b30c8aa0aa
md"""
### Comparing methods

Now we can plot all the three models. We can see below that both `numpy` and `GLM` models are exactly the same. The model that we build from scratch is slightly different.
"""

# ╔═╡ 9a5dcad1-9365-423c-b27a-a6c206850848
begin
	# plotting data
	scatter(xvals, yvals, color=:black, label="Data",
			xlabel="X values", ylabel="Y values", legend=:topleft)

	# plotting "from scratch" fitted line
	plot!(xvals, ŷvals₁, color=:red, label="From Scratch")

	# plotting numpy model
	plot!(xvals, ŷvals₂, color=:green, label="Numpy")

	# plotting GLM model
	plot!(xvals, ŷvals₃, color=:blue, label="GLM")
end

# ╔═╡ a89ec73a-5e35-49ad-8829-dc2ecc7e8aac
md"""
## Load real data

Now let's get some real data. We will use housing information from zillow, check out the file `zillow_data_download_april2020.xlsx` for a quick look of what the data looks like. Our goal will be to build a linear regression model between the number of houses listed vs the number of houses sold in a few states. Fitting these models can serve as a key real estate indicator.
"""

# ╔═╡ aa78c706-8860-46f1-9982-ea10bdc4fccb
# loading data
S = XLSX.readxlsx("data/zillow_data_download_april2020.xlsx")

# ╔═╡ 618ab9ee-ee9c-4dbb-b48f-ebeeb15549c6
md"""
**Note:** there are three different sheets in our data.

## Data wrangling

Let's create a `DataFrame` from the `Sale_counts_city` sheet...
"""

# ╔═╡ b530e842-bfa4-4c41-9cc9-4aa629e12f99
# selecting sheet
sale_counts = S["Sale_counts_city"][:]

# ╔═╡ 041e87c7-deee-46c1-85be-ea47df16217c
# converting Matrix into DataFrame
df_sale_counts = DataFrame(sale_counts[2:end,:], Symbol.(sale_counts[1,:]))

# ╔═╡ 19a44596-2367-4964-b21e-7d84ba147915
md"""
Now, we can do the same with the `MonthlyListings_City` sheet...
"""

# ╔═╡ ab8cf8d2-9ac6-4361-b168-243f781e7dfa
begin
	# selecting sheet
	month_list = S["MonthlyListings_City"][:]
	# converting Matrix into DataFrame
	df_month_list = DataFrame(month_list[2:end,:], Symbol.(month_list[1,:]))
end

# ╔═╡ 215483a1-6059-437d-a4b6-ca3cd897f13b
md"""
We may filter just the data from February 2020...
"""

# ╔═╡ b075d090-91aa-4d1f-bb9f-a93dd65413e7
begin
	# monthly listing data (Feb 2020)
	monthlist_2020_2 = df_month_list[!,[1,2,3,4,5,end]]
	rename!(monthlist_2020_2, Symbol("2020-02") => :listing)
end

# ╔═╡ 1b8f7b0d-a843-4af7-90b7-4a64d74edab0
begin
	# sale counts data (Feb 2020)
	salecounts_2020_2 = df_sale_counts[!,[1,end]]
	rename!(salecounts_2020_2, Symbol("2020-02") => :sales)
end

# ╔═╡ 39456bea-24dc-482b-9d8b-33a9f2ed715a
md"""
Next, we can (outer) join our two DataFrames and, then, remove missing values...
"""

# ╔═╡ f195db71-f2a4-4f5c-bfa4-2d9ec1f23a32
begin
	feb2020data = outerjoin(monthlist_2020_2, salecounts_2020_2, on=:RegionID)
	dropmissing!(feb2020data)
end

# ╔═╡ 2da58c4e-8db4-4d80-aa29-f3ec10952f73
md"""
Let's define variables for all the three columns we are interested in...
"""

# ╔═╡ 93179d9a-5bf4-488b-b710-28a5f56b2e26
begin
	sales = feb2020data[!,:sales]
	counts = feb2020data[!,:listing]
	states = feb2020data[!,:StateName]
end;

# ╔═╡ 560a134f-1175-4d56-83af-1ad6bbe1a99f
md"""
Now, we can select just the top 10 states (most cited states). We will use the `counter` function from the [DataStructures.jl](https://github.com/JuliaCollections/DataStructures.jl) package.
"""

# ╔═╡ 0c4a00f6-69f2-4fc3-a7f8-ca6e38404ac9
# counter object (Accumulator)
C = counter(states)

# ╔═╡ a02a9fdf-46bc-4ea1-bb29-b1731a66c057
# converting into a Dictionary
Cdict = C.map

# ╔═╡ 63b596be-9392-4efb-8d31-7f1d6e9c7254
begin
	# getting counts as an Array
	countvals = values(Cdict) |> collect

	# getting state names as an Array
	state_labels = keys(Cdict) |> collect
end;

# ╔═╡ feced955-2e69-44ae-bc2a-82d4f2136b38
# getting indexes from 10 top states
topstateindexes = sortperm(countvals, rev=true)[1:10]

# ╔═╡ 831484f5-6cdc-4641-8adf-117a33c518c2
md"""
**Note:** the built-in `sortperm` function returns a permutation vector `I` that puts `v[I]` in sorted order. The order is specified using the same keywords as `sort!`.
"""

# ╔═╡ b94019d8-a7d8-4781-819b-28eeac823812
# getting top 10 states names
topstates = state_labels[topstateindexes]

# ╔═╡ 4735fc87-68cb-4617-aa7c-6139ba7b7c0a
md"""
## Linear regression models

Firstly, we will plot the model for each state using the formula `Y ~ X`.
"""

# ╔═╡ ec7ae9e6-b191-4578-8a47-9d2245280316
begin
	# creating a blank Array of Plot
	all_plots₁ = Array{Plots.Plot,1}(undef,10)

	for (i,s) in enumerate(topstates)
		# finding indexes of all occurences of a given state
		ids = findall(feb2020data[!,:StateName] .== s)
		# getting data for a given state
		data = DataFrame(X=Float64.(counts[ids]), Y=Float64.(sales[ids]))
		# training the model
		ols = GLM.lm(@formula(Y ~ X), data)

		# plotting data by each state
		all_plots₁[i] = scatter(data.X, data.Y, markersize=2,
								xlims=(0,500), ylims=(0,500),
								color=i, aspect_ratio=:equal,
								legend=false, title=s)

		# plotting model by each state
		plot!(data.X, predict(ols), color=:black)
	end

	# display plots
	plot(all_plots₁..., layout=(2,5), size=(850,300))
end

# ╔═╡ 216dc17e-54c1-4520-890d-8663165b10dd
md"""
We can do the same thing, but now we will plot all the elements within the same figure.

**Note:** when we use the `Y ~ X + 0` formula, it forces the intercept to be zero.
"""

# ╔═╡ dddfb8fb-7047-4417-8dd2-569a93277f0c
begin
	plot()

	for (i,s) in enumerate(topstates)
	    ids = findall(feb2020data[!,:StateName] .== s)
	    data = DataFrame(X=Float64.(counts[ids]), Y=Float64.(sales[ids]))
	    ols = GLM.lm(@formula(Y ~ 0 + X), data)

	    scatter!(counts[ids], sales[ids], markersize=2, xlim=(0,500),
				 ylim=(0,500), color=i, aspect_ratio=:equal,
	        	 legend=false, marker=(3,3,stroke(0)), alpha=0.2)

	    if s ∈ ["NC","CA","FL"]
	        annotate!([(500-20,10+coef(ols)[1]*500, text(s,10))])
	    end

	    plot!(counts[ids], predict(ols), color=i, linewidth=2)
	end

	xlabel!("Listings")
	ylabel!("Sales")
end

# ╔═╡ 5116c117-75e4-4989-b06f-2883bba48330
md"""
## Logistic regression

So far, we have shown several ways to solve the linear regression problem in Julia. Here, we will first start with a motivating example of when you would want to use logistic regression. Let's assume that our predictor vector is binary (`0` or `1`), let's fit a linear regression model.
"""

# ╔═╡ 3a86069b-3f24-4eec-a7ae-f77189a4fd77
# generating random data
𝒮 = DataFrame(X=collect(1:7), Y=[1,0,1,1,1,1,1])

# ╔═╡ 453385dd-8bda-42b8-9b4a-6a4963b21346
# linear regression model
lin_reg = lm(@formula(Y ~ X), 𝒮)

# ╔═╡ bd31693d-dee9-4534-b4d4-3904997e8020
begin
	# plotting data
	scatter(𝒮.X, 𝒮.Y, legend=false, xlabel="X", ylabel="Y")
	#plotting model
	plot!(1:7, predict(lin_reg))
end

# ╔═╡ 0fb142cd-ac45-42e8-b5d1-24c6142392bb
md"""
What this plot quickly shows is that linear regression may end up predicting values outside the `[0,1]` interval. For an example like this, we will use logistic regression. Interestingly, a generalized linear model (https://en.wikipedia.org/wiki/Generalized_linear_model) unifies concepts like linear regression and logistic regression, and the [GLM.jl](https://github.com/JuliaStats/GLM.jl) package allows you to apply either of these regressions easily by specifying the `distribution family` and the `link` function.

To apply logistic regression via the `GLM.jl` package, you can readily use the `Binomial` family and the `LogitLink` link function.

### Load data

Firstly, let's load some data from the [RDatasets.jl](https://github.com/JuliaStats/RDatasets.jl) package...
"""

# ╔═╡ 3caa5331-2507-4299-b2e4-50e8a5cbbe46
cats = dataset("MASS", "cats")

# ╔═╡ 22c201af-9437-4dde-8653-508ebf514008
md"""
### Data wrangling

We will map the sex of each cat to a binary `0/1` value.
"""

# ╔═╡ 16367b7a-65ea-4fa0-abad-a1cc25c60ee6
# getting the map
sex_map = labelmap(cats.Sex)

# ╔═╡ 1d7ef9db-dd15-4ed4-841a-b4baf02302a1
# encoding Sex column
sex_enconde = labelencode(sex_map, cats.Sex)

# ╔═╡ f42bdc3a-943b-4157-853d-860af52e059c
# plotting data by Sex
scatter(cats.BWt, cats.HWt, color=sex_enconde, legend=false)

# ╔═╡ 310f2282-d6dd-405d-8399-e75fdc75c27c
md"""
Females (blue) seem to be more present in the lower left corner and Males (orange) seem to be present in the top right corner.

We will use the `HWt` variable as the only feature and the `encoded_sex` as the target. Let's convert labels from `1/2` into `0/1`...
"""

# ╔═╡ cc911a53-9139-4c90-b02b-15f30e6772f0
cats_df = DataFrame(X=cats[!,:HWt], Y=sex_enconde .- 1)

# ╔═╡ 79bbe213-c15a-4ecb-aef7-107846234f1d
md"""
### Modeling data

Let's run a logistic regression model on this data.
"""

# ╔═╡ 11b1f5d1-7ae2-43dd-bd35-5b0ba6fde6f4
# logistic model
probit = glm(@formula(Y ~ X), cats_df, Binomial(), LogitLink())

# ╔═╡ e3e9d570-af57-4065-84ec-b61400fa449d
begin
	# plotting ground truth data
	scatter(cats_df.X, cats_df.Y, label="Ground Truth", color=6)

	# plotting predictions
	scatter!(cats_df.X, predict(probit), label="Predictions", color=7,
			 xlabel="Heart Rate", ylabel="Pr(Sex=M)")
end

# ╔═╡ dfc76a4c-4e2b-4588-92e1-2209507fd670
md"""
**Note:** as you can see, contrary to the linear regression case, the predicted values do not go beyond `1`.
"""

# ╔═╡ 31308555-dfc9-4e28-9108-8a83d65b7770
md"""
## Non-linear regression

Finally, sometimes you may have a set of points and the goal is to fit a non-linear function (maybe a quadratic function, a cubic function, an exponential function...). The way we would solve such a problem is by minimizing the least square error between the fitted function and the observations we have. We will use the package [LsqFit.jl](https://github.com/JuliaNLSolvers/LsqFit.jl) for this task.

**Note:** this problem is usually modeled as a numerical optimizaiton problem.
"""

# ╔═╡ 4a7d2c25-04b1-4ccb-85cd-7276902a9679
md"""
### Generate data

We will first set up our data and plot it...
"""

# ╔═╡ 7203ba1b-ce21-4c53-9dce-358df13df21e
begin
	x = 0:0.05:10
	y = 1 * exp.(-x * 2) + 2 * sin.(0.8 * π * x) + 0.15 * randn(length(x));
	scatter(x, y, legend=false)
end

# ╔═╡ 775e6fbe-8c10-4f94-b53e-618e4dd551fc
md"""
### Modeling data

Then, we set up the model with `model(x,p)`. The vector `p` is what to be estimated given a set of values `x`.
"""

# ╔═╡ 5d75ebce-1181-4490-ae10-8bcc08456fdf
begin
	@. model(v, p) = p[1] * exp(-v * p[2]) + p[3] * sin(0.8 * π *v)
	p0 = [0.5, 0.5, 0.5]
	myfit = curve_fit(model, x, y, p0)
end

# ╔═╡ a94163fe-e1ac-483a-bc2d-5d7795cfe188
md"""
⚠ A note about `curve_fit`: this function can take multiple other inputs, for instance the Jacobian of what you are trying to fit. We don't dive into these details here, but be sure to check out the LsqFit.jl package to see what other things you can can pass to create a better fit.

Also note that julia has multiple packages that allow you to create Jacobians so you don't have to write them yourself. Two such packages are [FiniteDifferences.jl](https://github.com/JuliaDiff/FiniteDifferences.jl) or [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl).

↪ Back to our example. Let's define the paramenter variable and then calculate the predictions `ŷ`...
"""

# ╔═╡ 6a965e6d-4724-47c8-a3bb-bfec7d6706c9
begin
	p = myfit.param
	ŷ = p[1] * exp.(-x * p[2]) + p[3] * sin.(0.8 * π * x)
end

# ╔═╡ c7d07a47-3f3f-482e-9c06-78877b91e380
md"""
We are ready now to plot the curve we have generated...
"""

# ╔═╡ 6ca3290d-5980-4749-b309-fd84d73b7bb5
begin
	scatter(x, y, legend=false, xlabel="X", ylabel="Y")
	plot!(x, ŷ)
end

# ╔═╡ Cell order:
# ╟─d16725a0-6cd3-11ec-0998-fd4af26f10dd
# ╟─91471950-703f-4927-8bce-f4478f4a53ba
# ╟─a53d44c5-75c4-4508-aa19-33e838c9e0bb
# ╟─e9458f99-4ec0-4167-b771-3090dfec6320
# ╟─784c5b0e-59dc-4b2d-bf40-112261ac8e85
# ╠═65e22edc-ef5c-412f-82f4-1c84da209d6d
# ╟─db2c3459-2127-4b46-8e73-546850767dc0
# ╠═3bb8c382-f97d-48dc-8cbe-99bf77de0784
# ╟─75059ff4-6b9d-4e01-9aed-33ea8b0e9a71
# ╠═e83331a1-8655-4de0-a36f-b74b5e43674b
# ╟─d81e448b-b4e4-4d6f-8ca7-42c41569977d
# ╠═a047abfb-178c-40ac-8f8c-c3a8a73c1849
# ╠═208af138-5131-4034-885b-c41f6ccf36f5
# ╠═39884fc4-c41c-46e5-a92a-44d766243715
# ╟─13719feb-60ae-433b-bf1e-df6ac045215d
# ╠═10c072d4-b793-458c-b947-d5dc3aae5e25
# ╠═68a0574f-5b17-48ae-b8c4-a03c7e385561
# ╠═ecb1cb29-2b57-4197-9327-e762fc876281
# ╟─056b28f3-8da4-4a94-a33c-72b30c8aa0aa
# ╟─9a5dcad1-9365-423c-b27a-a6c206850848
# ╟─a89ec73a-5e35-49ad-8829-dc2ecc7e8aac
# ╠═aa78c706-8860-46f1-9982-ea10bdc4fccb
# ╟─618ab9ee-ee9c-4dbb-b48f-ebeeb15549c6
# ╠═b530e842-bfa4-4c41-9cc9-4aa629e12f99
# ╠═041e87c7-deee-46c1-85be-ea47df16217c
# ╟─19a44596-2367-4964-b21e-7d84ba147915
# ╠═ab8cf8d2-9ac6-4361-b168-243f781e7dfa
# ╟─215483a1-6059-437d-a4b6-ca3cd897f13b
# ╠═b075d090-91aa-4d1f-bb9f-a93dd65413e7
# ╠═1b8f7b0d-a843-4af7-90b7-4a64d74edab0
# ╟─39456bea-24dc-482b-9d8b-33a9f2ed715a
# ╠═f195db71-f2a4-4f5c-bfa4-2d9ec1f23a32
# ╟─2da58c4e-8db4-4d80-aa29-f3ec10952f73
# ╠═93179d9a-5bf4-488b-b710-28a5f56b2e26
# ╟─560a134f-1175-4d56-83af-1ad6bbe1a99f
# ╠═0c4a00f6-69f2-4fc3-a7f8-ca6e38404ac9
# ╠═a02a9fdf-46bc-4ea1-bb29-b1731a66c057
# ╠═63b596be-9392-4efb-8d31-7f1d6e9c7254
# ╠═feced955-2e69-44ae-bc2a-82d4f2136b38
# ╟─831484f5-6cdc-4641-8adf-117a33c518c2
# ╠═b94019d8-a7d8-4781-819b-28eeac823812
# ╟─4735fc87-68cb-4617-aa7c-6139ba7b7c0a
# ╠═ec7ae9e6-b191-4578-8a47-9d2245280316
# ╟─216dc17e-54c1-4520-890d-8663165b10dd
# ╠═dddfb8fb-7047-4417-8dd2-569a93277f0c
# ╟─5116c117-75e4-4989-b06f-2883bba48330
# ╠═3a86069b-3f24-4eec-a7ae-f77189a4fd77
# ╠═453385dd-8bda-42b8-9b4a-6a4963b21346
# ╠═bd31693d-dee9-4534-b4d4-3904997e8020
# ╟─0fb142cd-ac45-42e8-b5d1-24c6142392bb
# ╠═3caa5331-2507-4299-b2e4-50e8a5cbbe46
# ╟─22c201af-9437-4dde-8653-508ebf514008
# ╠═16367b7a-65ea-4fa0-abad-a1cc25c60ee6
# ╠═1d7ef9db-dd15-4ed4-841a-b4baf02302a1
# ╠═f42bdc3a-943b-4157-853d-860af52e059c
# ╟─310f2282-d6dd-405d-8399-e75fdc75c27c
# ╠═cc911a53-9139-4c90-b02b-15f30e6772f0
# ╟─79bbe213-c15a-4ecb-aef7-107846234f1d
# ╠═11b1f5d1-7ae2-43dd-bd35-5b0ba6fde6f4
# ╠═e3e9d570-af57-4065-84ec-b61400fa449d
# ╟─dfc76a4c-4e2b-4588-92e1-2209507fd670
# ╟─31308555-dfc9-4e28-9108-8a83d65b7770
# ╟─4a7d2c25-04b1-4ccb-85cd-7276902a9679
# ╠═7203ba1b-ce21-4c53-9dce-358df13df21e
# ╟─775e6fbe-8c10-4f94-b53e-618e4dd551fc
# ╠═5d75ebce-1181-4490-ae10-8bcc08456fdf
# ╟─a94163fe-e1ac-483a-bc2d-5d7795cfe188
# ╠═6a965e6d-4724-47c8-a3bb-bfec7d6706c9
# ╟─c7d07a47-3f3f-482e-9c06-78877b91e380
# ╠═6ca3290d-5980-4749-b309-fd84d73b7bb5

max_stars_repo_name: fnaghetini/DataScience

max_issues_repo_name: fnaghetini/DataScience

max_forks_repo_name: fnaghetini/DataScience


text:
module ExpectationPropagationTest

using Test
using ForneyLab
using ForneyLab: SoftFactor, generateId, addNode!, associate!, inferUpdateRule!, outboundType, isApplicable, setTargets!
using ForneyLab: EPProbitIn1PG, SPGaussianMeanVarianceOutNPP, SPClamp, VBGaussianMeanPrecisionOut, SPProbitOutNG, messagePassingSchedule

# Integration helper
mutable struct MockNode <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Int,Interface}

    function MockNode(vars::Vector{Variable}; id=generateId(MockNode))
        n_interfaces = length(vars)
        self = new(id, Array{Interface}(undef, n_interfaces), Dict{Int,Interface}())
        addNode!(currentGraph(), self)

        for idx = 1:n_interfaces
            self.i[idx] = self.interfaces[idx] = associate!(Interface(self), vars[idx])
        end

        return self
    end
end

@expectationPropagationRule(:node_type     => MockNode,
                            :outbound_type => Message{Gaussian},
                            :inbound_types => (Message{PointMass}, Message{Gaussian}),
                            :outbound_id   => 2,
                            :name          => EPMockIn1GP)

@testset "@expectationPropagationRule" begin
    @test EPMockIn1GP <: ExpectationPropagationRule{MockNode}
    @test isApplicable(EPMockIn1GP, [Message{PointMass}, Message{Gaussian}], 2)
    @test !isApplicable(EPMockIn1GP, [Message{PointMass}, Message{Gaussian}, Message{Gaussian}], 2)
end

@testset "inferUpdateRule!" begin
    FactorGraph()
    @RV m ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    nd = MockNode([constant(0.0), m])
    inferred_outbound_types = Dict(nd.i[2].partner => Message{Gaussian}, nd.i[1].partner => Message{PointMass})

    entry = ScheduleEntry(nd.i[2], ExpectationPropagationRule{MockNode})
    inferUpdateRule!(entry, entry.message_update_rule, inferred_outbound_types)

    @test entry.message_update_rule == EPMockIn1GP
end

@testset "messagePassingSchedule" begin
    fg = FactorGraph()
    m = Variable()
    nd_m = GaussianMeanVariance(m, constant(0.0), constant(1.0))
    z = Variable[]
    nd_z = FactorNode[]
    for i = 1:3
        z_i = Variable()
        nd_z_i = Probit(z_i, m)
        placeholder(z_i, :y, index=i)
        push!(z, z_i)
        push!(nd_z, nd_z_i)
    end

    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg)
    pf.ep_sites = Set{Interface}([nd_z_i.i[:in1] for nd_z_i in nd_z])
    setTargets!(pf, pfz, target_variables=Set{Variable}([m]))
    schedule = messagePassingSchedule(pf)

    @test length(schedule) == 15
    @test schedule[5] == ScheduleEntry(nd_z[1].i[:in1], EPProbitIn1PG)
    @test schedule[8] == ScheduleEntry(nd_z[2].i[:in1], EPProbitIn1PG)
    @test schedule[3] == ScheduleEntry(nd_m.i[:out], SPGaussianMeanVarianceOutNPP)
    @test schedule[11] == ScheduleEntry(nd_z[3].i[:in1], EPProbitIn1PG)
end

@testset "messagePassingSchedule" begin
    fg = FactorGraph()
    m = Variable()
    nd_m = GaussianMeanVariance(m, constant(0.0), constant(1.0))
    w = Variable()
    nd_w = Gamma(w, constant(0.01), constant(0.01))
    y = Variable()
    nd_y = GaussianMeanPrecision(y, m, w)
    z = Variable()
    nd_z = Probit(z, y)
    placeholder(z, :z)

    pfz = PosteriorFactorization()
    q_y_z = PosteriorFactor([y, z])
    q_y_z.ep_sites = Set{Interface}([nd_z.i[:in1]])
    setTargets!(q_y_z, pfz, external_targets=true)
    schedule = messagePassingSchedule(q_y_z)

    @test length(schedule) == 3
    @test ScheduleEntry(nd_y.i[:out], VBGaussianMeanPrecisionOut) in schedule
    @test ScheduleEntry(nd_z.i[:out].partner, SPClamp{Univariate}) in schedule
    @test ScheduleEntry(nd_z.i[:in1], EPProbitIn1PG) in schedule
end

@testset "messagePassingAlgorithm" begin
    fg = FactorGraph()
    m = Variable()
    nd_m = GaussianMeanVariance(m, constant(0.0), constant(1.0))
    z = Variable[]
    nd_z = FactorNode[]
    for i = 1:3
        z_i = Variable()
        nd_z_i = Probit(z_i, m)
        placeholder(z_i, :y, index=i)
        push!(z, z_i)
        push!(nd_z, nd_z_i)
    end

    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(m, ep_sites=[(:probit_*i, :in1) for i=1:3])

    @test isa(algo, InferenceAlgorithm)
end

@testset "messagePassingAlgorithm" begin
    fg = FactorGraph()
    m = Variable()
    nd_m = GaussianMeanVariance(m, constant(0.0), constant(1.0))
    w = Variable()
    nd_w = Gamma(w, constant(0.01), constant(0.01))
    y = Variable()
    nd_y = GaussianMeanPrecision(y, m, w)
    z = Variable()
    nd_z = Probit(z, y)
    placeholder(z, :z)

    pfz = PosteriorFactorization([y; z], m, w)
    algo = messagePassingAlgorithm(m, ep_sites=[(:probit_1, :in1)])

    @test isa(algo, InferenceAlgorithm)
end

end # module


text:
# This script implements a partial differential equation describing changes in
# concentration with respect to longitudinal distance along a stream channel
# due to advection. I try using both the symbolic PDE method (which works)
# and the finite-difference operators method (which doesn't work yet)

##### ----- Define Parameters, Initial Conditions, Boundary Conditions ----- #####
v = 1.0 # Velocity
p = [v] # parameter vector
tspan = (0.0, 10.0) # simulation timespan
xₙ = 10 # physical length of simulation
Δx = 0.01 # grid spacing along x dimension (longitudinal)
N = convert(Int64, xₙ/Δx - 1) # number of interior points
x = range(0, step = Δx, length = (N+2)) # x-points of solution, including boundaries
h(p,t) = zeros(N) # history function describing solute concentrations before the start of the simulation
C₀ = zeros(N) # initial condition
for i in 1:N
  u₀[i] = 1.0
end

##### ----- Automated Finite-Difference for Symbolic PDES Method ----- #####
using ModelingToolkit
@parameters t x v
@variables(u(..))
@derivatives Dt'~t
@derivatives Dx'~x

#1D PDE and boundary conditions
eq = Dt(u(t,x)) ~ -v*Dx(u(t,x))
bcs = [u(0,x) ~ 0,
        u(t,0) ~ 1.0,
        Dx(u(t,10.0)) ~ 0.0]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,convert(Float64, xₙ)),
           x ∈ IntervalDomain(tspan[1],tspan[2])]

# PDE system
pdesys = PDESystem(eq,bcs,domains,[t,x],[u])

# Method of lines discretization
discretization = MOLFiniteDifference(Δx,3) #Third order finite difference approximation with nodes spaced by 0.01

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

# Solve ODE problem
sol = solve(prob,Tsit5(),saveat=1.0)

# Plot results and compare with exact solution
t = sol.t
x = prob.space[2]

plt = plot()
for i in 1:length(t)
    plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]),label="t=$(t[i])")
end
display(plt)

##### ----- Finite Difference Operator Method ----- #####
ΔUpwind = DiffEqOperators.UpwindDifference(1,4,Δx,xₙ) # first derivative operator, fourth-order approximation
Dirichletbc = DiffEqOperators.DirichletBC(1.0,0.0) # boundary condition operator

@advect(ΔUpwind,Dirichletbc)
prob = ODEProblem(advect!, u₀, tspan, p)
alg = Rosenbrock23(autodiff=false)
maxdt = Δx/v
sol = solve(prob, alg, saveat = 1.0, dtmax = 0.8*maxdt)
t = sol.t
u = sol.u

plt = plot()
for i in 1:length(t)
    plot!(x, bc*u[[i]][1] ,label="t=$(t[i])")
end
display(plt)

max_stars_repo_name: elizabethjmohr/hyporheicBGC

max_issues_repo_name: elizabethjmohr/hyporheicBGC

max_forks_repo_name: elizabethjmohr/hyporheicBGC


text:
#BB side
function startbb()
    @async begin
        server = listen(2001)
        while true
            sock = accept(server)
            println("accepted")
            @async while isopen(sock)
                l = deserialize(sock);
                println(typeof(l))
                #print("Read: ");
                println(l);
                serialize(sock,l)
            end
        end
    end
end

#Computer side
clientside = connect(ip"192.168.7.2", 2001)
function getminor(t)
    println(Int64(t))
    t2 = UInt64(floor(UInt64,t/1e9)*1e9)
    tsmall = Int32(t-t2)
end

function runcomp(clientside)
    @async while isopen(clientside)
        l = deserialize(clientside)::Int32;
        #print("Read:");
        t = time_ns()
        #println(t)
        t2 = getminor(t)
        println((t2-l)/1e6)
    end
end
serialize(clientside, getminor(time_ns()))
end # module


text:
"""
`DatePicker` provides UI elements for picking a date, multiple dates, a range of dates,
  or multiple dates ranges. The `DatePicker` module automatically exports `Dates`.

## Example
```julia

```
"""
module DatePickers

using Genie, Stipple, StippleUI, StippleUI.API
using Stipple, Stipple.Reexport

@reexport using Dates

import Genie.Renderer.Html: HTMLString, normal_element, template

export datepicker, DateRange, DatePicker

function __init__()
  Genie.Renderer.Html.register_normal_element("q__date", context = Genie.Renderer.Html)
end

"""
    DateRange

Represents a date interval, between `start` and `stop`, with a 1 day step.
"""
Base.@kwdef mutable struct DateRange
  start::Date = today()
  stop::Date = today()
end

DateRange(dr::StepRange{Date,Day}) = DateRange(dr.start, dr.stop)

"""
    datepicker()

Renders a date picker (calendar) input element.
If the `fieldname` references a `Vector{Date}`, the `multiple` keyword parameter must be passed as `true`.
If the `fieldname` references a `DateRange`, the `range` keyword parameter must be passed as `true`.
If the `fieldname` references a `Vector{DateRange}`, both the `multiple` and the `range` keyword parameters must be passed as `true`.
"""
function datepicker(
                    fieldname::Union{Symbol,Nothing} = nothing,
                    args...;
                    wrap::Function = StippleUI.DEFAULT_WRAPPER,
                    mask::String = "YYYY-MM-DD",
                    content::Union{String,Vector,Function} = "",
                    kwargs...)

  wrap() do
    Genie.Renderer.Html.q__date([isa(content, Function) ? content() : join(content)],
            args...;
            attributes(
                  [ :fieldname => fieldname,
                    :mask => mask,
                    kwargs...
                  ], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
  end
end

datepicker( content::Union{Vector,Function},
            fieldname::Union{Symbol,Nothing} = nothing,
            args...;
            wrap::Function = StippleUI.DEFAULT_WRAPPER,
            mask::String = "YYYY-MM-DD",
            kwargs...) = datepicker(fieldname, args...; wrap = wrap, mask = mask, content = content, kwargs...)

mutable struct DatePicker
  fieldname
  args
  wrap
  mask
  content
  kwargs

  DatePicker( fieldname::Union{Symbol,Nothing} = nothing,
              args...;
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              mask::String = "YYYY-MM-DD",
              content::Union{String,Vector,Function} = "",
              kwargs...) = new(fieldname, args, wrap, mask, kwargs)
end

Base.string(dp::DatePicker) = datepicker(dp.fieldname, dp.args...; dp.wrap, dp.mask, dp.content, dp.kwargs...)

# internals

function Base.convert(::Type{Date}, d::String) :: Date
  Date(d)
end

function Stipple.render(dr::DateRange, _::Union{Symbol,Nothing} = nothing)
  Dict(:from => dr.start, :to => dr.stop)
end

function Stipple.render(vdr::Vector{DateRange}, _::Union{Symbol,Nothing} = nothing)
  [ Dict(:from => dr.start, :to => dr.stop) for dr in vdr ]
end

function Base.parse(::Type{DateRange}, d::Dict{String,Any})
  DateRange(d["from"], d["to"])
end

Base.convert(::Type{DateRange}, d::Dict{String,Any}) = parse(DateRange, d)

end