text:
using Lasso, FactCheck
DATADIR = joinpath(dirname(@__FILE__), "data")

function diffmatslow(order, n)
    D = spdiagm((fill(-1., n-1), fill(1., n-1)), (0, 1))
    for i = 1:order
        D = spdiagm((fill(-1., n-i-1), fill(1., n-i-1)), (0, 1))*D
    end
    D
end

# Test that DifferenceMatrix behaves as expected
facts("DifferenceMatrix") do
    for order in (1, 2, 3)
        context("order = $(order)") do
            D1 = diffmatslow(order, 100)
            D2 = Lasso.TrendFiltering.DifferenceMatrix{Float64}(order, 100)
            x = randn(100)
            @fact D2'D2 => D1'D1
            @fact D2*x => roughly(D1*x)
            @fact D2'*x[1:size(D1, 1)] => roughly(D1'*x[1:size(D1, 1)])
        end
    end
end

# Test against results from glmgen
lakehuron = readcsv(joinpath(DATADIR, "LakeHuron.csv"); header=true)[1][:, 3]
facts("TrendFiltering") do
    for order in (1, 2, 3)
        context("order = $(order)") do
            for lambda in (1., 10., 100.)
                context("lambda = $(lambda)") do
                    @fact coef(fit(TrendFilter, lakehuron, order, lambda; tol=1e-9)) =>
                          roughly(vec(readcsv(joinpath(DATADIR, "LakeHuron_order_$(order)_lambda_$(lambda).csv"))))
                end
            end
        end
    end
end


text:
using Base.Meta
using WebAssembly, WebAssembly.Instructions
using WebAssembly: WType, Func, Module, FuncType, Func, Table, Mem, Global, Elem, Data, Import, Export

# This struct stores state as a wasm module is built up.
struct ModuleState
  imports::Dict{Symbol, Import}
  exports::Dict{Symbol, Export}
  funcs::Dict{Symbol, Func}
  data::Dict{Symbol, Data}
end
ModuleState() = ModuleState(Dict{Symbol, Import}(), Dict{Symbol, Export}(), Dict{Symbol, Func}(), Dict{Symbol, Data}())

walk(x, inner, outer) = outer(x)

function walk(x::Expr, inner, outer)
  y = Expr(x.head, map(inner, x.args)...)
  y.typ = x.typ
  return outer(y)
end

prepostwalk(f, g, x) = walk(f(x), x -> prepostwalk(f, g, x), g)
prewalk(f, x)  = prepostwalk(f, identity, x)
postwalk(f, x) = prepostwalk(identity, f, x)

deref(x) = x
deref(x::GlobalRef) = getfield(x.mod, x.name)

exprtype(code::CodeInfo, x) = typeof(deref(x))
exprtype(code::CodeInfo, x::Expr) = x.typ
exprtype(code::CodeInfo, x::QuoteNode) = typeof(x.value)
exprtype(code::CodeInfo, x::SSAValue) = code.ssavaluetypes[x.id+1]
exprtype(code::CodeInfo, x::SlotNumber) = code.slottypes[x.id]

# We don't want SSAValues on a stack machine.
# If it's only used once, just inline it.
# If more than once, we have to turn it into a local.

function ssacounts(code)
  counts = Dict{SSAValue,Any}()
  inc(x) = counts[x] = get(counts, x, 0) + 1
  for c in code
    isexpr(c, :(=)) && (c = c.args[2])
    prewalk(x -> (x isa SSAValue && inc(x); x), c)
  end
  return counts
end

apply_rhs(f, x) =
    isexpr(x, :(=)) ?
      Expr(:(=), x.args[1], f(x.args[2])) :
      f(x)

function inlinessa(code)
  counts = ssacounts(code)
  values = Dict{SSAValue,Any}()
  code′ = []
  for c in code
    c = apply_rhs(c) do x
      prewalk(x -> get(values, x, x), x)
    end
    isexpr(c, :(=)) && c.args[1] isa SSAValue && counts[c.args[1]] == 1 ?
      (values[c.args[1]] = c.args[2]) :
      push!(code′, c)
  end
  return code′
end

function ssalocals(cinfo, ex)
  n = length(cinfo.slotnames)
  ls = Dict{SSAValue,SlotNumber}()
  ex = prewalk(ex) do x
    x isa SSAValue || return x
    haskey(ls, x) && return ls[x]
    push!(cinfo.slottypes, exprtype(cinfo, x))
    ls[x] = SlotNumber(n += 1)
  end
  return ex
end

function rmssa(c)
  is = copy(c.code)
  while is[1] isa NewvarNode shift!(is) end
  ex = Expr(:block, inlinessa(is)...)
  ssalocals(c, ex)
end

# Julia -> WASM expression

wasmfuncs = Dict()

int_unary_ops = [
  (:cttz_int, :ctz)
]

int_binary_ops = [
  (:add_int, :add),
  (:sub_int, :sub),
  (:mul_int, :mul),
  (:and_int, :and),
  (:or_int,  :or),
  (:xor_int, :xor),
  (:shl_int, :shl),
  (:lshr_int, :shr_u),
  (:ashr_int, :shr_s),
  (:(===),   :eq),
  (:slt_int, :lt_s),
  (:sle_int, :le_s),
  (:ult_int, :lt_u),
  (:ule_int, :le_u),
]

float_unary_ops = [
  (:abs_float, :abs),
  (:neg_float, :neg),
  (:rint_llvm, :nearest),
]

float_binary_ops = [
  (:add_float, :add),
  (:sub_float, :sub),
  (:mul_float, :mul),
  (:div_float, :div),
  (:eq_float,  :eq),
  (:ne_float,  :ne),
  (:lt_float,  :lt),
  (:le_float,  :le),
]

for (j, w) in vcat(int_unary_ops, float_unary_ops)
  wasmfuncs[GlobalRef(Base, j)] = function (T)
    Op(WType(T), w)
  end
end

for (j, w) in vcat(int_binary_ops, float_binary_ops)
  wasmfuncs[GlobalRef(Base, j)] = function (A,B)
    Op(WType(A), w)
  end
end

wasmfuncs[GlobalRef(Base, :not_int)] = function (A)
  Op(WType(A), :eqz)
end

wasmfunc(f, xs...) = wasmfuncs[f](xs...)

wasmcalls = Dict()

# Intrinsics without wasm equivalents

wasmcalls[GlobalRef(Base, :muladd_float)] = function (i, x, y, z)
  xy = Expr(:call, GlobalRef(Base, :mul_float), x, y)
  xy.typ = exprtype(i, x)
  Expr(:call, nop,
    Expr(:call, GlobalRef(Base, :add_float), xy, z))
end

wasmcalls[GlobalRef(Base, :neg_int)] = function (i, x)
  Expr(:call, nop,
    Expr(:call, GlobalRef(Base, :sub_int), exprtype(i, x)(0), x))
end

wasmcalls[GlobalRef(Base, :flipsign_int)] = function (i, x, y)
  T = exprtype(i, x)
  cond = Expr(:call, GlobalRef(Base, :slt_int), x, T(0))
  sign = Expr(:call, GlobalRef(Base, :select_value), cond, T(-1), T(1))
  sign.typ = T
  Expr(:call, nop,
    Expr(:call, GlobalRef(Base, :mul_int), sign, x))
end

# TODO: should probably not do this
wasmcalls[GlobalRef(Base, :check_top_bit)] = function (i, x)
  nop
end

# More complex intrinsics

wasmcalls[GlobalRef(Base, :select_value)] = function (i, c, a, b)
  return Expr(:call, Select(), a, b, c)
end

wasmcalls[GlobalRef(Base, :bitcast)] = function (i, T, x)
  T isa GlobalRef && (T = getfield(T.mod, T.name))
  X = exprtype(i, x)
  WType(T) == WType(X) && return Expr(:call, nop, x)
  @assert sizeof(T) == sizeof(X)
  Expr(:call, Convert(WType(T), WType(X), :reinterpret), x)
end

wasmcalls[GlobalRef(Base, :fptosi)] = function (i, T, x)
  T isa GlobalRef && (T = getfield(T.mod, T.name))
  X = exprtype(i, x)
  Expr(:call, Convert(WType(T), WType(X), :trunc_s), x)
end

wasmcalls[GlobalRef(Base, :sitofp)] = function (i, T, x)
  T isa GlobalRef && (T = getfield(T.mod, T.name))
  X = exprtype(i, x)
  Expr(:call, Convert(WType(T), WType(X), :convert_s), x)
end

wasmcall(i, f, xs...) =
  haskey(wasmcalls, f) ? wasmcalls[f](i, xs...) :
  Expr(:call, wasmfunc(f, exprtype.(i, xs)...), xs...)

isprimitive(x) = false
isprimitive(x::GlobalRef) =
  deref(x) isa Core.IntrinsicFunction ||
  deref(x) isa Core.Builtin

function lowercalls(m::ModuleState, c::CodeInfo, code)
  prewalk(code) do x
    if isexpr(x, :call) && deref(x.args[1]) == Base.throw
      unreachable
    elseif (isexpr(x, :call) && isprimitive(x.args[1]))
      wasmcall(c, x.args...)
    elseif isexpr(x, :invoke)
      lower_invoke(m, x.args)
    elseif isexpr(x, :foreigncall)
      lower_ccall(m, x.args)
    elseif isexpr(x, :(=)) && x.args[1] isa SlotNumber
      Expr(:call, SetLocal(false, x.args[1].id-2), x.args[2])
    elseif x isa SlotNumber
      Local(x.id-2)
    elseif isexpr(x, :gotoifnot)
      Expr(:call, Goto(true, x.args[2]),
           Expr(:call, GlobalRef(Base, :not_int), x.args[1]))
    elseif x isa GotoNode
      Goto(false, x.label)
    elseif x isa LabelNode
      Label(x.label)
    elseif isexpr(x, :return)
      Expr(:call, Return(), x.args[1])
    else
      x
    end
  end
end

function lower_invoke(m::ModuleState, args)
  # This lowers the function invoked. `args` is the Any[] from the :invoke Expr.
  # If the function has not been compiled, compile it.
  # Generate the WASM call.
  tt = argtypes(args[1])
  name = createfunname(args[1], tt)
  if !haskey(m.funcs, name)
    mi = args[1]
    ci = Base.uncompressed_ast(mi.def, mi.inferred)
    R = mi.rettype
    m.funcs[name] = code_wasm(m::ModuleState, name, tt, ci, R)
  end
  return Expr(:call, Call(name), args[3:end]...)
end

function lower_ccall(m::ModuleState, args)
  (fnname, env) = args[1]
  name = Symbol(env, :_, fnname)
  m.imports[name] = Import(env, fnname, :func, map(WType, args[3]), WType(args[2]))
  return Expr(:call, Call(name), args[4:2:end]...)
end

argtypes(x::Core.MethodInstance) = Tuple{x.specTypes.parameters[2:end]...}
argtypes(x::Method) = Tuple{x.sig.parameters[2:end]...}

createfunname(fun::Symbol, argtypes) = Symbol(fun, "_", join(collect(argtypes.parameters), "_"))
createfunname(fun::Function, argtypes) = createfunname(typeof(fun), argtypes)
createfunname(funtyp::DataType, argtypes) = Symbol(funtyp, "_", join(collect(argtypes.parameters), "_"))
createfunname(mi::Core.MethodInstance, argtypes) = createfunname(mi.def.sig.parameters[1], argtypes)

basename(f::Function) = Base.function_name(f)
basename(f::Core.IntrinsicFunction) = Symbol(unsafe_string(ccall(:jl_intrinsic_name, Cstring, (Core.IntrinsicFunction,), f)))
basename(x::GlobalRef) = x.name
basename(m::Core.MethodInstance) = basename(m.def)
basename(m::Method) = m.name == :Type ? m.sig.parameters[1].parameters[1].name.name : m.name

iscontrol(ex) = isexpr(ex, :while) || isexpr(ex, :if)

lower(m::ModuleState, c::CodeInfo) = lowercalls(m, c, rmssa(c))

# Convert to WASM instructions

towasm_(m::ModuleState, xs, is = Instruction[]) = (foreach(x -> towasm(m, x, is), xs); is)

function towasm(m::ModuleState, x, is = Instruction[])
  if x isa Instruction
    push!(is, x)
  elseif isexpr(x, :block)
    push!(is, Block(towasm_(m, x.args)))
  elseif isexpr(x, :call) && x.args[1] isa Instruction
    foreach(x -> towasm(m, x, is), x.args[2:end])
    push!(is, x.args[1])
  elseif isexpr(x, :if)
    towasm(m, x.args[1], is)
    push!(is, If(towasm_(m, x.args[2].args), towasm_(m, x.args[3].args)))
  elseif isexpr(x, :while)
    push!(is, Block([Loop(towasm_(m, x.args[2].args))]))
  elseif deref(x) isa Number
    push!(is, Const(deref(x)))
  elseif x isa LineNumberNode || isexpr(x, :inbounds) || isexpr(x, :meta)
  else
    error("Can't convert to wasm: $x")
  end
  return is
end

funname(f::Function) = Base.function_name(f)
funname(s::Symbol) = s

function code_wasm(m::ModuleState, ex, A)
  cinfo, R = code_typed(ex, A)[1]
  code_wasm(m, createfunname(ex, A), A, cinfo, R)
end

function code_wasm(m::ModuleState, name::Symbol, A, cinfo::CodeInfo, R)
  body = towasm_(m, lower(m, cinfo).args) |> Block |> WebAssembly.restructure |> WebAssembly.optimise
  Func(name,
       [WType(T) for T in A.parameters],
       [WType(R)],
       [WType(P) for P in cinfo.slottypes[length(A.parameters)+2:end]],
       body)
end

macro code_wasm(ex)
  isexpr(ex, :call) || error("@code_wasm f(xs...)")
  :(code_wasm(ModuleState(), $(esc(ex.args[1])), Base.typesof($(esc.(ex.args[2:end])...))))
end

"""
    wasm_module(funpairlist)

Return a compiled WebAssembly ModuleState that includes every function defined by `funpairlist`.

`funpairlist` is a vector of pairs. Each pair includes the function to include and
a tuple type of the arguments to that function. For example, here is the invocation to
return a wasm ModuleState with the `mathfun` and `anotherfun` included.

    m = wasm_module([mathfun => Tuple{Float64},
                     anotherfun => Tuple{Int, Float64}])
"""
function wasm_module(funpairlist)
  m = ModuleState()
  for (fun, tt) in funpairlist
    internalname = createfunname(fun, tt)
    exportedname = funname(fun)
    m.funcs[internalname] = Func(:name, [], [], [], Block([]))
    m.funcs[internalname] = code_wasm(m, fun, tt)
    m.exports[exportedname] = Export(exportedname, internalname, :func)
  end
  if length(m.data) > 0
    m.exports[:memory] = Export(:memory, :memory, :memory)
  end
  return Module(FuncType[], collect(values(m.funcs)), Table[], [Mem(:m, 1, nothing)], Global[], Elem[],
                collect(values(m.data)), Ref(0), collect(values(m.imports)), collect(values(m.exports)))
end

ext: jl


text:
# Function for FAVAR by Bernanke et al. (2015)
import Pkg;
using LinearAlgebra, Distributions, Statistics, MultivariateStats;
using ProgressMeter;
using Gadfly, Colors;

function func_std(mat)
    # Function to standardize data
    T, k = size(mat);
    μ = mean(mat; dims = 1);
    σ = std(mat; dims = 1);
    return (mat .- μ) ./ σ;
end

function func_favar(yₘ, yₚ, p, ζ)
    # yₘ: data - main block
    # yₚ: data - factor block - to be
    # p: lag term for VAR
    # ζ: Number of factors to be extracted

    T, kₘ = size(yₘ);

    # For principal component analysis (PCA), always use standardized data.
    res_pca = fit(PCA, func_std(yₚ); maxoutdim = ζ);
    𝔽ₒ = projection(res_pca);

    Q = yₚ';
    Z = [yₘ 𝔽ₒ]';
    ℾₒ = (Q * Z')/(Z * Z');
    ℾₒ𝑦 = ℾₒ[:, (1:kₘ)];
    x̃ = yₚ - (ℾₒ𝑦 * yₘ')';

    ssrₒ = sum((Q - ℾₒ * Z).^2, dims = 2)

    # Start the iteraction process
    𝔽₁ = 𝔽ₒ;
    ℾ₁ = ℾₒ; ℾ₁𝑦 = ℾₒ𝑦;
    ssr₁ = ssrₒ

    n_ϝ = 0; # Count iterations;
    while sum(abs(ssr₁ .- ssrₒ) > 10^(-6))
        # If no convergence, repeat the process
        # Update the criteria
        𝔽ₒ = 𝔽₁;
        ssrₒ = ssr₁;

        # Redraw component
        res_pca = fit(PCA, func_std(x̃); maxoutdim = ζ);
        𝔽₁ = projection(res_pca);

        Q = yₚ';
        Z = [yₘ 𝔽₁]';
        ℾ₁ = (Q * Z')/(Z * Z');
        ℾ₁𝑦 = ℾ₁[:, (1:kₘ)];
        x̃ = yₚ - (ℾ₁𝑦 * yₘ')';

        ssr₁ = sum((Q - ℾ₁ * Z).^2, dims = 2);

        n_ϝ = n_ϝ + 1;
    end

    # conditional check for factor loading (optional)
    #=
    #Check on elements of factor loading that belongs to first series in the
    # information set matrix - yₚ
    # If an element is negative, multiply whole factor loading with -1 and
    # then respective factor

    for iₚ = (kₘ + 1):size(Gamma_final, 2)
        if ℾ₁[1, iₚ] < 0
            ℾ₁[:, iₚ] = -1 .* ℾ₁[:, iₚ];
            𝔽₁[:, (iₚ - kₘ)] = -1 .* 𝔽₁[:, (iₚ - kₘ)];
        end
    end
    =#

    Q = yₚ';
    Z = [yₘ 𝔽₁]';
    ℾ₁ = (Q * Z')/(Z * Z');
    ℾ₁𝑦 = ℾ₁[:, (1:kₘ)];
    x̃ = yₚ - (ℾ₁𝑦 * yₘ')';

    # compute VAR(p)
    𝔇 = [yₘ 𝔽₁];

    𝚩, 𝞄, 𝝨 = func_VAR(𝔇, p);

    return 𝚩, 𝞄, 𝝨, 𝔽₁, ℾ₁;
end


text:
using ExactWrightFisher

function prepare_WF_dat_1D_2D(data::Dict{Float64,Array{Int64,2}})
    times = data |> keys |> collect |> sort
    return Dict(zip(times, [vec(data[t]) for t in times])), data
end
function prepare_WF_dat_1D_2D(data::Dict{Float64,Array{Int64,1}})
    times = data |> keys |> collect |> sort
    return data, zip(times, [collect(data[t]') for t in times]) |> Dict
end

function precompute_next_terms_ar!(last_sm_max::Integer, new_sm_max::Integer, log_ν_ar::Array{Float64,2}, log_Cmmi_ar::Array{Float64,2}, log_binomial_coeff_ar_offset::Array{Float64,2}, sα, Δt)
    if last_sm_max == 0
        log_binomial_coeff_ar_offset[1,1] = 0
    end
    for sm in (last_sm_max+1):new_sm_max
        for si in 1:sm
            log_ν_ar[sm, si] = logfirst_term_pmmi_no_alloc(si, sm, sα)
            log_Cmmi_ar[sm, si] = logCmmi_arb(sm, si, Δt, sα)
            log_binomial_coeff_ar_offset[sm+1,si+1] = log_binomial_safe_but_slow(sm, si)
        end
        for si in 0:sm
            log_binomial_coeff_ar_offset[sm+1,si+1] = log_binomial_safe_but_slow(sm, si)
        end
    end
end

function Λ_from_Λ_max(Λ_max) where U <: Integer
    # return Base.Iterators.product((0:Λi_max for Λi_max in Λ_max)...)
    return Base.Iterators.product(map(N -> 0:N, Λ_max)...)
end

indices_of_tree_below(m::Union{AbstractArray{U, 1}, Tuple}) where U <: Integer =  Λ_from_Λ_max(m)

function logfirst_term_pmmi_no_alloc(si::Int64, sm::Int64, sα::Number)
    ##### Equivalent to
    # logλm.((sm-si+1):sm, sα) |> sum
    # but the former version does one allocation
    res = 0
    for s in (sm-si+1):sm
        res += logλm.(s, sα)
    end
    return res
end

function logCmmi_arb(sm::Int64, si::Int64, t::Number, sα::Number)::Float64
    # tmp = [-λm(sm-k, sα)*t - log_denominator_Cmmi_nosign(si, k, sm, sα) for k in 0:si]
    # max_tmp = maximum(tmp)
    return ExactWrightFisher.signed_logsumexp_arb((-λm(sm-k, sα)*t - log_denominator_Cmmi_nosign(si, k, sm, sα) for k in 0:si), (sign_denominator_Cmmi(k) for k in 0:si))[2]
    # return ExactWrightFisher.signed_logsumexp_arb((-λm(sm-k, sα)*t - log_denominator_Cmmi_nosign(si, k, sm, sα) for k in 0:si), take(cycle([1,-1]), si+1))[2]
    # return ExactWrightFisher.signed_logsumexp_arb(-λm.(sm .- (0:si), sα) .* t - log_denominator_Cmmi_nosign.(si, 0:si, sm, sα) , sign_denominator_Cmmi.(0:si) )[2]
end

function logλm(sm::Int64, sα::Number)
    return log(sm) + log(sm + sα - 1) - log(2)
end

function λm(sm::Int64, sα::Number)
    return sm * (sm + sα - 1)/2
end

function log_denominator_Cmmi_nosign(si::Int64, k::Int64, sm::Int64, sα::Number)
    if k==0
        return SpecialFunctions.logfactorial(si) - si*log(2) + log_descending_fact_no0(2*sm + sα - 2, si)
    elseif k==si
        return SpecialFunctions.logfactorial(si) - si*log(2) + log_descending_fact_no0(2*sm + sα - si-1, si)
    else
        return -1.0 .* si * log(2) + SpecialFunctions.logfactorial(k) + SpecialFunctions.logfactorial(si-k) + log_descending_fact_no0(2*sm + sα - 2*k - 2, si-k) + log_descending_fact_no0(2*sm + sα - k - 1, k)
    end
end

function sign_denominator_Cmmi(k::Int64)
    if iseven(k)
        return 1
    else
        return -1
    end
end

function logpmmi_raw_precomputed(i, m, sm::Integer, si::Integer, t::Number, log_ν_dict::Dict{Tuple{Int64, Int64}, Float64}, log_Cmmi_dict::Dict{Tuple{Int64, Int64}, Float64}, log_binomial_coeff_dict::Dict{Tuple{Int64, Int64}, Float64}) where T <: Integer
    return log_ν_dict[(sm, si)] + log_Cmmi_dict[(sm, si)]  + loghypergeom_pdf_using_precomputed(i, m, si, sm, log_binomial_coeff_dict)
end

function logpmmi_raw_precomputed(i, m, sm::Integer, si::Integer, t::Number, log_ν_ar::Array{Float64,2}, log_Cmmi_ar::Array{Float64,2}, log_binomial_coeff_ar_offset::Array{Float64,2})
    #Would return an error when called on sm = 0, but this should never occur.
    return log_ν_ar[sm, si] + log_Cmmi_ar[sm, si]  + loghypergeom_pdf_using_precomputed(i, m, si, sm, log_binomial_coeff_ar_offset)
end

function logpmmi_precomputed(i, m, sm::Integer, si::Integer, t::Number, sα::Number, precomputed_log_ν, precomputed_log_Cmmi, precomputed_log_binomial_coeff) where T <: Integer
    if si==0
        return -λm(sm, sα)*t
    else
        return logpmmi_raw_precomputed(i, m, sm, si, t, precomputed_log_ν, precomputed_log_Cmmi, precomputed_log_binomial_coeff)
    end
end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


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

using Markdown
using InteractiveUtils

# ╔═╡ 30fe24e6-4584-11eb-1bbe-69bdf11577da
using DataFrames

# ╔═╡ 78006de0-4584-11eb-3990-d344428fb2ee
using CSV

# ╔═╡ 60da9828-4586-11eb-149b-93667e86c30f
md"""
# Data Analysis of Internet News
- [x] Import Data
- [x] Display column names
- [x] Plot the distrbution of sources, authors
- [ ] Make a word cloud.
"""

# ╔═╡ 59890220-4587-11eb-000d-17acdb41c67d
md"Importing the data file, \"articles_data.csv\""

# ╔═╡ 1d24d176-4585-11eb-33fe-a9d9abfc94ba
data = DataFrame(CSV.File("articles_data.csv"))

# ╔═╡ ae1a254c-4589-11eb-3361-8fffaa22f423
md"summary of data"

# ╔═╡ 9b855258-4589-11eb-1000-351ffd4e2a25
md"Removing the rows with missing Data"

# ╔═╡ 0c5a65c0-4589-11eb-1ed4-ebb3eb1a1874
DataFrames.dropmissing!(data)

# ╔═╡ d7941d3a-4587-11eb-0d16-6979760f882d
md"Importing Statistics and Ploting modules"

# ╔═╡ a72321de-458a-11eb-1c45-0ddd309b920b
names(data)

# ╔═╡ b6bf5f22-458a-11eb-0afd-a35685a58d7e
source_name_counts = Dict(title => count(x -> x == title, data.source_name) for title ∈ unique(data.source_name))

# ╔═╡ 9a93c46e-458f-11eb-3f1a-5bd1e2bae0d4
total_number_of_sources = source_name_counts |> values |> sum

# ╔═╡ 3e95df84-4590-11eb-3f87-87b652fc252d
(x -> ("$(x[1]) = $(round(x[2]/total_number_of_sources*100; digits = 2))%")).(Tuple(source_name_counts))

# ╔═╡ cfb5335e-459d-11eb-22a2-b3ead0a592bc
md"The titles' word distrobution."

# ╔═╡ Cell order:
# ╟─60da9828-4586-11eb-149b-93667e86c30f
# ╠═30fe24e6-4584-11eb-1bbe-69bdf11577da
# ╠═78006de0-4584-11eb-3990-d344428fb2ee
# ╟─59890220-4587-11eb-000d-17acdb41c67d
# ╟─1d24d176-4585-11eb-33fe-a9d9abfc94ba
# ╠═ae1a254c-4589-11eb-3361-8fffaa22f423
# ╟─9b855258-4589-11eb-1000-351ffd4e2a25
# ╠═0c5a65c0-4589-11eb-1ed4-ebb3eb1a1874
# ╟─d7941d3a-4587-11eb-0d16-6979760f882d
# ╠═a72321de-458a-11eb-1c45-0ddd309b920b
# ╟─b6bf5f22-458a-11eb-0afd-a35685a58d7e
# ╟─9a93c46e-458f-11eb-3f1a-5bd1e2bae0d4
# ╟─3e95df84-4590-11eb-3f87-87b652fc252d
# ╟─cfb5335e-459d-11eb-22a2-b3ead0a592bc

max_stars_repo_name: kawazugit/internet-news-data

max_issues_repo_name: kawazugit/internet-news-data

max_forks_repo_name: kawazugit/internet-news-data


text:
# This file was generated, do not modify it. # hide
rms(pred_dcrm, y[test])


text:
maxdims["dmrg_1d_qns"] = 200:200:1_000
descriptions["dmrg_1d_qns"] = "DMRG, 1D S=1 Heisenberg model\nN = 100, 5 sweeps\nConserve total spin projection symmetry"

function runbenchmark(::Val{:dmrg_1d_qns};
                      maxdim::Int, nsweeps::Int = 5,
                      outputlevel::Int = 0,
                      conserve_qns::Bool = true,
                      N::Int = 100, splitblocks::Bool = false)
  if !conserve_qns && splitblocks
    println("In benchmark dmrg_1d_qns, conserve_qns is $conserve_qns, cannot use with splitblocks $splitblocks.")
    return nothing
  end
  sites = siteinds("S=1", N; conserve_qns = conserve_qns)
  ampo = AutoMPO()
  for j in 1:N-1
    ampo .+= 0.5, "S+", j, "S-", j + 1
    ampo .+= 0.5, "S-", j, "S+", j + 1
    ampo .+=      "Sz", j, "Sz", j + 1
  end
  H = MPO(ampo,sites)
  if splitblocks
    H = ITensors.splitblocks(linkinds, H)
  end
  psi0 = productMPS(sites, n -> isodd(n) ? "↑" : "↓")
  sweeps = Sweeps(nsweeps)
  maxdims = min.(maxdim, [10, 20, 100, maxdim])
  maxdim!(sweeps, maxdims...)
  cutoff!(sweeps, 0.0)
  energy, ψ = dmrg(H, psi0, sweeps;
                   outputlevel = outputlevel)
  if outputlevel > 0
    @show nsweeps
    @show N
    @show cutoff
    @show energy
    @show flux(ψ)
  end
  return maxlinkdim(ψ)
end


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

# ––––––––
# Emphasis
# ––––––––

type Italic
    text
end

@trigger '*' ->
function asterisk_italic(stream::IO, md::MD)
    result = parse_inline_wrapper(stream, "*")
    return result == nothing ? nothing : Italic(parseinline(result, md))
end

type Bold
    text
end

@trigger '*' ->
function asterisk_bold(stream::IO, md::MD)
    result = parse_inline_wrapper(stream, "**")
    return result == nothing ? nothing : Bold(parseinline(result, md))
end

# ––––
# Code
# ––––

@trigger '`' ->
function inline_code(stream::IO, md::MD)
    result = parse_inline_wrapper(stream, "`"; rep=true)
    return result == nothing ? nothing : Code(result)
end

# ––––––––––––––
# Images & Links
# ––––––––––––––

type Image
    url::UTF8String
    alt::UTF8String
end

@trigger '!' ->
function image(stream::IO, md::MD)
    withstream(stream) do
        startswith(stream, "![") || return
        alt = readuntil(stream, ']', match = '[')
        alt ≡ nothing && return
        skipwhitespace(stream)
        startswith(stream, '(') || return
        url = readuntil(stream, ')', match = '(')
        url ≡ nothing && return
        return Image(url, alt)
    end
end

type Link
    text
    url::UTF8String
end

@trigger '[' ->
function link(stream::IO, md::MD)
    withstream(stream) do
        startswith(stream, '[') || return
        text = readuntil(stream, ']', match = '[')
        text ≡ nothing && return
        skipwhitespace(stream)
        startswith(stream, '(') || return
        url = readuntil(stream, ')', match = '(')
        url ≡ nothing && return
        return Link(parseinline(text, md), url)
    end
end

# –––––––––––
# Punctuation
# –––––––––––

type LineBreak end

@trigger '\\' ->
function linebreak(stream::IO, md::MD)
    if startswith(stream, "\\\n")
        return LineBreak()
    end
end

@trigger '-' ->
function en_dash(stream::IO, md::MD)
    if startswith(stream, "--")
        return "–"
    end
end

const escape_chars = "\\`*_#+-.!{[(\$"

@trigger '\\' ->
function escapes(stream::IO, md::MD)
    withstream(stream) do
        if startswith(stream, "\\") && !eof(stream) && (c = read(stream, Char)) in escape_chars
            return string(c)
        end
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ c769e6fe-59c2-11ec-1d28-23d959d24490
using CairoMakie

# ╔═╡ 91f6a6f4-274b-4f6f-81bc-e0313f78cdca
using XLSX

# ╔═╡ 1e4d24b2-a068-4d29-9061-c7196a2d1d1d
using OrdinaryDiffEq

# ╔═╡ a4ae4581-e7c7-4abb-a844-0e0bb34de766
using LsqFit

# ╔═╡ 4d2a70e0-6852-4125-a271-aba216ac5ed8
using DiffEqParamEstim

# ╔═╡ 2e61a3a8-8bb8-488f-be17-52e361a6c198
using BlackBoxOptim

# ╔═╡ 2d14ff52-5215-499f-aa4f-b4bcfaf1581f
const R = 8.31446261815324

# ╔═╡ a08e9b86-2d75-4eb0-b21e-f8e7fa0e45b8
getdatapaths(dirpath) = filter(readdir(dirpath; join=true)) do p
	return basename(p)[1:4] == "MK01" && split(p, '.')[end] == "xlsx"
end

# ╔═╡ 06c9e5d6-2e60-4713-8362-29761a55b929
datapaths = getdatapaths("../data")

# ╔═╡ 885005f0-60c7-4c72-b14c-2eafb4afb90a
function shiftzeros!(data)
	T = eltype(data)
	for i in eachindex(data)
		data[i] = max(eps(T), data[i])
	end

	return data
end

# ╔═╡ bc6c7f12-484a-4001-9b05-1a6f6c27efd5
function getdata(filepath)
    rawdata = XLSX.readtable(filepath, "Sheet1")

	tidx = findfirst(c -> c == :Time, rawdata[2])
	t = rawdata[1][tidx] |> Vector{Float64}
	t .-= first(t)  # Initial time set to 0

	Tidx = findfirst(c -> c == :Temperature_Cal, rawdata[2])
	T = rawdata[1][Tidx] |> Vector{Float64}
	T .+= 273.15  # to Kelvin

	αidx = findfirst(c -> c == :Conversion_DSC, rawdata[2])
	α = rawdata[1][αidx] |> Vector{Float64}
    return (T = T, α = α, t = t)
end

# ╔═╡ 05d49511-9430-494e-9713-336acef62525
filepath = datapaths[1]

# ╔═╡ 831a4e12-e2a8-4887-9d1c-1fa933adec92
pr=0.5

# ╔═╡ 8ba7c5cc-47ad-40b1-91e8-8bb4562b2f90
data = getdata(filepath)

# ╔═╡ 87d7f36d-c0f5-4e6e-99a6-822e72747dfd
function stripdata(data)
	idx1 = findfirst(x -> x > zero(x), data.α)
	idx2 = findlast(x -> x < one(x), data.α)

	(T = data.T[idx1:idx2], α = data.α[idx1:idx2], t = data.t[idx1:idx2])
end

# ╔═╡ e5d0e591-4bc2-4424-b30f-fe812a6f8a49
let fig = Figure()
	ax = Axis(fig[1, 1])

	lines!(ax, data.T, data.α)

	fig
end

# ╔═╡ 81f11d6e-b418-4b04-8653-4baea77a2087
function linfit(xs, ys)
    linmodel(x, p) = x .* p[1] .+ p[2]
	fit = curve_fit(linmodel, xs, ys, [1.0, 1.0])
	return (a = fit.param[1], b = fit.param[2], fit = fit)
end

# ╔═╡ b4c4af91-de7e-4a6a-bf20-54531ecd46f3
β, T0, Tfit = linfit(data.t, data.T)

# ╔═╡ d9b2b78c-052b-429f-a3b1-061eeeb7fd85
β * 60, T0 - 273.15

# ╔═╡ 216e2df2-8159-46b2-aca1-1b0eeb2b2eba
let fig = Figure()
	ax = Axis(fig[1, 1])

	lines!(ax, data.t, data.T, label="experiment")
	lines!(ax, data.t, β .* data.t .+ T0; label="fit")
	# lines!(ax, data.t, 4 / 60 .* data.t .+ (817.7268 + 273.15))

	axislegend(ax)
	fig
end

# ╔═╡ ac10a1f0-a6d2-4a35-9c5e-0e0e0e12634f
abstract type AbstractModel end

# ╔═╡ 417900e2-9fa0-483d-826a-e74d0b419b41
struct FxModel <: AbstractModel end

# ╔═╡ 58895114-1b00-41cf-a65c-6390e1c3d11d
function makeratefunction(β, T0, pr, ::FxModel)
	return function rate(α, p, t)
		A, B, Ea, ΔH, ΔS, C = p
		T = β * t + T0

		r1 = exp(A - Ea / R / T)
		r2 = 1 - pr * exp((ΔH / T - ΔS) / R)
		r3 = 1 - α
	return r1 * max(r2, 0)^B * max(r3, 0)^C
	end
end

# ╔═╡ 6c3e9753-a624-4d3f-a303-e21a45a39f9d
struct F1Model <: AbstractModel end

# ╔═╡ 035c8b21-0dbb-4652-b598-87ea8711d268
αrate(α, ::F1Model) = max(1 - α, 0)

# ╔═╡ b1cadb47-64d4-4d48-b9dc-e44d2cd7c43b
struct AE3Model <: AbstractModel end

# ╔═╡ f4a948fc-02ac-4111-8ed8-b0342d58a83f
αrate(α, ::AE3Model) = max(3 * (1 - α) * (-log(clamp(1 - α, eps(), 1 - eps())))^(2//3), 0)

# ╔═╡ 6eccb017-8abe-4c88-8afe-3d249f0ff0c0
struct PTModel <: AbstractModel end

# ╔═╡ 89723b6c-8ec7-4a16-9d1a-b93c1303f38e
αrate(α, ::PTModel) = (1 - α) * α

# ╔═╡ eb163472-c8b0-48af-95a0-cb9f6552fb48
struct DJ3Model <: AbstractModel end

# ╔═╡ 8ccb1799-c7f0-497d-8f67-bb3dc285839a
αrate(α, ::DJ3Model) = (c = cbrt(1 - α); 3/2 * c^2 / (1 - c))

# ╔═╡ 98d00c87-099a-4a3c-a5cf-e1300261f29b
struct GBModel <: AbstractModel end

# ╔═╡ 55760264-1ca1-43ed-a0ee-506d39098962
αrate(α, ::GBModel) = (c = cbrt(1 - α); 3/2 * c / (1 - c))

# ╔═╡ d73d9a7f-1f05-4355-9652-6b33bb27c78a
function makeratefunction(β, T0, pr, model::AbstractModel)
	return function rate(α, p, t)
		A, B, Ea, ΔH, ΔS = p
		T = β * t + T0

		r1 = exp(A - Ea / R / T)
		r2 = 1 - pr * exp((ΔH / T - ΔS) / R)
		r3 = αrate(α, model)
	return r1 * max(r2, 0)^B * r3
	end
end

# ╔═╡ 536a0f3d-70c3-4e7d-a136-d1262f491916
const MODELS = (FxModel(), F1Model(), AE3Model(), PTModel(), DJ3Model(), GBModel())

# ╔═╡ dbc79c4e-5ca2-47df-a7e4-19dc35149a18
makeratefunction(4, 500, 0.3, FxModel())(0.5, [20, 5, 2e5, 1.7e5, 150, 1.2], 456)

# ╔═╡ edab6163-e0f3-40bc-af8e-1d071ce9ec54
function fitkinetics(rate, ts, αs, p0, bounds, alg; losskwargs=(), solkwargs=(), optkwargs=())
	prob = ODEProblem(rate, first(αs), (first(ts), last(ts)), p0)
	loss = L2Loss(ts, αs)
	loss_objective = build_loss_objective(prob, alg, loss; solkwargs...)
	bboptimize(loss_objective, p0; SearchRange=bounds, optkwargs...)
end

# ╔═╡ 981584bd-3a65-4f99-83be-0095ffd61ea6
rate = makeratefunction(β, T0, pr, AE3Model())

# ╔═╡ fc1f8150-0b88-4d29-abb8-fceeb008bca0
p0 = [67.75974468855404, 10.891425746932146, 696081.8832423445, 131682.74431826902, 264.8105989317905]

# ╔═╡ 588951eb-d96e-4587-8eba-d78527aa2f97
let fig = Figure()
	ax = Axis(fig[1, 1])
	# ylims!(ax, -0.01, 0.015)

	dαs = [0; (data.α[2:end] .- data.α[1:end-1]) ./
		(data.t[2:end] .- data.t[1:end-1])]
	r0s = map(zip(data.α, data.t)) do (a, t)
		rate(a, p0, t)
	end

	lines!(ax, data.t, dαs; label="experiment")
	lines!(ax, data.t, r0s; label="model 0")

	axislegend(ax)

	fig
end

# ╔═╡ f432693b-0c8b-4f23-b4dc-99e2707fb9f8
res = fitkinetics(rate, data.t, data.α,
	p0, [(-50, 200), (-50, 50), (1e4, 2e6), (1e5, 2e5), (50, 300)],
	BS3(); losskwargs=(differ_weight=0,), solkwargs=(maxiters=1000, reltol=2e-3),
	optkwargs=(MaxTime=1, Method=:adaptive_de_rand_1_bin_radiuslimited)
)

# ╔═╡ 43f37494-6247-4b65-bdd4-2f5983981d9d
# C bounds = (0.0, 5.0)

# ╔═╡ b069488b-e485-46e1-966d-c1308514c4fa
bc = best_candidate(res)

# ╔═╡ e912818f-ee76-490d-934d-cd8cc2f52878
string(bc)

# ╔═╡ c6740e72-3870-4444-bdac-9d7529d37b36
best_fitness(res)

# ╔═╡ 2d4b3cbb-1f81-4962-bc64-55bb75e929aa
tspan = (first(data.t), last(data.t))

# ╔═╡ 58ad0dbe-94c7-4d31-b9fe-2c9ef99377d9
prob0 = ODEProblem(rate, first(data.α), tspan, p0)

# ╔═╡ 93fdc768-03cd-45ae-a1a7-69cfba984fdc
sol0 = solve(prob0, Vern7(); saveat=data.t,save_everystep=false,dense=false, reltol=1e-4);

# ╔═╡ db6e32dd-6c8a-4d01-a869-98d3d8d9a5a6
sum(abs2, sol0.(data.t) .- data.α)

# ╔═╡ 2b958eb1-e7d7-4abf-8537-34f5769d5482
bestprob = ODEProblem(rate, first(data.α), tspan, bc)

# ╔═╡ 1fa0c8a6-9e41-4840-a17b-8727445156b1
bestsol = solve(bestprob, Vern7(); saveat=data.t,save_everystep=false,dense=false, reltol=1e-4);

# ╔═╡ 7d5219c6-ad93-4dea-87e3-a092f2288e63
sum(abs2, bestsol.(data.t) .- data.α)

# ╔═╡ 7c0fe222-7413-4187-9bf2-e425f3af473c
let fig = Figure()
	ax = Axis(fig[1, 1])
	xlims!(ax, tspan)
	ylims!(ax, (-0.05, 1.05))

	xs = data.t
	ys = data.α

	lines!(ax, xs, ys; label="experiment")
	lines!(ax, xs, sol0.(xs); label="model p0")
	lines!(ax, xs, bestsol.(xs); label="best model")

	axislegend(ax)

	fig
end

# ╔═╡ 7fb59c6b-80cc-4e52-9cad-6a6ac8a635b9
let fig = Figure()
	ax = Axis(fig[1, 1])
	# ylims!(ax, -0.01, 0.015)

	dαs = [0; (data.α[2:end] .- data.α[1:end-1]) ./
		(data.t[2:end] .- data.t[1:end-1])]
	r0s = map(zip(data.α, data.t)) do (a, t)
		rate(a, p0, t)
	end
	brs = map(zip(data.α, data.t)) do (a, t)
		rate(a, bc, t)
	end

	lines!(ax, data.t, dαs; label="experiment")
	lines!(ax, data.t, r0s; label="model 0")
	lines!(ax, data.t, brs; label="best model")

	axislegend(ax)

	fig
end

# ╔═╡ 8dca68e6-b22a-47f0-b5c2-59936dbf40bd
let fig = Figure()
	ax = Axis(fig[1, 1])
	xlims!(ax, tspan)
	ylims!(ax, (-0.1, 0.1))

	xs = data.t
	ys = data.α

	lines!(ax, xs, zeros(length(xs)))
	lines!(ax, xs, sol0.(xs) .- ys)
	lines!(ax, xs, bestsol.(xs) .- ys)

	fig
end

# ╔═╡ 11045900-4a1f-430b-8412-f8dd63508818
const RATES = makeratefunction.(β, T0, 0.5, MODELS)

# ╔═╡ d7bf524a-0092-4775-bd9d-623e3a59019a
modeltostring(model::AbstractModel) = split(string(model), '.')[end]

# ╔═╡ 78d7ce10-9218-4a00-acce-40755a97328d
let fig = Figure()
	ax = Axis(fig[1, 1])
	ylims!(ax, -0.01, 0.02)

	lines!(ax, data.t, dαs; label="experiment")

	for model in MODELS
		rf = makeratefunction(β, T0, 0.5, model)
		rs = map(zip(data.α, data.t)) do (a, t)
			rf(a, p0, t)
		end
		lines!(ax, data.t, rs; label=modeltostring(model))
	end

	axislegend(ax)

	fig
end

# ╔═╡ Cell order:
# ╠═c769e6fe-59c2-11ec-1d28-23d959d24490
# ╠═91f6a6f4-274b-4f6f-81bc-e0313f78cdca
# ╠═1e4d24b2-a068-4d29-9061-c7196a2d1d1d
# ╠═a4ae4581-e7c7-4abb-a844-0e0bb34de766
# ╠═4d2a70e0-6852-4125-a271-aba216ac5ed8
# ╠═2e61a3a8-8bb8-488f-be17-52e361a6c198
# ╠═2d14ff52-5215-499f-aa4f-b4bcfaf1581f
# ╠═a08e9b86-2d75-4eb0-b21e-f8e7fa0e45b8
# ╠═06c9e5d6-2e60-4713-8362-29761a55b929
# ╠═885005f0-60c7-4c72-b14c-2eafb4afb90a
# ╠═bc6c7f12-484a-4001-9b05-1a6f6c27efd5
# ╠═05d49511-9430-494e-9713-336acef62525
# ╠═831a4e12-e2a8-4887-9d1c-1fa933adec92
# ╠═8ba7c5cc-47ad-40b1-91e8-8bb4562b2f90
# ╠═87d7f36d-c0f5-4e6e-99a6-822e72747dfd
# ╠═e5d0e591-4bc2-4424-b30f-fe812a6f8a49
# ╠═81f11d6e-b418-4b04-8653-4baea77a2087
# ╠═b4c4af91-de7e-4a6a-bf20-54531ecd46f3
# ╠═d9b2b78c-052b-429f-a3b1-061eeeb7fd85
# ╠═216e2df2-8159-46b2-aca1-1b0eeb2b2eba
# ╠═ac10a1f0-a6d2-4a35-9c5e-0e0e0e12634f
# ╠═417900e2-9fa0-483d-826a-e74d0b419b41
# ╠═58895114-1b00-41cf-a65c-6390e1c3d11d
# ╠═d73d9a7f-1f05-4355-9652-6b33bb27c78a
# ╠═6c3e9753-a624-4d3f-a303-e21a45a39f9d
# ╠═035c8b21-0dbb-4652-b598-87ea8711d268
# ╠═b1cadb47-64d4-4d48-b9dc-e44d2cd7c43b
# ╠═f4a948fc-02ac-4111-8ed8-b0342d58a83f
# ╠═6eccb017-8abe-4c88-8afe-3d249f0ff0c0
# ╠═89723b6c-8ec7-4a16-9d1a-b93c1303f38e
# ╠═eb163472-c8b0-48af-95a0-cb9f6552fb48
# ╠═8ccb1799-c7f0-497d-8f67-bb3dc285839a
# ╠═98d00c87-099a-4a3c-a5cf-e1300261f29b
# ╠═55760264-1ca1-43ed-a0ee-506d39098962
# ╠═536a0f3d-70c3-4e7d-a136-d1262f491916
# ╠═dbc79c4e-5ca2-47df-a7e4-19dc35149a18
# ╠═edab6163-e0f3-40bc-af8e-1d071ce9ec54
# ╠═981584bd-3a65-4f99-83be-0095ffd61ea6
# ╠═fc1f8150-0b88-4d29-abb8-fceeb008bca0
# ╟─588951eb-d96e-4587-8eba-d78527aa2f97
# ╠═f432693b-0c8b-4f23-b4dc-99e2707fb9f8
# ╠═43f37494-6247-4b65-bdd4-2f5983981d9d
# ╠═b069488b-e485-46e1-966d-c1308514c4fa
# ╠═e912818f-ee76-490d-934d-cd8cc2f52878
# ╠═c6740e72-3870-4444-bdac-9d7529d37b36
# ╠═2d4b3cbb-1f81-4962-bc64-55bb75e929aa
# ╠═58ad0dbe-94c7-4d31-b9fe-2c9ef99377d9
# ╠═93fdc768-03cd-45ae-a1a7-69cfba984fdc
# ╠═db6e32dd-6c8a-4d01-a869-98d3d8d9a5a6
# ╠═2b958eb1-e7d7-4abf-8537-34f5769d5482
# ╠═1fa0c8a6-9e41-4840-a17b-8727445156b1
# ╠═7d5219c6-ad93-4dea-87e3-a092f2288e63
# ╠═7c0fe222-7413-4187-9bf2-e425f3af473c
# ╟─7fb59c6b-80cc-4e52-9cad-6a6ac8a635b9
# ╠═8dca68e6-b22a-47f0-b5c2-59936dbf40bd
# ╠═11045900-4a1f-430b-8412-f8dd63508818
# ╠═d7bf524a-0092-4775-bd9d-623e3a59019a
# ╠═78d7ce10-9218-4a00-acce-40755a97328d

max_stars_repo_name: Luapulu/KineticFits

max_issues_repo_name: Luapulu/KineticFits

max_forks_repo_name: Luapulu/KineticFits


text:
module EMIRT
include("types.jl")
include("ios.jl")
include("Domains.jl")
include("common.jl")
include("sys.jl")
include("Images.jl")
include("Segmentations.jl")
include("AffinityMaps.jl")
include("evaluate.jl")
include("SegmentMSTs.jl")
include("parser.jl")

using .Types
using .Images
using .IOs
using .AffinityMaps
using .Segmentations
using .Evaluate
using .SegmentMSTs

export EMImage, Segmentation, AffinityMap, ParamDict, SegMST, SegmentPairs, SegmentPairAffities

end


text:
###############################################################################
#
# write messages to email file
#
###############################################################################
"""
Running program on the server mat take very long time, and one may not know the status without logging into the system. Thus, we provide a function to send
    email to notify the status:

$(METHODLIST)

"""
function send_email! end

"""
Here is the function to send email out:

    send_email!(subject::String,
                from_email::String,
                to_email::String,
                body::String)

Send out email, given
- `subject` Email subject
- `from_email` The outgoing email address
- `to_email` Email address to send out
- `body` Main body of the email

Note that you have to run this function on a machine with `sendmail` set up.

---
Example
```julia
send_email!("[DO NOT REPLY] Job status", "from.email", "to.email", "This is the main text");
```
"""
send_email!(subject::String,
            from_email::String,
            to_email::String,
            body::String) =
(
    @info tinfo("Sending email to $(to_email)...");

    # compose email and write the email to local drive
    _email_to_write = """
    Subject: $(subject)
    From: $(from_email)
    To: $(to_email)
    $(body)
    """;
    Base.write(".tmp_email", _email_to_write);

    # send the email through pipeline and delete the file
    run( pipeline(".tmp_email", `sendmail -t`) );
    rm(".tmp_email");

    return nothing
)


text:
module Biot

using LinearAlgebra, PyPlot, Interpolations

include("FieldCalc.jl")
include("Toroid_Inductance.jl")
include("MakeSolenoid.jl")
include("ToroidOptimizer.jl")
# Write your package code here.

end


text:
# TIPS:
# 1
const DEFAULT = 42

global x
y = f(x::Int + 1)

# 9
function with_keyword(x; name::String = "Smith")
   # ...
end

# 16
type Person
	name::String
	height::Float64
	weight::Float64
end
# instead of:
type Person
	name
	height
	weight
end

# 18
# Use
write(file, a, " ", b)
# instead of:
write(file, "$a $b")

# TOOLS:
# @time

# linter:
Pkg.add("Lint")
using Lint
lintfile("performance.jl")
#
# typechecker:
Pkg.add("TypeCheck")
using TypeCheck
foo1(x::Int) = isprime(x) ? x: false
checkreturntypes(foo1) #> foo1(Int64)::Union(Bool,Int64)

max_stars_repo_path: Module 1/Chapter09/Old/performance.jl

max_issues_repo_path: Module 1/Chapter09/Old/performance.jl

max_forks_repo_path: Module 1/Chapter09/Old/performance.jl


text:
# Soft threshold function
S(z, γ) = abs(z) <= γ ? zero(z) : ifelse(z > 0, z - γ, z + γ)

# Elastic net penalty with parameter α and given coefficients
function P{T}(α::T, β::SparseCoefficients{T})
    x = zero(T)
    @inbounds @simd for i = 1:nnz(β)
        x += (1 - α)/2*abs2(β.coef[i]) + α*abs(β.coef[i])
    end
    x
end

abstract CoordinateDescent{T,Intercept,M<:AbstractMatrix} <: LinPred
type NaiveCoordinateDescent{T,Intercept,M<:AbstractMatrix,S<:CoefficientIterator} <: CoordinateDescent{T,Intercept,M}
    X::M                          # original design matrix
    μy::T                         # mean of y at current weights
    μX::Vector{T}                 # mean of X at current weights (in predictor order)
    Xssq::Vector{T}               # weighted sum of squares of each column of X (in coefficient order)
    residuals::Vector{T}          # y - Xβ (unscaled with centered X)
    residualoffset::T             # offset of residuals (used only when X is sparse)
    weights::Vector{T}            # weights for each observation
    oldy::Vector{T}               # old y vector (for updating residuals without matrix multiplication)
    weightsum::T                  # sum(weights)
    coefitr::S                    # coefficient iterator
    dev::T                        # last deviance
    α::T                          # elastic net parameter
    maxiter::Int                  # maximum number of iterations
    maxncoef::Int                 # maximum number of coefficients
    tol::T                        # tolerance

    NaiveCoordinateDescent(X::M, α::T, maxncoef::Int, tol::T, coefitr::S) =
        new(X, zero(T), zeros(T, size(X, 2)), zeros(T, maxncoef), Array(T, size(X, 1)), zero(T),
            Array(T, size(X, 1)), Array(T, size(X, 1)), convert(T, NaN), coefitr, convert(T, NaN),
            α, typemax(Int), maxncoef, tol)
end

# Compute μX for all predictors
function computeμX!{T}(cd::CoordinateDescent{T})
    @extractfields cd X μX weights weightsum
    for ipred = 1:size(X, 2)
        μ = zero(T)
        @simd for i = 1:size(X, 1)
            @inbounds μ += X[i, ipred]*weights[i]
        end
        μX[ipred] = μ/weightsum
    end
    cd
end

function computeμX!{T,Intercept,M<:SparseMatrixCSC}(cd::CoordinateDescent{T,Intercept,M})
    @extractfields cd X μX weights weightsum
    @extractfields X rowval nzval colptr
    for ipred = 1:size(X, 2)
        μ = zero(T)
        @simd for i = colptr[ipred]:colptr[ipred+1]-1
            row = rowval[i]
            @inbounds μ += nzval[i]*weights[row]
        end
        μX[ipred] = μ/weightsum
    end
    cd
end

# Compute Xssq for given predictor, with mean subtracted if there is
# an intercept
function computeXssq{T,Intercept}(cd::NaiveCoordinateDescent{T,Intercept}, ipred::Int)
    @extractfields cd X weights
    μ = Intercept ? cd.μX[ipred] : zero(T)
    ssq = zero(T)
    @simd for i = 1:size(X, 1)
        @inbounds ssq += abs2(X[i, ipred] - μ)*weights[i]
    end
    ssq
end

function computeXssq{T,Intercept,M<:SparseMatrixCSC}(cd::NaiveCoordinateDescent{T,Intercept,M}, ipred::Int)
    @extractfields cd X weights
    @extractfields X rowval nzval colptr
    μ = Intercept ? cd.μX[ipred] : zero(T)
    ssq = zero(T)
    zeroweightsum = cd.weightsum
    @inbounds @simd for i = colptr[ipred]:colptr[ipred+1]-1
        row = rowval[i]
        ssq += abs2(nzval[i] - μ)*weights[row]
        zeroweightsum -= weights[row]
    end
    ssq + zeroweightsum*abs2(μ)
end

# Updates CoordinateDescent object with (possibly) new y vector and
# weights
function update!{T,Intercept}(cd::NaiveCoordinateDescent{T,Intercept}, coef::SparseCoefficients{T},
                              y::Vector{T}, wt::Vector{T})
    @extractfields cd residuals X Xssq weights oldy
    copy!(weights, wt)
    weightsum = cd.weightsum = sum(weights)
    weightsuminv = inv(weightsum)

    # Update residuals without recomputing X*coef
    if nnz(coef) == 0
        copy!(residuals, y)
        copy!(oldy, y)
    else
        cd.residualoffset = 0
        @inbounds @simd for i = 1:length(y)
            residuals[i] += y[i] - oldy[i]
            oldy[i] = y[i]
        end
    end

    if Intercept
        # Compute μy and μres
        μy = zero(T)
        μres = zero(T)
        @simd for i = 1:length(y)
            @inbounds μy += y[i]*weights[i]
            @inbounds μres += residuals[i]*weights[i]
        end
        μy *= weightsuminv
        μres *= weightsuminv
        cd.μy = μy

        # Center residuals
        @simd for i = 1:length(residuals)
            @inbounds residuals[i] -= μres
        end

        computeμX!(cd)
    end

    for icoef = 1:nnz(coef)
        Xssq[icoef] = computeXssq(cd, coef.coef2predictor[icoef])
    end
    cd
end

if VERSION < v"0.4.0-dev+707"
    # no-op inline macro for old Julia
    macro inline(x)
        esc(x)
    end
end

# Offset of each residual. This is used only for sparse matrices with
# an intercept, so that we can update only residuals for which a
# changed coefficient is non-zero instead of updating all of them. In
# the dense case, we need to update all coefficients anyway, so this
# strategy is unneeded.
residualoffset{T}(cd::NaiveCoordinateDescent{T}) = zero(T)
residualoffset{T,M<:SparseMatrixCSC}(cd::NaiveCoordinateDescent{T,true,M}) = cd.residualoffset

# Compute the gradient term (first term of RHS of eq. 8)
@inline function compute_grad{T}(::NaiveCoordinateDescent{T}, X::AbstractMatrix{T},
                                 residuals::Vector{T}, weights::Vector{T}, ipred::Int)
    v = zero(T)
    @simd for i = 1:size(X, 1)
        @inbounds v += X[i, ipred]*residuals[i]*weights[i]
    end
    v
end

@inline function compute_grad{T}(cd::NaiveCoordinateDescent{T}, X::SparseMatrixCSC{T},
                                 residuals::Vector{T}, weights::Vector{T}, ipred::Int)
    @extractfields X rowval nzval colptr
    @inbounds v = residualoffset(cd)*cd.weightsum*cd.μX[ipred]
    @inbounds @simd for i = colptr[ipred]:colptr[ipred+1]-1
        row = rowval[i]
        v += nzval[i]*residuals[row]*weights[row]
    end
    v
end

# Update coefficient and residuals, returning scaled squared difference
@inline function update_coef!{T,Intercept}(cd::NaiveCoordinateDescent{T,Intercept},
                                           coef::SparseCoefficients{T},
                                           newcoef::T, icoef::Int, ipred::Int)
    coefdiff = coef.coef[icoef] - newcoef
    if coefdiff != 0
        coef.coef[icoef] = newcoef
        μ = Intercept ? cd.μX[ipred] : zero(T)

        @extractfields cd X residuals
        @simd for i = 1:size(X, 1)
            @inbounds residuals[i] += coefdiff*(X[i, ipred] - μ)
        end
        abs2(coefdiff)*cd.Xssq[icoef]
    else
        zero(T)
    end
end

@inline function update_coef!{T,Intercept,M<:SparseMatrixCSC}(cd::NaiveCoordinateDescent{T,Intercept,M},
                                                              coef::SparseCoefficients{T},
                                                              newcoef::T, icoef::Int, ipred::Int)
    coefdiff = coef.coef[icoef] - newcoef
    if coefdiff != 0
        coef.coef[icoef] = newcoef
        @extractfields cd X residuals
        @extractfields X rowval nzval colptr
        if Intercept
            cd.residualoffset -= coefdiff*cd.μX[ipred]
        end
        @inbounds @simd for i = colptr[ipred]:colptr[ipred+1]-1
            row = rowval[i]
            residuals[row] += coefdiff*nzval[i]
        end
        abs2(coefdiff)*cd.Xssq[icoef]
    else
        zero(T)
    end
end

# Performs the cycle of all predictors
function cycle!{T}(coef::SparseCoefficients{T}, cd::NaiveCoordinateDescent{T}, λ::T, all::Bool)
    @extractfields cd residuals X weights Xssq α

    maxdelta = zero(T)
    @inbounds if all
        # Use all predictors for first and last iterations
        for ipred = 1:size(X, 2)
            v = compute_grad(cd, X, residuals, weights, ipred)

            icoef = coef.predictor2coef[ipred]
            if icoef != 0
                oldcoef = coef.coef[icoef]
                v += Xssq[icoef]*oldcoef
            else
                # Adding a new variable to the model
                abs(v) < λ*α && continue
                oldcoef = zero(T)
                nnz(coef) > cd.maxncoef &&
                    error("maximum number of coefficients $(cd.maxncoef) exceeded at λ = $λ")
                icoef = addcoef!(coef, ipred)
                cd.coefitr = addcoef(cd.coefitr, icoef)
                Xssq[icoef] = computeXssq(cd, ipred)
            end
            newcoef = S(v, λ*α)/(Xssq[icoef] + λ*(1 - α))

            maxdelta = max(maxdelta, update_coef!(cd, coef, newcoef, icoef, ipred))
        end
    else
        # Iterate over only the predictors already in the model
        for icoef = cd.coefitr
            oldcoef = coef.coef[icoef]
            oldcoef == 0 && continue
            ipred = coef.coef2predictor[icoef]

            v = Xssq[icoef]*oldcoef + compute_grad(cd, X, residuals, weights, ipred)
            newcoef = S(v, λ*α)/(Xssq[icoef] + λ*(1 - α))

            maxdelta = max(maxdelta, update_coef!(cd, coef, newcoef, icoef, ipred))
        end
    end
    maxdelta
end

# Sum of squared residuals. Residuals are always up to date
function ssr{T}(coef::SparseCoefficients{T}, cd::NaiveCoordinateDescent{T})
    @extractfields cd residuals weights
    roffset = residualoffset(cd)
    s = zero(T)
    @simd for i = 1:length(residuals)
        @inbounds s += abs2(residuals[i] + roffset)*weights[i]
    end
    s
end

# Value of the intercept
intercept{T}(coef::SparseCoefficients{T}, cd::CoordinateDescent{T,false}) = zero(T)
function intercept{T}(coef::SparseCoefficients{T}, cd::NaiveCoordinateDescent{T,true})
    μX = cd.μX
    v = cd.μy
    for icoef = 1:nnz(coef)
        v -= μX[coef.coef2predictor[icoef]]*coef.coef[icoef]
    end
    v
end

# Value of the linear predictor
function linpred!{T}(mu::Vector{T}, cd::NaiveCoordinateDescent{T}, coef::SparseCoefficients{T}, b0::T)
    @extractfields cd oldy residuals
    roffset = residualoffset(cd)
    @simd for i = 1:length(mu)
        @inbounds mu[i] = oldy[i] - (residuals[i] + roffset)
    end
    mu
end

type CovarianceCoordinateDescent{T,Intercept,M<:AbstractMatrix,S<:CoefficientIterator} <: CoordinateDescent{T,Intercept,M}
    X::M                          # original design matrix
    μy::T                         # mean of y at current weights
    μX::Vector{T}                 # mean of X at current weights
    yty::T                        # y'y (scaled by weights)
    Xty::Vector{T}                # X'y (scaled by weights)
    Xssq::Vector{T}               # weighted sum of squares of each column of X
    XtX::Matrix{T}                # X'X (scaled by weights, in order of coefficients)
    tmp::Vector{T}                # scratch used when computing X'X
    weights::Vector{T}            # weights for each observation
    weightsum::T                  # sum(weights)
    coefitr::S                    # coefficient iterator
    dev::T                        # last deviance
    α::T                          # elastic net parameter
    maxiter::Int                  # maximum number of iterations
    maxncoef::Int                 # maximum number of coefficients
    tol::T                        # tolerance

    function CovarianceCoordinateDescent(X::M, α::T, maxncoef::Int, tol::T, coefiter::S)
        new(X, zero(T), zeros(T, size(X, 2)), convert(T, NaN), Array(T, size(X, 2)),
            Array(T, size(X, 2)), Array(T, maxncoef, size(X, 2)), Array(T, size(X, 1)),
            Array(T, size(X, 1)), convert(T, NaN), coefiter, convert(T, NaN), α,
            typemax(Int), maxncoef, tol)
    end
end

# Compute y'y
function computeyty{T,Intercept,M}(cd::CovarianceCoordinateDescent{T,Intercept,M}, y::Vector{T})
    @extractfields cd X weights
    μy = Intercept ? cd.μy : zero(T)
    yty = zero(T)
    @inbounds @simd for i = 1:length(y)
        yty += abs2(y[i] - μy)*weights[i]
    end
    yty
end

# Compute Xssq and Xty for given predictor, with mean subtracted if
# there is an intercept.
function computeXssqXty{T,Intercept}(cd::CovarianceCoordinateDescent{T,Intercept}, y::Vector{T}, ipred::Int)
    @extractfields cd X μy weights
    μ = Intercept ? cd.μX[ipred] : zero(T)
    ssq = zero(T)
    ty = zero(T)
    @inbounds @simd for i = 1:size(X, 1)
        ssq += abs2(X[i, ipred] - μ)*weights[i]
        ty += (y[i] - μy)*(X[i, ipred] - μ)*weights[i]
    end
    (ssq, ty)
end

function computeXssqXty{T,Intercept,M<:SparseMatrixCSC}(cd::CovarianceCoordinateDescent{T,Intercept,M}, y::Vector{T}, ipred::Int)
    @extractfields cd X μy weights weightsum
    @extractfields X rowval nzval colptr
    μ = Intercept ? cd.μX[ipred] : zero(T)
    ssq = zero(T)
    ty = zero(T)

    zeroweightsum = weightsum
    zeroweighty = zero(T)
    @inbounds @simd for i = colptr[ipred]:colptr[ipred+1]-1
        row = rowval[i]
        ssq += abs2(nzval[i] - μ)*weights[row]
        ty += (y[row] - μy)*weights[row]*(nzval[i] - μ)
        zeroweightsum -= weights[row]
        zeroweighty -= (y[row] - μy)*weights[row]
    end

    # Correct for zero values
    ssq += abs2(μ)*zeroweightsum
    ty -= μ*zeroweighty

    (ssq, ty)
end

# Compute:
# XtX = (X .- μX')'*(weights.*(X[:, icoef] - μX[icoef]))
#     = X'*(weights.*(X[:, icoef] - μX[icoef])) - μX*weights'*(X[:, icoef] - μX[icoef])
#     = X'*(weights.*(X[:, icoef] - μX[icoef]))
# since μX'*weights'*(X[:, icoef] - μX[icoef])
#     = \sum_{i=1}^n w*(X_{i, icoef} - \frac{\sum_{j=1}^n X_{j, icoef} w_j}{\sum_{j=1}^n w})
#     = 0
#
# This must be called for all icoef < icoef before it may be called
# for a given icoef
function computeXtX!{T,Intercept}(cd::CovarianceCoordinateDescent{T,Intercept},
                                  coef::SparseCoefficients{T}, icoef::Int, ipred::Int)
    @extractfields cd X tmp weights XtX μX

    @simd for i = 1:size(X, 1)
        @inbounds tmp[i] = X[i, ipred]*weights[i]
    end

    @inbounds for jpred = 1:size(X, 2)
        if jpred == ipred
            XtX[icoef, jpred] = cd.Xssq[jpred]
        else
            jcoef = coef.predictor2coef[jpred]
            if 0 < jcoef < icoef
                XtX[icoef, jpred] = XtX[jcoef, ipred]
            else
                s = zero(T)
                μ = Intercept ? μX[jpred] : zero(T)
                @simd for i = 1:size(X, 1)
                    s += tmp[i]*(X[i, jpred] - μ)
                end
                XtX[icoef, jpred] = s
            end
        end
    end
    cd
end

# In the sparse case, we actually compute:
# XtX = X'*(weights.*X[:, icoef]) - μX*μX[icoef]*sum(weights)
#
# TODO: Can this be done efficiently without making X[:, icoef] dense?
# TODO: Can this be done efficiently with a sparse XtX, e.g by
#       saving only the first term above? Storage may be a problem.
function computeXtX!{T,Intercept,M<:SparseMatrixCSC}(cd::CovarianceCoordinateDescent{T,Intercept,M},
                                                     coef::SparseCoefficients{T}, icoef::Int, ipred::Int)
    @extractfields cd X tmp weights weightsum XtX μX
    @extractfields X colptr rowval nzval

    fill!(tmp, zero(T))
    μi = Intercept ? μX[ipred] : zero(T)
    @inbounds @simd for i = colptr[ipred]:colptr[ipred+1]-1
        tmp[rowval[i]] = nzval[i]*weights[rowval[i]]
    end

    @inbounds for jpred = 1:size(X, 2)
        if jpred == ipred
            XtX[icoef, jpred] = cd.Xssq[jpred] # + abs2(μi)*weightsum
        else
            jcoef = coef.predictor2coef[jpred]
            if 0 < jcoef < icoef
                XtX[icoef, jpred] = XtX[jcoef, ipred]
            else
                s = zero(T)
                @simd for i = colptr[jpred]:colptr[jpred+1]-1
                    @inbounds s += tmp[rowval[i]]*nzval[i]
                end
                μj = Intercept ? μX[jpred] : zero(T)
                XtX[icoef, jpred] = s - μi*μj*weightsum
            end
        end
    end
    cd
end

# Updates CoordinateDescent object with (possibly) new y vector and
# weights
function update!{T,Intercept,M}(cd::CovarianceCoordinateDescent{T,Intercept,M},
                                coef::SparseCoefficients{T}, y::Vector{T}, wt::Vector{T})
    @extractfields cd X Xty μX Xssq weights

    copy!(weights, wt)
    weightsum = cd.weightsum = sum(weights)
    weightsuminv = inv(weightsum)

    # Compute μy
    μy = zero(T)
    if Intercept
        @simd for i = 1:length(y)
            @inbounds μy += y[i]*weights[i]
        end
        μy *= weightsuminv
        cd.μy = μy
    end

    cd.yty = computeyty(cd, y)

    if Intercept
        # Compute weighted mean
        computeμX!(cd)
    end

    # Compute Xssq and X'y
    for j = 1:size(X, 2)
        Xssq[j], Xty[j] = computeXssqXty(cd, y, j)
    end

    # Compute X'X for predictors in model
    for icoef = 1:nnz(coef)
        computeXtX!(cd, coef, icoef, coef.coef2predictor[icoef])
    end

    cd
end

# This is abstracted out in case someday we want to make XtX sparse
getXtX(::CovarianceCoordinateDescent, XtX, jcoef, ipred) = XtX[jcoef, ipred]

function compute_gradient{T}(cd::CovarianceCoordinateDescent{T}, XtX, coef, ipred)
    s = zero(T)
    @simd for jcoef = 1:nnz(coef)
        @inbounds s += getXtX(cd, XtX, jcoef, ipred)*coef.coef[jcoef]
    end
    s
end

# Performs the cycle of all predictors
function cycle!{T}(coef::SparseCoefficients{T}, cd::CovarianceCoordinateDescent{T}, λ::T, all::Bool)
    @extractfields cd X Xty XtX Xssq α

    maxdelta = zero(T)
    if all
        @inbounds for ipred = 1:length(Xty)
            # Use all predictors for first and last iterations
            s = Xty[ipred] - compute_gradient(cd, XtX, coef, ipred)

            icoef = coef.predictor2coef[ipred]
            if icoef != 0
                oldcoef = coef.coef[icoef]
                s += getXtX(cd, XtX, icoef, ipred)*oldcoef
            else
                oldcoef = zero(T)
            end

            newcoef = S(s, λ*α)/(Xssq[ipred] + λ*(1 - α))
            if oldcoef != newcoef
                if icoef == 0
                    # Adding a new variable to the model
                    nnz(coef) > cd.maxncoef &&
                        error("maximum number of coefficients $(cd.maxncoef) exceeded at λ = $λ")
                    icoef = addcoef!(coef, ipred)
                    cd.coefitr = addcoef(cd.coefitr, icoef)

                    # Compute cross-product with predictors
                    computeXtX!(cd, coef, icoef, ipred)
                end
                maxdelta = max(maxdelta, abs2(oldcoef - newcoef)*Xssq[ipred])
                coef.coef[icoef] = newcoef
            end
        end
    else
        # Iterate over only the predictors already in the model
        @inbounds for icoef = cd.coefitr
            ipred = coef.coef2predictor[icoef]
            oldcoef = coef.coef[icoef]
            oldcoef == 0 && continue
            s = Xty[ipred] + getXtX(cd, XtX, icoef, ipred)*oldcoef - compute_gradient(cd, XtX, coef, ipred)
            newcoef = coef.coef[icoef] = S(s, λ*α)/(Xssq[ipred] + λ*(1 - α))
            maxdelta = max(maxdelta, abs2(oldcoef - newcoef)*Xssq[ipred])
        end
    end
    maxdelta
end

# y'y - 2β'X'y + β'X'Xβ
function ssr{T}(coef::SparseCoefficients{T}, cd::CovarianceCoordinateDescent{T})
    XtX = cd.XtX
    Xty = cd.Xty
    v = cd.yty
    @inbounds for icoef = 1:nnz(coef)
        ipred = coef.coef2predictor[icoef]
        s = -2*Xty[ipred] + compute_gradient(cd, XtX, coef, ipred)
        v += coef.coef[icoef]*s
    end
    v
end

# Value of the intercept
intercept{T}(coef::SparseCoefficients{T}, cd::CovarianceCoordinateDescent{T,true}) =
    cd.μy - dot(cd.μX, coef)

# Value of the linear predictor
function linpred!{T}(mu::Vector{T}, cd::CovarianceCoordinateDescent{T},
                     coef::SparseCoefficients{T}, b0::T)
    A_mul_B!(mu, cd.X, coef)
    if b0 != 0
        @simd for i = 1:length(mu)
            @inbounds mu[i] += b0
        end
    end
    mu
end

function cdfit!{T}(coef::SparseCoefficients{T}, cd::CoordinateDescent{T}, λ, criterion)
    maxiter = cd.maxiter
    tol = cd.tol
    n = size(cd.X, 1)

    obj = convert(T, Inf)
    objold = convert(T, Inf)
    dev = convert(T, Inf)
    prev_converged = false
    converged = true
    b0 = intercept(coef, cd)

    iter = 0
    for iter = 1:maxiter
        oldb0 = b0
        maxdelta = cycle!(coef, cd, λ, converged)
        b0 = intercept(coef, cd)
        maxdelta = max(maxdelta, abs2(oldb0 - b0)*cd.weightsum)

        # Test for convergence
        prev_converged = converged
        if criterion == :obj
            objold = obj
            dev = ssr(coef, cd)
            obj = dev/2 + λ*P(α, coef)
            converged = objold - obj < tol*obj
        elseif criterion == :coef
            converged = maxdelta < tol
        end

        # Require two converging steps to return. After the first, we
        # will iterate through all variables
        prev_converged && converged && break
    end

    (!prev_converged || !converged) &&
        error("coordinate descent failed to converge in $maxiter iterations at λ = $λ")

    cd.dev = criterion == :coef ? ssr(coef, cd) : dev
    iter
end

# Compute working response based on eta and working residuals
function wrkresp!(out, eta, wrkresid, offset)
    if isempty(offset)
        @simd for i = 1:length(out)
            @inbounds out[i] = eta[i] + wrkresid[i]
        end
    else
        @simd for i = 1:length(out)
            @inbounds out[i] = eta[i] + wrkresid[i] - offset[i]
        end
    end
    out
end

# Fits GLMs (outer and middle loops)
function StatsBase.fit!{S<:GeneralizedLinearModel,T}(path::LassoPath{S,T}; verbose::Bool=false, irls_maxiter::Int=30,
                                                     cd_maxiter::Int=100000, cd_tol::Real=1e-7, irls_tol::Real=1e-7,
                                                     criterion=:coef, minStepFac::Real=0.001)
    irls_maxiter >= 1 || error("irls_maxiter must be positive")
    0 < minStepFac < 1 || error("minStepFac must be in (0, 1)")
    criterion == :obj || criterion == :coef || error("criterion must be obj or coef")

    @extractfields path nulldev λ autoλ Xnorm m
    nλ = length(λ)
    m = path.m
    r = m.rr
    cd = m.pp
    cd.maxiter = cd_maxiter
    cd.tol = cd_tol

    if criterion == :coef
        cd_tol *= path.nulldev/2
        irls_tol *= path.nulldev/2
    end

    @extractfields cd X α
    @extractfields r offset eta wrkresid
    coefs = spzeros(T, size(X, 2), nλ)
    b0s = zeros(T, nλ)
    oldcoef = SparseCoefficients{T}(size(X, 2))
    newcoef = SparseCoefficients{T}(size(X, 2))
    pct_dev = zeros(nλ)
    dev_ratio = convert(T, NaN)
    dev = convert(T, NaN)
    b0 = zero(T)
    scratchmu = Array(T, size(X, 1))
    objold = convert(T, Inf)

    if autoλ
        # No need to fit the first model
        coefs[:, 1] = zero(T)
        b0s[1] = path.nullb0
        i = 2
    else
        i = 1
    end

    dev = NaN
    niter = 0
    if nλ == 0
        i = 0
    else
        while true # outer loop
            obj = convert(T, Inf)
            last_dev_ratio = dev_ratio
            curλ = λ[i]

            converged = false

            for iirls=1:irls_maxiter # middle loop
                copy!(oldcoef, newcoef)
                oldb0 = b0

                # Compute working response
                wrkresp!(scratchmu, eta, wrkresid, offset)
                wrkwt = wrkwt!(r)

                # Run coordinate descent inner loop
                niter += cdfit!(newcoef, update!(cd, newcoef, scratchmu, wrkwt), curλ, criterion)
                b0 = intercept(newcoef, cd)

                # Update GLM and get deviance
                updatemu!(r, linpred!(scratchmu, cd, newcoef, b0))

                # Compute Elastic Net objective
                objold = obj
                dev = deviance(r)
                obj = dev/2 + curλ*P(α, newcoef)

                if obj > objold + length(scratchmu)*eps(objold)
                    f = 1.0
                    b0diff = b0 - oldb0
                    while obj > objold
                        f /= 2.; f > minStepFac || error("step-halving failed at beta = $(newcoef)")
                        for icoef = 1:nnz(newcoef)
                            oldcoefval = icoef > nnz(oldcoef) ? zero(T) : oldcoef.coef[icoef]
                            newcoef.coef[icoef] = oldcoefval+f*(newcoef.coef[icoef] - oldcoefval)
                        end
                        b0 = oldb0+f*b0diff
                        updatemu!(r, linpred!(scratchmu, cd, newcoef, b0))
                        dev = deviance(r)
                        obj = dev/2 + curλ*P(α, newcoef)
                    end
                end

                # Determine if we have converged
                if criterion == :obj
                    converged = objold - obj < irls_tol*obj
                elseif criterion == :coef
                    maxdelta = zero(T)
                    Xssq = cd.Xssq
                    converged = abs2(oldb0 - b0)*cd.weightsum < irls_tol
                    if converged
                        for icoef = 1:nnz(newcoef)
                            oldcoefval = icoef > nnz(oldcoef) ? zero(T) : oldcoef.coef[icoef]
                            j = isa(cd, NaiveCoordinateDescent) ? icoef : newcoef.coef2predictor[icoef]
                            if abs2(oldcoefval - newcoef.coef[icoef])*Xssq[j] > irls_tol
                                converged = false
                                break
                            end
                        end
                    end
                end

                converged && break
            end
            converged || error("IRLS failed to converge in $irls_maxiter iterations at λ = $(curλ)")

            dev_ratio = dev/nulldev

            pct_dev[i] = 1 - dev_ratio
            addcoefs!(coefs, newcoef, i)
            b0s[i] = b0

            # Test whether we should continue
            if i == nλ || (autoλ && last_dev_ratio - dev_ratio < MIN_DEV_FRAC_DIFF ||
                           pct_dev[i] > MAX_DEV_FRAC)
                break
            end

            i += 1
        end
    end

    path.λ = path.λ[1:i]
    path.pct_dev = pct_dev[1:i]
    path.coefs = coefs[:, 1:i]
    path.b0 = b0s[1:i]
    path.niter = niter
    if !isempty(Xnorm)
        scale!(Xnorm, path.coefs)
    end
end

# Fits linear models (just the outer loop)
function StatsBase.fit!{S<:LinearModel,T}(path::LassoPath{S,T}; verbose::Bool=false,
                                          cd_maxiter::Int=10000, cd_tol::Real=1e-7, irls_tol::Real=1e-7,
                                          criterion=:coef, minStepFac::Real=eps())
    0 < minStepFac < 1 || error("minStepFac must be in (0, 1)")
    criterion == :obj || criterion == :coef || error("criterion must be obj or coef")

    @extractfields path nulldev λ autoλ Xnorm m
    nλ = length(λ)
    r = m.rr
    cd = m.pp
    cd.maxiter = cd_maxiter
    cd.tol = cd_tol

    X = cd.X
    coefs = spzeros(T, size(X, 2), nλ)
    b0s = zeros(T, nλ)
    newcoef = SparseCoefficients{T}(size(X, 2))
    pct_dev = zeros(nλ)
    dev_ratio = convert(T, NaN)
    niter = 0

    update!(cd, newcoef, r.mu, r.wts)

    if autoλ
        # No need to fit the first model
        b0s[1] = path.nullb0
        i = 2
    else
        i = 1
    end

    while true # outer loop
        last_dev_ratio = dev_ratio
        curλ = λ[i]

        # Run coordinate descent
        niter += cdfit!(newcoef, cd, curλ, criterion)

        dev_ratio = cd.dev/nulldev
        pct_dev[i] = 1 - dev_ratio
        addcoefs!(coefs, newcoef, i)
        b0s[i] = intercept(newcoef, cd)

        # Test whether we should continue
        if i == nλ || (autoλ && last_dev_ratio - dev_ratio < MIN_DEV_FRAC_DIFF ||
                       pct_dev[i] > MAX_DEV_FRAC)
            break
        end

        i += 1
    end

    path.λ = path.λ[1:i]
    path.pct_dev = pct_dev[1:i]
    path.coefs = coefs[:, 1:i]
    path.b0 = b0s[1:i]
    path.niter = niter
    if !isempty(Xnorm)
        scale!(Xnorm, path.coefs)
    end
end


text:
module Render
using Images
using ..Grids
using ImageView
using OffsetArrays

export to_img, pack_imgs, to_img_diff, to_img_diffgrid

ImageView.imshow(io::Union{ARCIO,ARCTask,ARCDiff,ARCDiffGrid,ARCGrid}) = imshow(to_img(io))

@inline scale_grid(img::AbstractMatrix,scale::Int) = repeat(img; inner=(scale,scale))

# surprisingly verbose function to scale from a grid pixel location to a view of a range of pixels on a larger grid
@inline function view_px(img::AbstractArray, scale::Int, x::Int, y::Int)
    top_left_x = (x-1)*scale + 1
    top_left_y = (y-1)*scale + 1
    view(img,top_left_x:top_left_x+scale-1,top_left_y:top_left_y+scale-1)
end
@inline view_px(img::AbstractArray, scale::Int, I::CartesianIndex) = view_px(img,scale,I[1],I[2])

# function gridlines!(img,cell_sz; color=colorant"gray")
#     @assert mod.(size(img),cell_sz) == (0,0)
#     img[1:cell_sz:end,:] .= color
#     img[:,1:cell_sz:end] .= color
#     img[end,:] .= color
#     img[:,end] .= color
#     img
# end

function borders!(img, up=true, down=true, left=true, right=true; color=colorant"gray")
    up && (img[1,:] .= color)
    down && (img[end,:] .= color)
    left && (img[:,1] .= color)
    right && (img[:,end] .= color)
    img
end

function to_img(task::ARCTask; px_sz=20)
    pack_imgs(
    pack_imgs(dim=1,to_img.(task.ios,px_sz=px_sz)...),
    pack_imgs(dim=1,to_img.(task.test_ios,px_sz=px_sz)...)
    )
end

function to_img(io::ARCIO; px_sz=20)
    pack_imgs(to_img(io.i,px_sz=px_sz),to_img(io.o,px_sz=px_sz))
end

function to_img(grid::ARCGrid; px_sz=20)
    img = map(color->colors[color+1], grid)
    img = scale_grid(img,px_sz)
    if px_sz > 3
        for I in CartesianIndices(grid)
            borders!(view_px(img,px_sz,I))
        end
    end
    img
end

function to_img(diff::ARCDiff; edges=true)

    function grayscale(img)
        img /= 9 # scale to 0.0-1.0
        img ./= 2 # scale to 0.0-0.5
        img[img.>0] .+= .3 # extra boost for nonzero (so black background stays dark)
        img .+ .2 # everyone gets a slight boost
    end

    # R G B are all indexed like A but B has some negative indices in it too
    R = grayscale(diff.A)
    G = grayscale(diff.B)
    B = Float64.(diff.match)

    img = collect(colorview(RGB, paddedviews(0, R, G, B)...))

    if edges
        img = scale_grid(img,5)
        for IA in CartesianIndices(diff.A)
            borders!(view_px(img,5,IA+diff.padded_offsets),diff.edges_match[IA,:]..., color=colorant"white")
        end
    end

    parent(img) # remove offsets for rendering
end

function to_img(diff_grid::ARCDiffGrid; edges=true)
    img_grid = to_img.(diff_grid.grid,edges=edges)
    img_grid = pack_imgs(
        [pack_imgs(imgs...,dim=2) for imgs in eachrow(img_grid)]...,
        dim=1
        )
    pack_imgs(to_img(diff_grid.A),to_img(diff_grid.B),img_grid,dim=2)
end

function pack_imgs(imgs::Matrix...; dim=2, pad=colorant"light gray", npad=5)
    @assert dim in [1,2]
    if dim == 2
        dim1 = maximum(size.(imgs,1))
        dim2 = sum(size.(imgs,2)) + (npad*(length(imgs)-1))
    else
        dim2 = maximum(size.(imgs,2))
        dim1 = sum(size.(imgs,1)) + (npad*(length(imgs)-1))
    end
    out = fill(pad, dim1, dim2)
    offset = 1
    if dim == 2
        for img in imgs
            d1,d2 = size(img)
            out[1:d1, offset:offset + d2 - 1] = img
            offset += d2 + npad
        end
    else
        for img in imgs
            d1,d2 = size(img)
            out[offset:offset + d1 - 1, 1:d2] = img
            offset += d1 + npad
        end
    end
    out
end

end


text:
using Compat
using Grabbit
using Grabbit: parse_config
using Base.Test

# write your own tests here

@testset "Grabbit" begin
    config_fn = joinpath(Pkg.dir("Grabbit"), "test", "specs", "test.json")
    root = joinpath(Pkg.dir("Grabbit"), "test", "data", "7t_trt")

    config = parse_config(root, config_fn)

    @testset "Domain" begin
        d = Domain(config)
    end

    @testset "Layout" begin
        layout = Layout(config)

        session1files = get(layout, session=1)
        @test all(haskey(f.tags, "session") for f in session1files)
        @test all(f.tags["session"] == "1" for f in session1files)
        @test unique(f.tags["subject"] for f in session1files) ==
            [@sprintf("%02d", i) for i in 1:10]
        @test any(!haskey(f.tags, "acquisition") for f in session1files)

        # use second arg Dict with queries equivalent to kw args
        @test get(layout, Dict("session"=>1)) == session1files

        subj12 = get(layout, subject = [1,2])
        @test all(haskey(f.tags, "subject") for f in subj12)
        @test unique(f.tags["subject"] for f in subj12) == ["01", "02"]

        @test isempty(setdiff(subj12, vcat(get(layout, subject=1), get(layout, subject=2))))

        @test subj12 == get(layout, subject = ["01", "02"])
        @test subj12 == get(layout, subject = [1, "02"])
        @test subj12 == get(layout, subject = [2, 1])
        @test subj12 == get(layout, subject = [1, 2, 1])

        fs = get(layout, Dict("type"=>["magnitude1", "bold"]), subject="01", session=1)
        @test isempty(setdiff(unique(get.(fs, "type")), ["magnitude1", "bold"]))
        @test all(get.(fs, "subject") .== "01")
        @test all(get.(fs, "session") .== "1")
        @test all(get.(fs, "extension") .== "nii.gz")
    end

    @testset "Default entities" begin
        @testset "extension" begin
            layout = Layout(config)
            @test !isempty(get(layout, extension="json"))
            @test length(get(layout, subject=1, session=1, run=1, extension="json")) == 1
            @test length(get(layout, subject=1, extension="nii.gz")) == 20
            # needs to be exact match
            @test isempty(get(layout, extension="nii"))
            # ...unless provided a regex
            @test length(get(layout, subject=1, extension="nii.*")) == 20
            @test isempty(get(layout, extension=["txt", "rtf"]))
            @test length(get(layout, subject=1, extension="json")) == 4
            @test length(get(layout, subject=1, extension=["nii.gz", "json"])) == 24
        end
    end

    @testset "excludes/includes" begin
        layout = Layout(root, config_fn)
        @test !any(occursin.(r"derivatives", dirname.(get(layout))))

        config_incl = parse_config(root,
                                   joinpath(Pkg.dir("Grabbit"),
                                            "test",
                                            "specs",
                                            "test_include.json"))
        layout_incl = Layout(config_incl)
        @test !any(occursin.(r"models/excluded_model.json", path.(get(layout_incl))))
    end

end


text:
# Keep these in sync with the Symbols in src/core/definitions.
get_store_container_type(::AuxVarKey) = STORE_CONTAINER_AUX_VARIABLES
get_store_container_type(::ConstraintKey) = STORE_CONTAINER_DUALS
get_store_container_type(::ExpressionKey) = STORE_CONTAINER_EXPRESSIONS
get_store_container_type(::ParameterKey) = STORE_CONTAINER_PARAMETERS
get_store_container_type(::VariableKey) = STORE_CONTAINER_VARIABLES

# Aliases used for clarity in the method dispatches so it is possible to know if writing to
# DecisionModel data or EmulationModel data
const DecisionModelIndexType = Dates.DateTime
const EmulationModelIndexType = Int

function write_results!(
    store,
    model::OperationModel,
    index::Union{DecisionModelIndexType, EmulationModelIndexType},
    update_timestamp::Dates.DateTime;
    exports=nothing,
)
    if exports !== nothing
        export_params = Dict{Symbol, Any}(
            :exports => exports,
            :exports_path => joinpath(exports.path, string(get_name(model))),
            :file_type => get_export_file_type(exports),
            :resolution => get_resolution(model),
            :horizon => get_horizon(get_settings(model)),
        )
    else
        export_params = nothing
    end

    write_model_dual_results!(store, model, index, update_timestamp, export_params)
    write_model_parameter_results!(store, model, index, update_timestamp, export_params)
    write_model_variable_results!(store, model, index, update_timestamp, export_params)
    write_model_aux_variable_results!(store, model, index, update_timestamp, export_params)
    write_model_expression_results!(store, model, index, update_timestamp, export_params)
    return
end

function write_model_dual_results!(
    store,
    model::T,
    index::Union{DecisionModelIndexType, EmulationModelIndexType},
    update_timestamp::Dates.DateTime,
    export_params::Union{Dict{Symbol, Any}, Nothing},
) where {T <: OperationModel}
    container = get_optimization_container(model)
    model_name = get_name(model)
    if export_params !== nothing
        exports_path = joinpath(export_params[:exports_path], "duals")
        mkpath(exports_path)
    end

    for (key, constraint) in get_duals(container)
        !should_write_resulting_value(key) && continue
        write_result!(store, model_name, key, index, update_timestamp, constraint)

        if export_params !== nothing &&
           should_export_dual(export_params[:exports], index, model_name, key)
            horizon = export_params[:horizon]
            resolution = export_params[:resolution]
            file_type = export_params[:file_type]
            df = axis_array_to_dataframe(constraint, key)
            time_col = range(index, length=horizon, step=resolution)
            DataFrames.insertcols!(df, 1, :DateTime => time_col)
            export_result(file_type, exports_path, key, index, df)
        end
    end
    return
end

function calaculate_parameter_values(
    param_array::DenseAxisArray,
    multiplier_array::DenseAxisArray,
)
    return jump_value.(param_array) .* multiplier_array
end

function calaculate_parameter_values(
    param_array::SparseAxisArray,
    multiplier_array::SparseAxisArray,
)
    p_array = jump_value.(to_matrix(param_array))
    m_array = to_matrix(multiplier_array)
    return p_array .* m_array
end

function write_model_parameter_results!(
    store,
    model::T,
    index::Union{DecisionModelIndexType, EmulationModelIndexType},
    update_timestamp::Dates.DateTime,
    export_params::Union{Dict{Symbol, Any}, Nothing},
) where {T <: OperationModel}
    container = get_optimization_container(model)
    model_name = get_name(model)
    if export_params !== nothing
        exports_path = joinpath(export_params[:exports_path], "parameters")
        mkpath(exports_path)
    end

    horizon = get_horizon(get_settings(model))

    parameters = get_parameters(container)
    for (key, container) in parameters
        !should_write_resulting_value(key) && continue
        param_array = get_parameter_array(container)
        multiplier_array = get_multiplier_array(container)
        # @assert_op length(axes(param_array)) == 2
        data = calaculate_parameter_values(param_array, multiplier_array)
        write_result!(store, model_name, key, index, update_timestamp, data)

        if export_params !== nothing &&
           should_export_parameter(export_params[:exports], index, model_name, key)
            resolution = export_params[:resolution]
            file_type = export_params[:file_type]
            df = axis_array_to_dataframe(data, key)
            time_col = range(index, length=horizon, step=resolution)
            DataFrames.insertcols!(df, 1, :DateTime => time_col)
            export_result(file_type, exports_path, key, index, df)
        end
    end
    return
end

function write_model_variable_results!(
    store,
    model::T,
    index::Union{DecisionModelIndexType, EmulationModelIndexType},
    update_timestamp::Dates.DateTime,
    export_params::Union{Dict{Symbol, Any}, Nothing},
) where {T <: OperationModel}
    container = get_optimization_container(model)
    model_name = get_name(model)
    if export_params !== nothing
        exports_path = joinpath(export_params[:exports_path], "variables")
        mkpath(exports_path)
    end

    if !isempty(container.primal_values_cache)
        variables = container.primal_values_cache.variables_cache
    else
        variables = get_variables(container)
    end

    for (key, variable) in variables
        !should_write_resulting_value(key) && continue
        write_result!(store, model_name, key, index, update_timestamp, variable)

        if export_params !== nothing &&
           should_export_variable(export_params[:exports], index, model_name, key)
            horizon = export_params[:horizon]
            resolution = export_params[:resolution]
            file_type = export_params[:file_type]
            df = axis_array_to_dataframe(variable, key)
            time_col = range(index, length=horizon, step=resolution)
            DataFrames.insertcols!(df, 1, :DateTime => time_col)
            export_result(file_type, exports_path, key, index, df)
        end
    end
    return
end

function write_model_aux_variable_results!(
    store,
    model::T,
    index::Union{DecisionModelIndexType, EmulationModelIndexType},
    update_timestamp::Dates.DateTime,
    export_params::Union{Dict{Symbol, Any}, Nothing},
) where {T <: OperationModel}
    container = get_optimization_container(model)
    model_name = get_name(model)
    if export_params !== nothing
        exports_path = joinpath(export_params[:exports_path], "aux_variables")
        mkpath(exports_path)
    end

    for (key, variable) in get_aux_variables(container)
        !should_write_resulting_value(key) && continue
        write_result!(store, model_name, key, index, update_timestamp, variable)

        if export_params !== nothing &&
           should_export_aux_variable(export_params[:exports], index, model_name, key)
            horizon = export_params[:horizon]
            resolution = export_params[:resolution]
            file_type = export_params[:file_type]
            df = axis_array_to_dataframe(variable, key)
            time_col = range(index, length=horizon, step=resolution)
            DataFrames.insertcols!(df, 1, :DateTime => time_col)
            export_result(file_type, exports_path, key, index, df)
        end
    end
    return
end

function write_model_expression_results!(
    store,
    model::T,
    index::Union{DecisionModelIndexType, EmulationModelIndexType},
    update_timestamp::Dates.DateTime,
    export_params::Union{Dict{Symbol, Any}, Nothing},
) where {T <: OperationModel}
    container = get_optimization_container(model)
    model_name = get_name(model)
    if export_params !== nothing
        exports_path = joinpath(export_params[:exports_path], "expressions")
        mkpath(exports_path)
    end

    if !isempty(container.primal_values_cache)
        expressions = container.primal_values_cache.expressions_cache
    else
        expressions = get_expressions(container)
    end

    for (key, expression) in expressions
        !should_write_resulting_value(key) && continue
        write_result!(store, model_name, key, index, update_timestamp, expression)

        if export_params !== nothing &&
           should_export_expression(export_params[:exports], index, model_name, key)
            horizon = export_params[:horizon]
            resolution = export_params[:resolution]
            file_type = export_params[:file_type]
            df = axis_array_to_dataframe(expression, key)
            time_col = range(index, length=horizon, step=resolution)
            DataFrames.insertcols!(df, 1, :DateTime => time_col)
            export_result(file_type, exports_path, key, index, df)
        end
    end
    return
end


text:
using DelimitedFiles, ProgressMeter, PyPlot, Gurobi

include("./util.jl")
include("./parsers.jl")
include("./Algorithms/TEAVAR.jl")

function timeScenarios(cutoffs, iterations, weibull=true, shape=.8, scale=.0001)
    x_vals = [18, 22, 26, 30, 34, 38, 42, 46, 50, 54, 58, 62, 66, 70, 74, 78, 82, 86, 90, 94, 100, 104, 108, 112]
    y_vals = [[] for i=1:length(cutoffs) + 1]
    p = Progress(length(x_vals)*iterations*(length(cutoffs)+1), .1, "Scenario Times...", 50)
    for k in 1:length(x_vals)
        vals = zeros(length(cutoffs) + 1)
        for i in 1:iterations
            if weibull
                probabilities = weibullProbs(x_vals[k], shape=.8, scale=.0001)
            else
                probabilities = map(n -> rand(n) .* 2 ./ 10, iterations)
            end
            cutoff = (sum(probabilities)/length(probabilities))^3
            next = @elapsed subScenariosRecursion(probabilities, cutoff)
            next!(p, showvalues = [(:x_vals,"$(k)/$(length(x_vals))"), (:iterations,"$(i)/$(iterations)"), (:cutoffs,cutoff), (:last, next)])
            vals[1] += next
            for j in 1:length(cutoffs)
                next = @elapsed subScenariosRecursion(probabilities, cutoffs[j])
                next!(p, showvalues = [(:x_vals,"$(k)/$(length(x_vals))"), (:iterations,"$(i)/$(iterations)"), (:cutoff,cutoffs[j]), (:last, next)])
                vals[j+1] += next
            end
        end
        vals ./ iterations
        for v in 1:length(vals)
            push!(y_vals[v], vals[v])
        end
    end

    z = zeros(length(x_vals), length(y_vals) + 1)
    z[:,1] = x_vals
    for i in 1:length(y_vals)
        z[:,i+1] = y_vals[i]
    end
    writedlm("./data/raw/time_scenarios/x_vals", x_vals)
    writedlm("./data/raw/time_scenarios/y_vals", y_vals)
    writedlm("./data/raw/time_scenarios/z", z)
end

function timeOptimizer(topologies,
                       cutoffs,
                       iterations;
                       weibull=true,
                       shape=.8,
                       scale=.0001,
                       paths="SMORE",
                       plot=true,
                       beta=.9,
                       demand_num=1,
                       dirname="./data/raw/time_optimizer")

    env = Gurobi.Env()

    y_vals = [[] for i=1:length(cutoffs) + 1]
    x_vals = []

    p = Progress(length(topologies)*iterations*(length(cutoffs)+1), .1, "Scenario Times...", 50)
    for t in 1:length(topologies)
        vals = zeros(length(cutoffs) + 1)
        links, capacity, link_probs, nodes = readTopology(topologies[t])
        demand, flows = readDemand("$(topologies[t])/demand", length(nodes), demand_num, matrix=true)
        T, Tf, k = parsePaths("$(topologies[t])/paths/$(paths)", links, flows)
        push!(x_vals, length(links))

        for i in 1:iterations
            if weibull
                probabilities = weibullProbs(length(links), shape=.8, scale=.0001)
            else
                probabilities = map(n -> rand(n) .* 2 ./ 10, iterations)
            end

            ## OPTIMAL
            cutoff = (sum(probabilities)/length(probabilities))^2
            scenarios, scenario_probs = subScenariosRecursion(probabilities, cutoff)
            next = @elapsed TEAVAR(env, links, capacity, flows, demand, beta, k, T, Tf, scenarios, scenario_probs)
            next!(p, showvalues = [(:topologies,"$(t)/$(length(topologies))"), (:iterations,"$(i)/$(iterations)"), (:cutoffs,cutoff), (:last, next)])
            vals[1] += next

            ## OTHERS
            for j in 1:length(cutoffs)
                scenarios, scenario_probs = subScenariosRecursion(probabilities, cutoffs[j])
                next = @elapsed TEAVAR(env, links, capacity, flows, demand, beta, k, T, Tf, scenarios, scenario_probs)
                next!(p, showvalues = [(:topologies,"$(t)/$(length(topologies))"), (:iterations,"$(i)/$(iterations)"), (:cutoff,cutoffs[j]), (:last, next)])
                vals[j+1] += next
            end
        end
        vals ./ iterations
        for v in 1:length(vals)
            push!(y_vals[v], vals[v])
        end
    end

    # LOG RESULTS
    dir = nextRun(dirname)
    z = zeros(length(x_vals), length(y_vals) + 1)
    z[:,1] = x_vals
    for i in 1:length(y_vals)
        z[:,i+1] = y_vals[i]
    end
    writedlm("$(dir)/x_vals", x_vals)
    writedlm("$(dir)/y_vals", y_vals)
    writedlm("$(dir)/z", z)

    # PLOT
    if plot
        PyPlot.clf()
        for i in 1:length(y_vals)
            PyPlot.plot(x_vals, y_vals[i])
        end
        PyPlot.xlabel("Number of Edges", fontweight="bold")
        PyPlot.ylabel("Time (s)", fontweight="bold")
        PyPlot.legend(pushfirst!(map(elt -> string(elt), cutoffs), "near optimal"), loc="lower right")
        PyPlot.show()
    end
end

# timeOptimizer(["./data/B4", "./data/IBM"], [.001, .0001, .00001, .000001, .0000001], 2)
# timeScenarios([.001, .0001, .00001, .000001, .0000001], 200)


text:
t = -3.1:0.01:3.1
u = 16*sin.(t).^3
v = 13*cos.(t) .- 5*cos.(2t) .- 2*cos.(3t) .- cos.(4t)
plot(u, v, c = :black, leg=false, fill=(0, 0.8,:black),
    xlab = L"u", ylab = L"v", xlim = (-20,20), ylim = (-20,15),
    size = (400, 400))
savefig("figheart.pdf")

max_stars_repo_name: lazarusA/helpfulJulia

max_issues_repo_name: lazarusA/helpfulJulia

max_forks_repo_name: lazarusA/helpfulJulia


text:
# Autogenerated wrapper script for NativeFileDialog_jll for aarch64-apple-darwin
export libnfd

JLLWrappers.@generate_wrapper_header("NativeFileDialog")
JLLWrappers.@declare_library_product(libnfd, "@rpath/libnfd.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libnfd,
        "lib/libnfd.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
`solve_system(f,u0,n)` #code

Solves the dynamical system

``u_{n+1} = f(u_n)`` #latex

for N steps. Returns the solution at step `n` with parameters `p`.

"""
function solve_system(f,u0,p,n)
  u = u0
  for i in 1:n-1
    u = f(u,p)
  end
  u
end

f(u,p) = u^2 - p*u
typeof(f)

function lorenz(u,p)
  α,σ,ρ,β = p
  du1 = u[1] + α*(σ*(u[2]-u[1]))
  du2 = u[2] + α*(u[1]*(ρ-u[3]) - u[2])
  du3 = u[3] + α*(u[1]*u[2] - β*u[3])
  [du1,du2,du3]
end
p = (0.02,10.0,28.0,8/3)

function solve_system_save(f,u0,p,n)
  u = Vector{typeof(u0)}(undef,n)
  u[1] = u0
  for i in 1:n-1
    u[i+1] = f(u[i],p)
  end
  u
end

a = solve_system_save(lorenz,[1.0,0.0,0.0],p,1000)

a[1]

reduce(hcat,a)

using BenchmarkTools
function solve_system_save_matrix(f,u0,p,n)
  u = Matrix{eltype(u0)}(undef,length(u0),n)
  u[:,1] = u0
  for i in 1:n-1
    u[:,i+1] = f(u[:,i],p)
  end
  u
end
@btime solve_system_save_matrix(lorenz,[1.0,0.0,0.0],p,1000)
@btime solve_system_save(lorenz,[1.0,0.0,0.0],p,1000)

function solve_system_save_matrix_view(f,u0,p,n)
  u = Matrix{eltype(u0)}(undef,length(u0),n)
  u[:,1] = u0
  for i in 1:n-1
    u[:,i+1] = f(@view(u[:,i]),p)
  end
  u
end
@btime  solve_system_save_matrix_view(lorenz,[1.0,0.0,0.0],p,1000)

function lorenz(du,u,p)
  α,σ,ρ,β = p
  du[1] = u[1] + α*(σ*(u[2]-u[1]))
  du[2] = u[2] + α*(u[1]*(ρ-u[3]) - u[2])
  du[3] = u[3] + α*(u[1]*u[2] - β*u[3])
end
p = (0.02,10.0,28.0,8/3)
function solve_system_save(f,u0,p,n)
  u = Vector{typeof(u0)}(undef,n)
  du = similar(u0) # create cache array, modify it, add to output
  u[1] = u0
  for i in 1:n-1
    f(du,u[i],p)
    u[i+1] = du
  end
  u
end
@btime solve_system_save(lorenz,[1.0,0.0,0.0],p,1000)

function solve_system_save(f,u0,p,n)
  u = Vector{typeof(u0)}(undef,n)
  du = similar(u0) # create cache array, modify it, add to output
  u[1] = u0
  for i in 1:n-1
    f(du,u[i],p)
    u[i+1] = du[:]
  end
  u
end
@btime solve_system_save(lorenz,[1.0,0.0,0.0],p,1000)

b = solve_system_save(lorenz,[1.0,0.0,0.0],p,1000)


text:
function compare(E, idx, ∇x, ∇y, ∇z)
    N = size(E, 5)
    println("Checking X")
    for i = 1:N
        f = true
        for j = i:N
            for k = 1:N
                for l = k:N
                    if GaussianBasis.index2(i-1,j-1) > GaussianBasis.index2(k-1,l-1)
                        continue
                    end
                    A = 0.0
                    try
                        Aidx = findfirst(x->x==(i-1, j-1, k-1, l-1), idx)
                        A = ∇x[Aidx]
                    catch ArgumentError
                        A = 0.0
                    end
                    B = E[1,i,j,k,l]
                    test = A ≈ B
                    if !test
                        println("$i $j $k $l")
                        f = false
                    end
                end
                !f ? break : nothing
            end
            !f ? break : nothing
        end
        println(f)
        !f ? break : nothing
    end

    println("Checking Y")
    for i = 1:N
        f = true
        for j = i:N
            for k = 1:N
                for l = k:N
                    if GaussianBasis.index2(i-1,j-1) > GaussianBasis.index2(k-1,l-1)
                        continue
                    end
                    A = 0.0
                    try
                        Aidx = findfirst(x->x==(i-1, j-1, k-1, l-1), idx)
                        A = ∇y[Aidx]
                    catch ArgumentError
                        A = 0.0
                    end
                    B = E[2,i,j,k,l]
                    test = A ≈ B
                    if !test
                        println("$i $j $k $l")
                        f = false
                    end
                end
                !f ? break : nothing
            end
            !f ? break : nothing
        end
        println(f)
        !f ? break : nothing
    end

    println("Checking Z")
    for i = 1:N
        f = true
        for j = i:N
            for k = 1:N
                for l = k:N
                    if GaussianBasis.index2(i-1,j-1) > GaussianBasis.index2(k-1,l-1)
                        continue
                    end
                    A = 0.0
                    try
                        Aidx = findfirst(x->x==(i-1, j-1, k-1, l-1), idx)
                        A = ∇z[Aidx]
                    catch ArgumentError
                        A = 0.0
                    end
                    B = E[3,i,j,k,l]
                    test = A ≈ B
                    if !test
                        println("$i $j $k $l")
                        f = false
                    end
                end
                !f ? break : nothing
            end
            !f ? break : nothing
        end
        println(f)
        !f ? break : nothing
    end
end