text:
function red_reg(::Type{TT}, sg::Spinglass{LT,T}, Ly::Int, k::Int, hk::Int; usecuda=false) where {LT,T,TT}
    reg = _init_reg(TT, Ly*4, Val(usecuda))
    for i=1:Ly*4
        hk += 1
        sg.hs[hk] != 0 && (reg |> put(4*Ly, i=>Gh(vertextensor(TT, sg, hk+i))))
    end
    for i=1:Ly-1
        for j = 1:4
            k += 1
            reg |> put(4*Ly, (4*(i-1)+j,4*i+j)=>Gvb(bondtensor(TT, sg, k)))
        end
    end
    for i=1:Ly
        reg |> put(4*Ly, (4i-3,4i-2,4i-1,4i)=>G16(TT, sg.Js[k+1:k+16]))
        k += 16
    end
    return reg
end

function solve(::Type{TT}, sg::Spinglass{LT,T}; usecuda::Bool) where {LT<:ChimeraLattice,T, TT}
    _, Lx, Ly = size(sg.lattice)
    Js = sg.Js
    hs = sg.hs
    println("Running Chimera, Lx = $Lx, Ly = $Ly, eltype = $T, usecuda = $usecuda.")
    nj_red = (Ly-1)*4 + 16 * Ly
    println("Layer 1/$Lx")
    k = 0
    reg = red_reg(TT, sg, Ly, k, 0; usecuda=usecuda)
    k += nj_red
    for j=1:Ly*4
        hs[j]!= 0 && (reg |> put(4*Ly, j=>Gh(vertextensor(TT, sg, j))))
    end
    for i=2:Lx
        hk = (i-1)*Ly*8
        println("Layer $i/$Lx")
        # BLACK
        for j=1:Ly*4
            k += 1
            reg |> put(Ly*4, j=>Ghb(bondtensor(TT, sg, k)))
        end
        for j=1:Ly*4 # `hs` interated in red->black order
            hs[j+hk+4Ly] != 0 && (reg |> put(4*Ly, j=>Gh(vertextensor(TT, sg, j+hk+4Ly))))
        end

        # Contract with RED
        rr = red_reg(TT, sg, Ly, k, hk; usecuda=usecuda)
        k += nj_red
        reg.state .*= rr.state
    end
    if length(Js) !== k
        @warn "length of parameters is $(length(Js)), used $k"
    end
    sum(statevec(reg))
end

# NOTE: repeat with the definition in Viznet.
function sgbonds(lt::ChimeraLattice)
    bonds = Tuple{Int,Int}[]
    red_bond!(lt, 1, bonds)
    redblack_bond!(lt, 1, bonds)
    for i=2:lt.Nx
        # BLACK
        black_bond!(lt, i-1, bonds)
        # Contract with RED
        red_bond!(lt, i, bonds)
        redblack_bond!(lt, i, bonds)
    end
    return bonds
end

function red_bond!(lt, i::Int, bonds)
    CI = LinearIndices(lt |> size)
    for j=1:lt.Ny-1
        for k = 1:4
            push!(bonds, (CI[k, i, j], CI[k, i, j+1]))
        end
    end
    return bonds
end

function redblack_bond!(lt, i::Int, bonds)
    CI = LinearIndices(lt |> size)
    for j=1:lt.Ny
        for k2=5:8
            for k1=1:4
                push!(bonds, (CI[k2, i, j], CI[k1, i, j]))
            end
        end
    end
    return bonds
end

function black_bond!(lt, i::Int, bonds)
    CI = LinearIndices(lt |> size)
    for j=1:lt.Ny
        for k = 5:8
            push!(bonds, (CI[k,i,j], CI[k,i+1,j]))
        end
    end
    return bonds
end

# to index `h`.
function sgvertices(lt::ChimeraLattice)
    v = Int[]
    LI = LinearIndices(lt)
    for i=1:lt.Nx
        for ki = 1:2
            for j=1:lt.Ny
                for kj = 1:4
                    k = (ki-1)*4 + kj
                    push!(v, LI[k,i,j])
                end
            end
        end
    end
    return v
end

regsize(lt::ChimeraLattice) = lt.Ny*4

function _init_reg(::Type{T}, lt::ChimeraLattice, usecuda) where T<:TropicalTypes
    _init_reg(T, lt.Ny*4, usecuda)
end


text:
using HTTP
using JSON

include("constants.jl")

include("src/telegram_api.jl")
include("src/update_handling.jl")

function run_bot()
    offset = 0

    while true
        updates = try
            get_updates(TOKEN, offset)
        catch
            @warn "Something went wrong while getting updates from Telegram."
            continue
        end

        if length(updates) == 0
            continue
        else
            offset = maximum([u["update_id"] for u in updates]) + 1
        end

        for message in updates
            @async handle_update(TOKEN, message)
        end

        sleep(1)
    end
end

run_bot()

ext: jl

max_stars_repo_name: iwsylit/spotify-utils-bot

max_issues_repo_name: iwsylit/spotify-utils-bot

max_forks_repo_name: iwsylit/spotify-utils-bot


text:
#####
##### General halo filling functions
#####

fill_halo_regions!(::Nothing, args...) = nothing

"""
    fill_halo_regions!(fields, bcs, arch, grid)

Fill halo regions for each field in the tuple `fields` according
to the single instance of `FieldBoundaryConditions` in `bcs`, possibly recursing into
`fields` if it is a nested tuple-of-tuples.
"""
function fill_halo_regions!(fields::Union{Tuple, NamedTuple}, bcs::FieldBoundaryConditions, arch, grid, args...)
    for field in fields
        fill_halo_regions!(field, bcs, arch, grid, args...)
    end
    return nothing
end

"""
    fill_halo_regions!(fields, bcs, arch, grid)

Fill halo regions for all fields in the `NamedTuple` `fields` according
to the corresponding `NamedTuple` of `bcs`.
"""
function fill_halo_regions!(fields::NamedTuple{S}, bcs::NamedTuple{S}, arch, grid, args...) where S
    for (field, fieldbcs) in zip(fields, bcs)
        fill_halo_regions!(field, fieldbcs, arch, grid, args...)
    end
    return nothing
end

"Fill halo regions in x, y, and z for a given field."
function fill_halo_regions!(c::AbstractArray, fieldbcs, arch, grid, args...)
      fill_west_halo!(c, fieldbcs.x.left,  arch, grid, args...)
      fill_east_halo!(c, fieldbcs.x.right, arch, grid, args...)

     fill_south_halo!(c, fieldbcs.y.left,  arch, grid, args...)
     fill_north_halo!(c, fieldbcs.y.right, arch, grid, args...)

     fill_bottom_halo!(c, fieldbcs.z.bottom, arch, grid, args...)
        fill_top_halo!(c, fieldbcs.z.top,    arch, grid, args...)
    return nothing
end


text:
struct DiGraph
    n :: Int32 # |V|
    m :: Int32 # |E|
    V :: Array{String, 1} # V[i] = Real Label of node i
    E :: Array{Tuple{Int32, Int32, Float64}, 1} # <u, v, w> in Edge Set
end

function readDiGraph(fileName) # Read a directed graph from a file
    n = 0
    V = Array{String, 1}()
    label = Dict{String, Int32}()
    W = Dict{Tuple{Int32, Int32}, Float64}()

    getid(x) = begin
        if !haskey(label, x)
            n += 1
            label[x] = n
            push!(V, x)
            return n
        end
        return label[x]
    end

    open(fileName) do f1
        for line in eachline(f1)
            buf = split(line)
            u = getid(buf[1])
            v = getid(buf[2])
            if u == v
                continue
            end
            if haskey(W, (u, v))
                W[(u, v)] += 1.0
            else
                W[(u, v)] = 1.0
            end
        end
    end

    m = length(W)
    E = Array{Tuple{Int32, Int32, Float64}, 1}()
    for (u, v) in keys(W)
        push!(E, (u, v, W[(u, v)]))
    end

    return DiGraph(n, m, V, E)
end

function getAdjacentList(G :: DiGraph)
    g = Vector{Vector{Int32}}(undef, G.n)
    foreach(i -> g[i] = Vector{Int32}(), 1 : G.n)
    for (u, v, w) in G.E
        push!(g[u], v)
    end
    return g
end

function SubDiGraph(G :: DiGraph, C, xx)
    n = 0
    h = zeros(Int32, G.n)
    V = Vector{String}()
    for i = 1 : G.n
        if C[i] == xx
            n += 1
            h[i] = n
            push!(V, G.V[i])
        end
    end

    E = Vector{Tuple{Int32, Int32, Float64}}()
    for (u, v, w) in G.E
        if h[u] != 0 && h[v] != 0
            push!(E, (h[u], h[v], w))
        end
    end

    return DiGraph(n, size(E, 1), V, E)
end

function LSCC(G :: DiGraph)
    g = getAdjacentList(G)
    dfn = zeros(Int32, G.n)
    low = zeros(Int32, G.n)
    index = 0
    S = Vector{Int32}()
    C = zeros(Int32, G.n)
    inStack = zeros(Bool, G.n)

    dfs(u) = begin
        index += 1
        dfn[u] = index
        low[u] = index
        push!(S, u)
        inStack[u] = true
        for v in g[u]
            if dfn[v] == 0
                dfs(v)
                low[u] = min(low[u], low[v])
            elseif inStack[v]
                low[u] = min(low[u], dfn[v])
            end
        end
        if dfn[u] == low[u]
            while true
                v = pop!(S)
                inStack[v] = false
                C[v] = u
                if u == v
                    break
                end
            end
        end
    end

    for i = 1 : G.n
        if dfn[i] == 0
            dfs(i)
        end
    end

    sz = zeros(Int32, G.n)
    foreach(x -> sz[x] += 1, C)
    xx = argmax(sz)

    return SubDiGraph(G, C, xx)
end


text:
function standardize(data::Matrix)
  m,n = size(data)
  standardized = zeros(m,n)
  mean_vector = zeros(n)
  std_vector = zeros(n)
  for i = 1:n
    mean_vector[i] = mean(data[:,i])
    std_vector[i] = std(data[:,i])
    if std_vector[i] != 0
      standardized[:,i] = (data[:,i]-mean_vector[i])./std_vector[i]
    elseif mean_vector[i] < 1 && mean_vector[i] >=0
      standardized[:,i] = data[:,i]
    else
      standardized[:,i] = 1
    end
  end
  return standardized::Matrix{Float64}, mean_vector::Vector{Float64}, std_vector::Vector{Float64}
end

function standardize(data::Matrix, mean_vector::Vector{Float64}, std_vector::Vector{Float64})
  m,n = size(data)
  standardized = zeros(m,n)
  for i = 1:n
    if std_vector[i] != 0
      standardized[:,i] = (data[:,i]-mean_vector[i]) ./ std_vector[i]
    elseif mean_vector[i] < 1 && mean_vector[i] >=0
      standardized[:,i] = data[:,i]
    else
      standardized[:,i] = 1
    end
  end
  return standardized
end

function normalize(X::Matrix)
  # normalize to 0 1
  r_nrm = 1.0   # range
  shift = 0.0
  X_max = maximum(X, 1)
  X_min = minimum(X, 1)
  X_nrm = (r_nrm * broadcast(-, X, X_min) ./ broadcast(-, X_max, X_min)) + shift

  return X_nrm
end

function k_fold(n::Int, k::Int)
  idx = randperm(n)

  # allocate folds
  folds = Vector[]
  n_f = round(Int, floor(n/k))
  add_f = n % k
  j = 1
  for i=1:k
    if i <= add_f
      push!(folds, idx[j:j+n_f])
      j += n_f+1
    else
      push!(folds, idx[j:j+n_f-1])
      j += n_f
    end
  end

  return folds::Vector{Vector}
end

function jaakkola_heuristic(X::Matrix, y::Vector)
  n = size(X,1)
  Xt = X'   # transpose is more efficient
  min_dist = zeros(n)
  for i=1:n
    ds = Inf
    for j=1:n
      if y[i] != y[j]
        d = norm(view(Xt,:,i) - view(Xt,:,j))
        if d < ds
          ds = d
        end
      end
    end
    min_dist[i] = ds
  end

  sigma = median(min_dist)
  gamma = 0.5 / (sigma * sigma)

  return gamma
end

function count_constraints(model::MultiAdversarialModel, y_train::Vector)
  w = model.w
  alpha = model.alpha
  cs = model.constraints

  n_train = length(y_train)

  # via dec or full
  if cs[1][1] == 1 && cs[1][2] == [y_train[1]]
    optim = :dec
    cs_added = cs[n_train+1:end]
  else
    optim = :full
    cs_added = cs
  end

  n_cs_added = length(cs_added)
  const_num = zeros(n_train)
  for i=1:n_train
    const_num[i] = length(find(x::Tuple{Integer, Vector} -> x[1] == i, cs_added))
  end
  max_cs_added = maximum(const_num)

  zero_thereshold = 1e-3
  active_const_num = zeros(n_train)
  for i=1:n_train
    idx = find(x::Tuple{Integer, Vector} -> x[1] == i, cs_added)
    active_const_num[i] = sum(alpha[idx] .> zero_thereshold)
  end

  n_cs_active = sum(active_const_num)
  max_cs_active = maximum(active_const_num)
  n_sv = sum(active_const_num .> 0)

  return n_cs_added, max_cs_added, n_cs_active, max_cs_active, n_sv
end

function count_constraints(model::KernelMultiAdversarialModel, y_train::Vector)
  alpha = model.alpha
  cs = model.constraints

  n_train = length(y_train)

  # via dec or full
  if cs[1][1] == 1 && cs[1][2] == [y_train[1]]
    optim = :dec
    cs_added = cs[n_train+1:end]
  else
    optim = :full
    cs_added = cs
  end

  n_cs_added = length(cs_added)
  const_num = zeros(n_train)
  for i=1:n_train
    const_num[i] = length(find(x::Tuple{Integer, Vector} -> x[1] == i, cs_added))
  end
  max_cs_added = maximum(const_num)

  zero_thereshold = 1e-3
  active_const_num = zeros(n_train)
  for i=1:n_train
    idx = find(x::Tuple{Integer, Vector} -> x[1] == i, cs_added)
    active_const_num[i] = sum(alpha[idx] .> zero_thereshold)
  end

  n_cs_active = sum(active_const_num)
  max_cs_active = maximum(active_const_num)
  n_sv = sum(active_const_num .> 0)

  return n_cs_added, max_cs_added, n_cs_active, max_cs_active, n_sv
end

max_stars_repo_name: Koprvhdix/adversarial-multiclass

max_issues_repo_name: Koprvhdix/adversarial-multiclass

max_forks_repo_name: Koprvhdix/adversarial-multiclass


text:
# export
#     ParamsEP,
#     epNP,
#     epdNP,
#     epMP,
#     epSNEP

immutable ParamsEP
    prior::EFamily
    logfactors::Vector{Function} # NOTE these functions need to be able to handle matrices of size pxN
    nEP::Int
    nIS::Int    # NOTE would have to generalise this --> could be SAMPLER obj
    damp::Float # NOTE would have to generalise this (stepping)
    minvar::Float
    # alpha::Float
    # ------------
    dim::Int
    nfactors::Int
    function ParamsEP(prior,factors;
                        nEP=100,nIS=100,damp=0.1,minvar=1e-5)
        @assert nIS>1 "Need at least two samples per estimator"
        #@assert 0.0<damp<=1.0 "Damping parameter should be in (0,1]"
        new(prior, factors, nEP, nIS, damp, minvar,
            length(mean(prior)), length(factors))
    end
end

"""
    epNP

Serial EP updates with damping in the Natural Parameter space.
This is the "SMS" style update.
"""
function epNP(p::ParamsEP)
    globapprox_mem = []

    # For each factor, we start with a very flat spherical gaussian
    locapprox  = [ones(GaussianNatParam, p.dim)/100 for i in 1:p.nfactors]
    globapprox = p.prior + sum(locapprox)

    push!(globapprox_mem, globapprox)

    for iter in 1:p.nEP
        for i in randperm(p.nfactors)
            cavity      = globapprox - locapprox[i]
            lltilted(x) = p.logfactors[i](x) + uloglik(cavity, x)
            (iss,w)     = impsampling(lltilted, globapprox, p.nIS)
            muhat       = sum( w[j]*suffstats(GaussianNatParam, iss[:,j])
                                for j in 1:p.nIS )
            # Cast it as mean parameter
            muhat = GaussianMeanParam(muhat[1], muhat[2])
            # Transform to natural parameter, project in case out of NP space
            np_muhat     = natparam(project(muhat,minvar=p.minvar))
            globapprox   = (1.0-p.damp)*globapprox + p.damp*np_muhat
            locapprox[i] = globapprox-cavity
        end
        push!(globapprox_mem, globapprox)
    end
    (globapprox, globapprox_mem)
end

"""
    epNP

Serial EP updates with damping in the Natural Parameter space.
"""
function epdNP(p::ParamsEP)
    locapprox  = [ones(DiagGaussianNatParam, p.dim)/100 for i in 1:p.nfactors]
    globapprox = p.prior + sum(locapprox)
    updatecounter = 0
    for iter in 1:p.nEP
        for i in randperm(p.nfactors)
            cavity      = globapprox - locapprox[i]
            lltilted(x) = p.logfactors[i](x) + uloglik(cavity, x)
            (iss,w)     = impsampling(lltilted, globapprox, p.nIS)
            muhat       = sum( w[j]*suffstats(DiagGaussianNatParam, iss[:,j])
                                for j in 1:p.nIS )
            # Cast it as mean parameter
            muhat = DiagGaussianMeanParam(muhat[1], muhat[2])
            # Transform to natural parameter, project in case out of NP space
            np_muhat = natparam(muhat)
            # Perform the tentative damping step
            t_globapprox = (1.0-p.damp)*globapprox+p.damp*np_muhat
            if maximum(-t_globapprox.theta2) < 1.0/p.minvar
                globapprox   = t_globapprox
                locapprox[i] = globapprox-cavity
                updatecounter += 1
            end
        end
    end
    println("Ratio of accepted updates: $updatecounter/$(p.nEP*p.nfactors)")
    globapprox
end

"""
    epMP

Serial EP updates with damping in the Mean Parameter space.
"""
function epMP(p::ParamsEP)
    globapprox_mem = []

    locapprox  = [ones(GaussianNatParam, p.dim)/100 for i in 1:p.nfactors]
    globapprox = p.prior + sum(locapprox)

    push!(globapprox_mem, globapprox)

    for iter in 1:p.nEP
        for i in randperm(p.nfactors)
            cavity      = globapprox - locapprox[i]
            lltilted(x) = p.logfactors[i](x) + uloglik(cavity, x)
            (iss,w)     = impsampling(lltilted, globapprox, p.nIS)
            muhat       = sum( w[j]*suffstats(GaussianNatParam, iss[:,j])
                                for j in 1: p.nIS )
            globapprox = natparam( (1.0-p.damp)*meanparam(globapprox)+
                                    p.damp*GaussianMeanParam(muhat[1], muhat[2])
                                  )
            locapprox[i] = globapprox-cavity
        end
        push!(globapprox_mem, globapprox)
    end
    (globapprox, globapprox_mem)
end

"""
    epSNEP

Serial SNEP updates attempt
"""
function epSNEP(p::ParamsEP)
    globapprox_mem = []

    locapprox  = [ones(GaussianNatParam, p.dim)/100 for i in 1:p.nfactors]
    globapprox = p.prior + sum(locapprox)

    push!(globapprox_mem, globapprox)

    for iter in 1:p.nEP
        for i in randperm(p.nfactors)
            cavity      = globapprox - locapprox[i]
            lltilted(x) = p.logfactors[i](x) + uloglik(cavity, x)
            (iss,w)     = impsampling(lltilted, globapprox, p.nIS)
            muhat       = sum( w[j]*suffstats(GaussianNatParam, iss[:,j])
                                for j in 1: p.nIS )
            locapprox[i] = natparam( meanparam(locapprox[i]) + p.damp*(
                                GaussianMeanParam(muhat[1], muhat[2]) -
                                meanparam(globapprox)
                                )
                              )
            globapprox = cavity + locapprox[i]
        end
        push!(globapprox_mem, globapprox)
    end
    (globapprox, globapprox_mem)
end


text:
"""
Macro that inputs an expression corresponding to a reaction network and output a Reaction Network Structure that can be used as input to generation of SDE and ODE and Jump problems.
Most arrows accepted (both right, left and bi drectional arrows).
Note that while --> is a correct arrow, neither <-- nor <--> works.
Using non-filled arrows (⇐, ⟽, ⇒, ⟾, ⇔, ⟺) will disable mass kinetics and lets you cutomize reaction rates yourself.
Use 0 or ∅ for degradation/creation to/from nothing.
Example systems:
    ### Basic Usage ###
    rn = @reaction_network rType begin #Creates a reaction network of type rType.
        2.0, X + Y --> XY                  #This will have reaction rate corresponding to 2.0*[X][Y]
        2.0, XY ← X + Y                    #Identical to 2.0, X + Y --> XY
    end

    ### Manipulating Reaction Rates ###
    rn = @reaction_network rType begin
        2.0, X + Y ⟾ XY                   #Ignores mass kinetics. This will have reaction rate corresponding to 2.0.
        2.0X, X + Y --> XY                 #Reaction rate needs not be constant. This will have reaction rate corresponding to 2.0*[X]*[X]*[Y].
        XY+log(X)^2, X + Y --> XY          #Reaction rate accepts quite complicated expressions (user defined functions must first be registered using the @reaction_func macro).
        hill(XY,2,2,2), X + Y --> XY       #Reaction inis activated by XY according to a hill function. hill(x,v,K,N).
        mm(XY,2,2), X + Y --> XY           #Reaction inis activated by XY according to a michaelis menten function. mm(x,v,K).
    end

    ### Multipple Reactions on a Single Line ###
    rn = @reaction_network rType begin
        (2.0,1.0), X + Y ↔ XY              #Identical to reactions (2.0, X + Y --> XY) and (1.0, XY --> X + Y).
        2.0, (X,Y) --> 0                   #This corresponds to both X and Y degrading at rate 2.0.
        (2.0, 1.0), (X,Y) --> 0            #This corresponds to X and Y degrading at rates 2.0 and 1.0, respectively.
        2.0, (X1,Y1) --> (X2,Y2)           #X1 and Y1 becomes X2 and Y2, respectively, at rate 2.0.
    end

    ### Adding Parameters ###
    kB = 2.0; kD = 1.0
    p = [kB, kD]
    p = []
    rn = @reaction_network type begin
        (kB, kD), X + Y ↔ XY            #Lets you define parameters outside on network. Parameters can be changed without recalling the network.
    end kB, kD

    ### Defining New Functions ###
    @reaction_func my_hill_repression(x, v, k, n) = v*k^n/(k^n+x^n)     #Creates and adds a new function that the @reaction_network macro can see.
    r = @reaction_network MyReactionType begin
        my_hill_repression(x, v_x, k_x, n_x), 0 --> x                       #After it has been added in @reaction_func the function can be used when defining new reaction networks.
    end v_x k_x n_x

    ### Simulating Reaction Networks ###
    probODE = ODEProblem(rn, args...; kwargs...)        #Using multiple dispatch the reaction network can be used as input to create ODE, SDE and Jump problems.
    probSDE = SDEProblem(rn, args...; kwargs...)
    probJump = JumpProblem(prob,aggregator::Direct,rn)
"""

"""
    @reaction_network

Generates a subtype of an `AbstractReactionNetwork` that encodes a chemical
reaction network, and complete ODE, SDE and jump representations of the system.
See the [Chemical Reaction Model
docs](http://docs.juliadiffeq.org/latest/models/biological.html) for details on
parameters to the macro.
"""
macro reaction_network(name, ex::Expr, p...)
    coordinate(name, MacroTools.striplines(ex), p, :no___noise___scaling)
end

#Macro to create a reaction network model. Multiple dispatch is used to allow for SDE noise scalling.
macro reaction_network(name, scale_noise, ex::Expr, p...)
    in(scale_noise, p) || (p = (p..., scale_noise))
    coordinate(name, MacroTools.striplines(ex), p, scale_noise)
end

#If no type name is given, creates a network with a default name.
macro reaction_network(ex::Expr, p...)
    coordinate(:reaction_network, MacroTools.striplines(ex), p, :no___noise___scaling)
end

################# query-based macros:

"""
    @min_reaction_network

Generates a subtype of an `AbstractReactionNetwork` that only encodes a chemical
reaction network. Use [`addodes!`](@ref), [`addsdes!`](@ref) or
[`addjumps!`](@ref) to complete the network for specific problem types. It accepts
the same arguments as [`@reaction_network`](@ref).
"""
macro min_reaction_network(name, ex::Expr, p...)
    min_coordinate(name, MacroTools.striplines(ex), p, :no___noise___scaling)
end

#Macro to create a reaction network model. Multiple dispatch is used to allow for SDE noise scalling.
macro min_reaction_network(name, scale_noise, ex::Expr, p...)
    in(scale_noise, p) || (p = (p..., scale_noise))
    min_coordinate(name, MacroTools.striplines(ex), p, scale_noise)
end

#If no type name is given, creates a network with a default name.
macro min_reaction_network(ex::Expr, p...)
    min_coordinate(:min_reaction_network, MacroTools.striplines(ex), p, :no___noise___scaling)
end

"""
    @empty_reaction_network networktype

Generates a subtype of an `AbstractReactionNetwork` that encodes an empty
chemical reaction network. `networktype` is an optional parameter that specifies
the type of the generated network. Use [`addspecies!`](@ref), [`addparam!`](@ref)
and [`addreaction!`](@ref) to extend the network.  Use [`addodes!`](@ref),
[`addsdes!`](@ref) or [`addjumps!`](@ref) to complete the network for specific
problem types.
"""
macro empty_reaction_network(name::Symbol=:min_reaction_network)
    min_coordinate(name, MacroTools.striplines(:(begin end)), (), :no___noise___scaling)
end

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

#Declare various arrow types symbols used for the empty set (also 0).
empty_set = Set{Symbol}([:∅])
fwd_arrows = Set{Symbol}([:>, :→, :↣, :↦, :⇾, :⟶, :⟼, :⥟, :⥟, :⇀, :⇁, :⇒, :⟾])
bwd_arrows = Set{Symbol}([:<, :←, :↢, :↤, :⇽, :⟵, :⟻, :⥚, :⥞, :↼, :↽, :⇐, :⟽])
double_arrows = Set{Symbol}([:↔, :⟷, :⇄, :⇆, :⇔, :⟺])
no_mass_arrows = Set{Symbol}([:⇐, :⟽, :⇒, :⟾, :⇔, :⟺])      #Using this arrows will disable the program from multiplying reaction rates with the substrate concentrations. Gives user full control of reaction rates.

funcdict = Dict{Symbol, Function}()                             #Stores user defined functions.

#Coordination function, actually does all the work of the macro.
function coordinate(name, ex::Expr, p, scale_noise)

    # minimal reaction network components
    (reactions, reactants, parameters, syms, params) = get_minnetwork(ex, p)

    # expressions for ODEs
    (f_expr, f, f_rhs, symjac, jac, paramjac, f_symfuncs, jac_prototype) = genode_exprs(reactions, reactants, parameters, syms)
    odefun = :(ODEFunction(f; jac=$jac, jac_prototype=$jac_prototype, paramjac=$paramjac, syms=$syms))

    # expressions for SDEs
    (g_expr, g, g_funcs, p_matrix) = gensde_exprs(reactions, reactants, parameters, scale_noise)
    sdefun = :(SDEFunction(f, g; jac=$jac, jac_prototype=$jac_prototype, paramjac=$paramjac, syms=$syms))

    # expressions for jumps
    (jump_rate_expr, jump_affect_expr, jumps) = get_jumps(reactions, reactants, parameters)

    # expressions for regular jumps
    regular_jumps = get_regularjumps(reactions, reactants, parameters)

    # Build the type
    exprs = Vector{Expr}(undef,0)
    typeex,constructorex = maketype(DiffEqBase.AbstractReactionNetwork, name, f, f_rhs, f_symfuncs, g, g_funcs, jumps, regular_jumps, Meta.quot(jump_rate_expr), Meta.quot(jump_affect_expr), p_matrix, syms, scale_noise; params=params, reactions=reactions, jac=jac, paramjac=paramjac, symjac=symjac, syms_to_ints=reactants, params_to_ints=parameters, odefun=odefun, sdefun=sdefun)
    push!(exprs,typeex)
    push!(exprs,constructorex)

    # add type functions
    append!(exprs, gentypefun_exprs(name))

    # return as one expression block
    expr_arr_to_block(exprs)
end

# min_reaction_network coordination function, actually does all the work of the macro.
function min_coordinate(name, ex::Expr, p, scale_noise)

    # minimal reaction network components
    (reactions, reactants, parameters, syms, params) = get_minnetwork(ex, p)

    # Build the type
    exprs = Vector{Expr}(undef,0)
    typeex,constructorex = maketype(DiffEqBase.AbstractReactionNetwork, name, nothing, nothing,
                                    nothing, nothing, nothing, nothing, nothing, nothing,
                                    nothing, nothing, syms, scale_noise; params=params,
                                    reactions=reactions, symjac=nothing, syms_to_ints=reactants,
                                    params_to_ints=parameters)
    push!(exprs,typeex)
    push!(exprs,constructorex)

    # add type functions
    append!(exprs, gentypefun_exprs(name, gen_inplace=false, gen_outofplace=false))

    # return as one expression block
    expr_arr_to_block(exprs)
end

# ODE expressions
function genode_exprs(reactions, reactants, parameters, syms; build_jac=true,
                                                              build_paramjac=true,
                                                              sparse_jac=false,
                                                              build_symfuncs=true,
                                                              zeroout_jac=false)
    f_expr     = get_f(reactions, reactants)
    f          = make_func(f_expr, reactants, parameters)
    f_rhs      = ExprValues[element.args[2] for element in f_expr]
    f_symfuncs = build_symfuncs ? hcat([SymEngine.Basic(f) for f in f_rhs]) : nothing
    symjac, jac, jac_protoype, paramjac = get_jacs(f_rhs, syms, reactions, reactants, parameters; build_jac=build_jac, sparse_jac=sparse_jac, build_paramjac=build_paramjac, zeroout_jac=zeroout_jac)

    (f_expr,f,f_rhs,symjac,jac,paramjac,f_symfuncs, jac_protoype)
end

# SDE expressions
function gensde_exprs(reactions, reactants, parameters, scale_noise)
    g_expr   = get_g(reactions, reactants, scale_noise)
    g        = make_func(g_expr, reactants, parameters)
    g_funcs  = ExprValues[element.args[2] for element in g_expr]
    p_matrix = zeros(length(reactants), length(reactions))

    (g_expr,g,g_funcs,p_matrix)
end

# generate the minimal network components
function get_minnetwork(ex::Expr, p)
    reactions = get_reactions(ex)           ::Vector{ReactionStruct}
    reactants = get_reactants(reactions)    ::OrderedDict{Symbol,Int}
    parameters = get_parameters(p)          ::OrderedDict{Symbol,Int}

    syms = collect(keys(reactants))
    params = collect(keys(parameters))
    (in(:t,union(syms,params))) && error("t is reserved for the time variable and may neither be used as a reactant nor a parameter")

    update_reaction_info(reactions,syms)

    (reactions,reactants,parameters,syms,params)
end

#Generates a vector containing a number of reaction structures, each containing the infromation about one reaction.
function get_reactions(ex::Expr, reactions = Vector{ReactionStruct}(undef,0))
    for line in ex.args
        (line.head != :tuple) && (continue)
        (rate,r_line) = line.args
        if r_line.head  == :-->
            r_line = Expr(:call,:→,r_line.args[1],r_line.args[2])
        end

        arrow = r_line.args[1]  ::Symbol
        if in(arrow,double_arrows)
            (typeof(rate) == Expr && rate.head == :tuple) || error("Error: Must provide a tuple of reaction rates when declaring a bi-directional reaction.")
            push_reactions(reactions, r_line.args[2], r_line.args[3], rate.args[1], !in(arrow,no_mass_arrows))
            push_reactions(reactions, r_line.args[3], r_line.args[2], rate.args[2], !in(arrow,no_mass_arrows))
        elseif in(arrow,fwd_arrows)
            push_reactions(reactions, r_line.args[2], r_line.args[3], rate, !in(arrow,no_mass_arrows))
        elseif in(arrow,bwd_arrows)
            push_reactions(reactions, r_line.args[3], r_line.args[2], rate, !in(arrow,no_mass_arrows))
        else
            throw("malformed reaction")
        end
    end
    return reactions
end

#Structure containing information about one reactant in one reaction.
struct ReactantStruct
    reactant::Symbol
    stoichiometry::Int
end

#Structure containing information about one Reaction. Contain all its substrates and products as well as its rate. Contains an specialized constructor.
struct ReactionStruct
    substrates::Vector{ReactantStruct}
    products::Vector{ReactantStruct}
    netstoich::Vector{ReactantStruct}
    rate_org::ExprValues
    rate_DE::ExprValues
    rate_SSA::ExprValues
    dependants::Vector{Symbol}
    is_pure_mass_action::Bool

    function ReactionStruct(s::Vector{ReactantStruct}, p::Vector{ReactantStruct},
                            ns::Vector{ReactantStruct},
                            ro::ExprValues, rde::ExprValues, rssa::ExprValues,
                            dep::Vector{Symbol}, isma::Bool)
        new(s,p,ns,ro,rde,rssa,dep,isma)
    end

    function ReactionStruct(sub_line::ExprValues, prod_line::ExprValues, rate::ExprValues, use_mass_kin::Bool)
        sub = add_reactants!(sub_line,1,Vector{ReactantStruct}(undef,0))
        prod = add_reactants!(prod_line,1,Vector{ReactantStruct}(undef,0))
        ns = netstoich(sub, prod)

        rate_DE = isempty(sub) ? rate : mass_rate_DE(sub, use_mass_kin, rate)
        rate_SSA =  isempty(sub) ? rate : mass_rate_SSA(sub, use_mass_kin, rate)
        new(sub, prod, ns, rate, rate_DE, rate_SSA, [], use_mass_kin)
    end
    function ReactionStruct(r::ReactionStruct, syms::Vector{Symbol})
        deps = unique!(recursive_content(r.rate_DE,syms,Vector{Symbol}()))
        is_ma = r.is_pure_mass_action && (length(recursive_content(r.rate_org,syms,Vector{Symbol}()))==0)
        new(r.substrates, r.products, r.netstoich, r.rate_org, r.rate_DE, r.rate_SSA, deps, is_ma)
    end
end

# calculates the net stoichiometry of a reaction
function netstoich(substrates::Vector{ReactantStruct}, products::Vector{ReactantStruct})
    nsdict = Dict{Symbol,Int}(s.reactant => -s.stoichiometry for s in substrates)

    for prod in products
        rsym = prod.reactant
        rcoef = prod.stoichiometry
        @inbounds nsdict[rsym] = haskey(nsdict, rsym) ? nsdict[rsym] + rcoef : rcoef
    end

    net_stoich = Vector{ReactantStruct}()
    for (specsym,coef) in nsdict
        (coef != zero(coef)) && push!(net_stoich, ReactantStruct(specsym,coef))
    end

    net_stoich
end

#Calculates the rate used by ODEs and SDEs. If we want to use masskinetics we have to include substrate concentration, taking higher order terms into account.
function mass_rate_DE(substrates::Vector{ReactantStruct}, use_mass_kin::Bool, old_rate::ExprValues)
    rate = Expr(:call, :*, old_rate)
    coef = one(typeof(substrates[1].stoichiometry))
    if use_mass_kin
        for sub in substrates
            stoich = sub.stoichiometry
            coef *= factorial(stoich)
            push!(rate.args, stoich > one(stoich) ? Expr(:call, :^, sub.reactant, stoich) : sub.reactant)
        end
        (coef > one(coef)) && (rate.args[2] = :($old_rate / $coef))
    end
    return rate
end

#Calculates the rate used by SSAs. If we want to use masskinetics we have to include substrate concentration, taking higher order terms into account.
function mass_rate_SSA(substrates::Vector{ReactantStruct}, use_mass_kin::Bool, old_rate::ExprValues)
    rate = Expr(:call, :*, old_rate)
    use_mass_kin && foreach(sub -> push!(rate.args, :(binomial($(sub.reactant),$(sub.stoichiometry)))), substrates)
    return rate
end

#Takes a reaction line and creates reactions from it and pushes those to the reaction array. Used to creat multiple reactions from e.g. 1.0, (X,Y) --> 0.
function push_reactions(reactions::Vector{ReactionStruct}, sub_line::ExprValues, prod_line::ExprValues, rate::ExprValues, use_mass_kin::Bool)
    lengs = [tup_leng(sub_line), tup_leng(prod_line), tup_leng(rate)]
    (count(lengs.==1) + count(lengs.==maximum(lengs)) < 3) && (throw("malformed reaction"))
    for i = 1:maximum(lengs)
        push!(reactions, ReactionStruct(get_tup_arg(sub_line,i), get_tup_arg(prod_line,i), get_tup_arg(rate,i), use_mass_kin))
    end
end

#Recursive function that loops through the reactants in an reaction line and finds the reactants and their stoichiometry. Recursion makes it able to handle e.g. 2(X+Y+3(Z+XY)) (probably one will not need it though).
function add_reactants!(ex::ExprValues, mult::Int, reactants::Vector{ReactantStruct})
    if typeof(ex)!=Expr
        (ex == 0 || in(ex,empty_set)) && (return reactants)
        if in(ex, getfield.(reactants,:reactant))
            idx = findall(x -> x==ex ,getfield.(reactants,:reactant))[1]
            reactants[idx] = ReactantStruct(ex,mult+reactants[idx].stoichiometry)
        else
            push!(reactants, ReactantStruct(ex,mult))
        end
    elseif ex.args[1] == :*
        add_reactants!(ex.args[3],mult*ex.args[2],reactants)
    elseif ex.args[1] == :+
        for i = 2:length(ex.args)
            add_reactants!(ex.args[i],mult,reactants)
        end
    else
        throw("malformed reaction")
    end
    return reactants
end

#For each reaction, sets its dependencies and whenever it is a pure mass action reaction.
function update_reaction_info(reactions::Vector{ReactionStruct},syms::Vector{Symbol})
    for i = 1:length(reactions)
        reactions[i] = ReactionStruct(reactions[i],syms)
    end
end

#From the vector with all reactions, generates a dictionary with all reactants. Each reactant will point to a number so that X --> means X will be replaced with u[1] in the equations.
function get_reactants(reactions::Vector{ReactionStruct})
    reactants = OrderedDict{Symbol,Int}()
    r_count = 0    ::Int
    for reaction in reactions
        for sub in reaction.substrates
            (!haskey(reactants,sub.reactant)) && (reactants[sub.reactant] = r_count += 1)
        end
        for prod in reaction.products
            (!haskey(reactants,prod.reactant)) && (reactants[prod.reactant] = r_count += 1)
        end
    end
    return reactants
end

#Generates a dictionary with all parameters.
function get_parameters(p)
    parameters = OrderedDict{Symbol,Int}()
    p_count = 0    ::Int
    for parameter in p
        (!haskey(parameters,parameter)) && (parameters[parameter] = p_count += 1)
    end
    return parameters
end

#Produces an array of expressions. Each entry corresponds to a line in the function f, which constitutes the deterministic part of the system. The Expressions can be used for debugging, making LaTex code, or creating the real f function for simulating the network.
function get_f(reactions::Vector{ReactionStruct}, reactants::OrderedDict{Symbol,Int})
    f = Vector{Expr}(undef,length(reactants))
    @inbounds for i = 1:length(f)
        f[i] = :(internal_var___du[$i] = $(Expr(:call, :+)))
    end
    @inbounds for reaction in reactions
        rl = recursive_clean!(reaction.rate_DE)
        @inbounds for r in reaction.netstoich
            sidx = reactants[r.reactant]
            scoef = r.stoichiometry
            ex = recursive_clean!(:($scoef * $(deepcopy(reaction.rate_DE))))
            !(ex isa Number && iszero(ex)) && push!(f[sidx].args[2].args, ex)
        end
    end

    @inbounds for line in f
        if length(line.args[2].args) == 1
            @assert line.args[2].args[1] == :+
            line.args[2] = 0
        else
            line.args[2] = clean_subtractions(line.args[2])
        end
    end

    return f
end

#Produces an array of expressions. Each entry corresponds to a line in the function g, which constitutes the stochastic part of the system. Uses the Guillespie Approach for creating Langevin equations.  The Expressions can be used for debugging, making LaTex code, or creating the real f function for simulating the network.
function get_g(reactions::Vector{ReactionStruct}, reactants::OrderedDict{Symbol,Int}, scale_noise::Symbol)
    numspec = length(reactants)
    numrxs = length(reactions)
    g = Vector{Expr}(undef,numspec*numrxs)
    idx = 0
    @inbounds for k = 1:numrxs
        # initialize to zero for all reactions
        @inbounds for (ssym,sidx) in reactants
            g[idx += 1] = :(internal_var___du[$sidx,$k] = 0)
        end

        ns = reactions[k].netstoich
        @inbounds for r in ns
            sidx = reactants[r.reactant]
            scoef = r.stoichiometry
            ex = recursive_clean!( :($scoef * $scale_noise * sqrt(abs($(deepcopy(reactions[k].rate_DE))))) )
            g[sidx + (k-1)*numspec] = :(internal_var___du[$sidx,$k] = $ex)
        end
    end

    return g
end

#Makes the various jacobian elements required.
function get_jacs(f_rhs::Vector{ExprValues}, syms::Vector{Symbol}, reactions::Vector{ReactionStruct},
                                                                   reactants::OrderedDict{Symbol,Int},
                                                                   parameters::OrderedDict{Symbol,Int};
                                                                   build_jac = true,
                                                                   build_paramjac = true,
                                                                   sparse_jac = false,
                                                                   zeroout_jac = false)
    # jacobian logic
    if build_jac
        if sparse_jac
            jac_prototype, jac = calculate_sparse_jac(reactions, reactants, parameters)
            symjac = nothing
        else
            symjac = calculate_symjac(reactions, reactants)
            jac = calculate_jac(deepcopy(symjac), reactants, parameters, zeroout_jac)
            jac_prototype = nothing
        end
    else
        jac_prototype = jac = symjac = nothing
    end

    paramjac = build_paramjac ? calculate_paramjac(deepcopy(f_rhs), reactants, parameters) : nothing

    if isnothing(symjac)
        return (symjac, jac, jac_prototype, paramjac)
    else
        return (Expr(:quote, symjac), jac, jac_prototype, paramjac)
    end
end

# assumes subs is non-empty!
function massaction_ratelaw_deriv(spec::Symbol, scoef::Int, rate::ExprValues, subs::Vector{ReactantStruct})
    # only case where we should get a constant derivative
    ((length(subs) == 1) && (subs[1].stoichiometry == one(scoef))) && (return rate)

    # otherwise have non-constant derivative
    rl = Expr(:call, :*, rate)
    denom = one(scoef)
    for sub in subs
        name = sub.reactant
        stoich = sub.stoichiometry
        if name == spec
            expon = stoich - one(stoich)
            denom *= factorial(expon)
            if expon == one(expon)
                push!(rl.args, name)
            elseif expon > one(expon)
                push!(rl.args, Expr(:call, :^, name, expon))
            end
        else
            denom *= factorial(stoich)
            push!(rl.args, stoich > one(stoich) ? Expr(:call, :^, name, stoich) : name)
        end
    end
    (denom != one(denom)) && (rl.args[2] = :($(rl.args[2])/$denom))

    rl
end

# given an expression representing a running sum, a ratelaw like expression,
# and a stoichiometric coefficient, update the running sum with the coef
# times the rate law. specialized for mass action reactions
@inline function addstoich_to_exprsum(ex::ExprValues, rl::ExprValues, scoef)
    if (ex isa Number) && iszero(ex)
        if scoef < zero(scoef)
            newex = (scoef == -one(scoef)) ? :(-$rl) : :($scoef * $rl)
        else
            newex = (scoef == one(scoef)) ? rl : :($scoef * $rl)
        end
    else
        if scoef < zero(scoef)
            scoef = -scoef
            newex = (scoef == one(scoef)) ? rl : :($scoef * $rl)
            newex = :($ex - $newex)
        else
            newex = (scoef == one(scoef)) ? rl : :($scoef * $rl)
            newex = :($ex + $newex)
        end
    end
    newex
end

# given an expression representing a running sum, a ratelaw like expression,
# and a stoichiometric coefficient, update the running sum with the coef
# times the rate law. Specialized for SymEngine ratelaws.
@inline function addstoich_to_exprsum(ex::ExprValues, rl::SymEngine.Basic, scoef)
    if (ex isa Number) && iszero(ex)
        newex = Meta.parse(string(scoef*rl))
    else
        if scoef < zero(scoef)
            newex = :($ex - $(Meta.parse(string(-scoef*rl))))
        else
            newex = :($ex + $(Meta.parse(string(scoef*rl))))
        end
    end
    newex
end

# get the correct index type for a dict mapping (i,j) to ExprValues,
# or for a Matrix of ExprValues
@inline getkeytype(d::Dict)                       = keytype(d)
@inline getkeytype(d::Matrix)                     = CartesianIndex{2}
@inline resolveindex(T::Type,i,j)                 = T(i,j)
@inline resolveindex(T::Type{Tuple{Int,Int}},i,j) = (i,j)
@inline checkforkey(d::Dict, key)                 = haskey(d, key)
@inline checkforkey(d::Matrix, key)               = true

# generate Jacobian. Using preceding functions it supports jacexprs as
# a Dict mapping (i,j) => ExprValues or as a dense matrix of ExprValues
function jac_as_exprvalues!(jacexprs, reactions, reactants)

   @inbounds for rx in reactions
        if rx.is_pure_mass_action
            for sub in rx.substrates
                spec  = sub.reactant
                scoef = sub.stoichiometry
                @inbounds j = reactants[spec]

                # derivative with respect to this species
                dratelaw = massaction_ratelaw_deriv(spec, scoef, rx.rate_org, rx.substrates)

                # determine stoichiometric coefficent to multiply by
                @inbounds for ns in rx.netstoich
                    key = resolveindex(getkeytype(jacexprs), reactants[ns.reactant], j)
                    if checkforkey(jacexprs, key)
                        jacexprs[key] = addstoich_to_exprsum(jacexprs[key], deepcopy(dratelaw), ns.stoichiometry)
                    else
                        jacexprs[key] = addstoich_to_exprsum(0, deepcopy(dratelaw), ns.stoichiometry)
                    end
                end
            end
        else
            ratelaw = SymEngine.Basic(recursive_clean!(deepcopy(rx.rate_DE)))
            @inbounds for dep in rx.dependants
                dratelaw = diff(ratelaw, dep)
                j = reactants[dep]

                @inbounds for ns in rx.netstoich
                    key = resolveindex(getkeytype(jacexprs), reactants[ns.reactant], j)
                    if checkforkey(jacexprs, key)
                        jacexprs[key] = addstoich_to_exprsum(jacexprs[key], dratelaw, ns.stoichiometry)
                    else
                        jacexprs[key] = addstoich_to_exprsum(0, dratelaw, ns.stoichiometry)
                    end
                end
            end
        end
    end
    nothing
end

# generate a dense matrix of ExprValues representing the Jacobian
function calculate_symjac(reactions, reactants)
    nspecs   = length(reactants)
    jacexprs = Matrix{ExprValues}(undef, nspecs, nspecs)
    fill!(jacexprs, 0)
    jac_as_exprvalues!(jacexprs, reactions, reactants)
    jacexprs
end

# Generate the Jacobian evaluation function as a single expression.
# zeroout_jac kwarg controls whether to first fill the matrix with zeros and then only fill in non-zero entries
# vs looping through and putting expressions for entries that are always zero.
# This is needed for sufficiently large dense matrices to allow compilation. (e.g. 1000x1000)
function calculate_jac(symjac::Matrix{ExprValues}, reactants::OrderedDict{Symbol,Int}, parameters::OrderedDict{Symbol,Int}, zeroout_jac=false)
    func_body = Expr(:block)
    zeroout_jac && push!(func_body.args, :(fill!(internal___var___J, zero(eltype(internal___var___J)))))
    @inbounds for j = 1:size(symjac)[2], i = 1:size(symjac)[1]
        if symjac[i,j] isa Number
            ex = symjac[i,j]
            (!zeroout_jac || !iszero(ex)) && push!(func_body.args, :(internal___var___J[$i,$j] = $(ex)))
        else
            ex = recursive_replace!(symjac[i,j],(reactants,:internal___var___u), (parameters, :internal___var___p))
            splitplus!(ex)
            push!(func_body.args, :(internal___var___J[$i,$j] = $(ex)))
        end
    end
    push!(func_body.args,:(return nothing))

    return :((internal___var___J,internal___var___u,internal___var___p,t) -> @inbounds $func_body)
end

# convert Dict mapping (i,j) => val to sparse matrix
function dict_to_sparsemat(d; vals=nothing)
    V   = isnothing(vals) ? collect(values(d)) : vals
    len = length(d)
    I   = Vector{Int}(undef,len)
    J   = Vector{Int}(undef,len)
    for (i,k) in enumerate(keys(d))
        I[i] = k[1]
        J[i] = k[2]
    end
    return sparse(I,J,V)
end

# create a sparse Jacobian
function calculate_sparse_jac(reactions, reactants, parameters)

    # get the elements and structure of sparse matrix
    jacexprs = Dict{Tuple{Int,Int},ExprValues}()
    jac_as_exprvalues!(jacexprs, reactions, reactants)
    jac_prototype = dict_to_sparsemat(jacexprs, vals=ones(length(jacexprs)))

    # build the function
    jfun = Expr(:block)
    nspecs = length(reactants)
    push!(jfun.args, :(internal___var___Jvals = nonzeros(internal___var___J)))
    rows = rowvals(jac_prototype)
    @inbounds for j = 1:nspecs
        @inbounds for k in nzrange(jac_prototype,j)
            i = rows[k]
            ex = (jacexprs[(i,j)] isa Number) ? jacexprs[(i,j)] : recursive_replace!(jacexprs[(i,j)],(reactants,:internal___var___u), (parameters,:internal___var___p))
            push!(jfun.args, :(internal___var___Jvals[$k] = $ex))
        end
    end
    push!(jfun.args, :(return nothing))
    jfunex = :((internal___var___J,internal___var___u,internal___var___p,t) -> @inbounds $jfun)

    return jac_prototype, jfunex
end

#Makes the Jacobian, with respect to parameter values.
function calculate_paramjac(f_rhs::Vector{ExprValues}, reactants::OrderedDict{Symbol,Int}, parameters::OrderedDict{Symbol,Int})
    func_body = Expr(:block)
    for j = 1:length(parameters), i = 1:length(reactants)
        paramjac_entry = Meta.parse(string(diff(SymEngine.Basic(f_rhs[i]), parameters.keys[j])))
        push!(func_body.args, :(internal___var___pJ[$i,$j] = $(recursive_replace!(paramjac_entry,(reactants,:internal___var___u), (parameters, :internal___var___p)))))
    end
    push!(func_body.args,:(return internal___var___pJ))
    return :((internal___var___pJ,internal___var___u,internal___var___p,t) -> @inbounds $func_body)
end

### Jump related functions ###

#Creates expressions for jump affects and rates. Also creates and array with MassAction, ConstantRate and VariableRate Jumps.
function get_jumps(reactions::Vector{ReactionStruct}, reactants::OrderedDict{Symbol,Int}, parameters::OrderedDict{Symbol,Int}; minimal_jumps=false)
    rates = Vector{ExprValues}()
    affects = Vector{Vector{Expr}}()
    jumps = Expr(:tuple)
    idx = 0
    for reaction in reactions
        hastime = recursive_contains(:t,reaction.rate_SSA)
        if !minimal_jumps || !reaction.is_pure_mass_action || hastime
            idx += 1
            push!(rates, recursive_clean!(deepcopy(reaction.rate_SSA)))
            push!(affects, Expr[:(integrator.u[$(reactants[r.reactant])] += $(r.stoichiometry)) for r in reaction.netstoich])
            syntax_rate = recursive_replace!(deepcopy(rates[idx]), (reactants, :internal_var___u), (parameters, :internal_var___p))
            hastime ? push!(jumps.args,Expr(:call,:VariableRateJump)) : push!(jumps.args,Expr(:call,:ConstantRateJump))
            push!(jumps.args[idx].args, :((internal_var___u,internal_var___p,t) -> @inbounds $syntax_rate))
            push!(jumps.args[idx].args, :(integrator -> @inbounds $(expr_arr_to_block(deepcopy(affects[idx])))))
        end
    end
    return (Tuple(rates),Tuple(affects),jumps)
end

# generates regular jumps representation
function get_regularjumps(reactions::Vector{ReactionStruct}, reactants::OrderedDict{Symbol,Int}, parameters::OrderedDict{Symbol,Int}; minimal_jumps=false)
    reg_rates = Expr(:block)
    reg_c = Expr(:block)
    @inbounds for (idx,reaction) in enumerate(reactions)
        rate = recursive_clean!(deepcopy(reaction.rate_SSA))
        syntax_rate = recursive_replace!(rate, (reactants,:internal_var___u), (parameters, :internal_var___p))
        push!(reg_rates.args, :(@inbounds internal_var___out[$idx]= $syntax_rate))
        foreach(r -> push!(reg_c.args,:(@inbounds internal_var___dc[$(reactants[r.reactant]),$idx]=$(r.stoichiometry))), reaction.netstoich)
    end
    reg_jumps = :(RegularJump((internal_var___out,internal_var___u,internal_var___p,t)->$reg_rates,(internal_var___dc,internal_var___u,internal_var___p,t,internal_var___mark)->$reg_c,zeros($(length(reactants)),$(length(reactions)));constant_c=true))
    return reg_jumps
end

### Pre Defined Functions that can be inserted into the reaction rates ###

#Hill function made avaiable
hill_name = Set{Symbol}([:hill, :Hill, :h, :H, :HILL])
function hill(expr::Expr)
    return :($(expr.args[3])*($(expr.args[2])^$(expr.args[5]))/($(expr.args[4])^$(expr.args[5])+$(expr.args[2])^$(expr.args[5])))
end
hillR_name = Set{Symbol}([:hill_repressor, :hillr, :hillR, :HillR, :hR, :hR, :Hr, :HR, :HILLR])
function hillR(expr::Expr)
    return :($(expr.args[3])*($(expr.args[4])^$(expr.args[5]))/($(expr.args[4])^$(expr.args[5])+$(expr.args[2])^$(expr.args[5])))
end

#Michaelis menten function made avaiable.
mm_name = Set{Symbol}([:MM, :mm, :Mm, :mM, :M, :m])
function mm(expr::Expr)
    return :($(expr.args[3])*$(expr.args[2])/($(expr.args[4])+$(expr.args[2])))
end
#Michaelis menten function made avaiable.
mmR_name = Set{Symbol}([:mm_repressor, :MMR, :mmr, :mmR, :MmR, :mMr, :MR, :mr, :Mr, :mR])
function mmR(expr::Expr)
    return :($(expr.args[3])*$(expr.args[4])/($(expr.args[4])+$(expr.args[2])))
end

#Allows the user to define new function and enable the @reaction_network macro to see them.
macro reaction_func(expr)
    name = expr.args[1].args[1]
    args = expr.args[1].args[2:end]
    maths = expr.args[2].args[2]

    funcdict[name]  = x -> replace_names(maths, args, x)
end


text:
using Revise

using POMDPs
using POMDPModels
using POMDPModels: LDNormalStateDist, LightDark1DState
using POMDPSimulators
using ParticleFilters
using POMCPOW

m = LightDark1D()
ds0 = POMDPs.initialstate_distribution(m)
s0 = rand(ds0)

up = BootstrapFilter(m, 1000);
pomdp = LightDark1D()
N=1000

up = BootstrapFilter(pomdp, N)
true_b0 = LDNormalStateDist(s0.y, 5.0)
sample_b0 = LDNormalStateDist(s0.y, 10.0)

function state_weight(d::LDNormalStateDist, s::LightDark1DState)
    dist = POMDPModels.Normal(d.mean, d.std)
    return pdf(dist, s.y)
end

function POMCPOW.importance_weight(s::LightDark1DState, rb::LDNormalStateDist, tb::LDNormalStateDist)
    q = state_weight(rb, s)
    p = state_weight(tb, s)
    return p/q
end

solver = POMCPOWSolver(tree_queries=10,
                        check_repeat_obs=false,
                        check_repeat_act=true,
                        k_observation=2,
                        alpha_observation=0.1
                        )
planner = POMDPs.solve(solver, m)

a, info = POMCPOW.action_info(planner, sample_b0, tree_in_info=true, true_b=true_b0)

# println("Entering Simulation...")
# for (s, a, r, bp, t) in stepthrough(m, planner, up, b0, s0, "s,a,r,bp,t", max_steps=50)
#     @show a
#     @show r
#     @show t
# end


text:
__precompile__()
module LShapedSolvers

# Standard library
using LinearAlgebra
using SparseArrays
using Distributed
using Printf

# External libraries
using TraitDispatch
using Parameters
using JuMP
using StochasticPrograms
using MathProgBase
using RecipesBase
using ProgressMeter
using Gurobi
using Gurobi: updatemodel!

import Base: show, put!, wait, isready, take!, fetch, zero, +, length
import StochasticPrograms: StructuredModel, internal_solver, optimize_structured!, fill_solution!, solverstr

const MPB = MathProgBase

export
    LShapedSolver,
    Crash,
    StructuredModel,
    add_params!,
    optimsolver,
    optimize_structured!,
    fill_solution!,
    get_decision,
    get_objective_value,
    LShaped,
    DLShaped,
    Regularized,
    DRegularized,
    TrustRegion,
    DTrustRegion,
    LevelSet,
    DLevelSet

# Include files
include("LQSolver.jl")
include("subproblem.jl")
include("hyperplane.jl")
include("AbstractLShaped.jl")
include("regularization.jl")
include("parallel.jl")
include("solvers/LShaped.jl")
include("solvers/DLShaped.jl")
include("solvers/Regularized.jl")
include("solvers/DRegularized.jl")
include("solvers/TrustRegion.jl")
include("solvers/DTrustRegion.jl")
include("solvers/LevelSet.jl")
include("solvers/DLevelSet.jl")
include("crash.jl")
include("spinterface.jl")

end # module


text:
__precompile__()

module AdaptiveLASSO

    ### exports
    export
    # main function
    LSIAF,
    # input names
    CoordinateDescent3, CoordinateDescentTM

    ### include source files

    # core algorithm
    include("utils.jl")
    include("ComposeMatrices.jl")
    include("CoordinateDescentTM.jl")
    include("CoordinateDescent3.jl")
    include("adaLASSO.jl")
    include("LSIAF.jl")

end


text:
using DataStructures: OrderedDict

@testset "Datasets" begin
  data = [rand(4,5,12), rand(4,5,12), rand(4,5)]
  axlist1 = [RangeAxis("XVals",1.0:4.0), CategoricalAxis("YVals",[1,2,3,4,5]), RangeAxis("Time",Date(2001,1,15):Month(1):Date(2001,12,15))]
  axlist2 = [RangeAxis("XVals",1.0:4.0), CategoricalAxis("YVals",[1,2,3,4,5])]
  props = [Dict("att$i"=>i) for i=1:3]
  c1,c2,c3 = (YAXArray(axlist1, data[1], props[1]),
  YAXArray(axlist1, data[2], props[2]),
  YAXArray(axlist2, data[3], props[3])
  )
  ds = Dataset(avar = c1, something = c2, smaller = c3)
  @testset "Basic functions" begin
    b = IOBuffer()
    show(b,ds)
    s = split(String(take!(b)),"\n")
    s2 = """
    YAXArray Dataset
    Dimensions:
    XVals               Axis with 4 Elements from 1.0 to 4.0
    YVals               Axis with 5 elements: 1 2 3 4 5
    Time                Axis with 12 Elements from 2001-01-15 to 2001-12-15
    Variables: avar something smaller """
    s2 = split(s2,"\n")
#     @test s[[1,2,6]] == s2[[1,2,6]]
#     @test all(i->in(i,s2), s[3:5])
    for n in [:avar, :something, :smaller, :XVals, :Time, :YVals]
      @test n in propertynames(ds)
      @test n in propertynames(ds, true)
    end
    @test :axes ∉ propertynames(ds)
    @test :cubes ∉ propertynames(ds)
    @test :axes ∈ propertynames(ds, true)
    #Test getproperty
    @test all(i->in(i,values(ds.axes)),axlist1)
    @test collect(keys(ds.cubes)) == [:avar, :something, :smaller]
    @test collect(values(ds.cubes)) == [c1,c2,c3]
    @test ds.avar === c1
    @test ds.something === c2
    @test ds.smaller === c3
    @test ds[:avar] === c1
    ds2 = ds[[:avar, :smaller]]
    @test collect(keys(ds2.cubes)) == [:avar, :smaller]
    @test collect(values(ds2.cubes)) == [c1,c3]
    @test YAXArrays.Datasets.fuzzyfind("hal", ["hallo","bla","something"]) == 1
    ds3 = ds[["av", "some"]]
    @test collect(keys(ds3.cubes)) == [:av, :some]
    @test collect(values(ds3.cubes)) == [c1,c2]
    @test ds["avar"] === c1
  end
  @testset "Dataset interface" begin
    struct MockDataset
      vars
      dims
      attrs
      path
    end
    Base.getindex(d::MockDataset,i) = d.vars[i]
    Base.haskey(d::MockDataset,i) = haskey(d.vars,i)
    YAXArrayBase.get_varnames(d::MockDataset) = (keys(d.vars)...,)
    YAXArrayBase.get_var_dims(d::MockDataset,name) = d.dims[name]
    YAXArrayBase.get_var_attrs(d::MockDataset, name) = d.attrs[name]
    YAXArrayBase.allow_missings(d::MockDataset) = !occursin("nomissings",d.path)
    function YAXArrayBase.create_empty(::Type{MockDataset},path)
      mkpath(dirname(path))
      open(_->nothing, path,"w")
      MockDataset(Dict(),Dict(),Dict(),path)
    end
    function YAXArrayBase.add_var(ds::MockDataset, T, name, s, dimlist, atts;kwargs...)
      data = Array{T}(undef,s...)
      ds.vars[name] = data
      ds.dims[name] = dimlist
      ds.attrs[name] = atts
      data
    end
    YAXArrayBase.backendlist[:mock] = MockDataset
    push!(YAXArrayBase.backendregex,r".mock$"=>MockDataset)
    data1,data2,data3,d1,d2,d3 = (rand(12,5,2),rand(12,5),rand(12,5,2),1:12, 0.1:0.1:0.5, ["One","Two"])
    att1 =  Dict("att1"=>5,"_ARRAY_OFFSET"=>(2,0,0))
    att2 =  Dict("att2"=>6,"_ARRAY_OFFSET"=>(2,0))
    attd1 = Dict("_ARRAY_OFFSET"=>2, "units"=>"days since 2001-01-01", "calendar"=>"gregorian")
    attd2 = Dict("attd"=>"d")
    attd3 = Dict("attd"=>"d")
    function MockDataset(p)
        MockDataset(
        OrderedDict("Var1"=>data1, "Var2"=>data2, "Var3"=>data3, "time"=>d1,"d2"=>d2, "d3"=>d3),
        Dict("Var1"=>("time","d2","d3"),"Var2"=>("time","d2"),"Var3"=>("time","d2","d3"),"time"=>("time",),"d2"=>["d2"],"d3"=>["d3"]),
        Dict("Var1"=>att1,"Var2"=>att2,"Var3"=>att1,"time"=>attd1,"d2"=>attd2,"d3"=>attd3),
        p
        )
    end
    m = MockDataset("testpath.mock")
    @testset "collectdims" begin
      dcollect = YAXArrays.Datasets.collectdims(m)
      @test dcollect["time"].ax isa RangeAxis
      @test YAXArrays.Cubes.Axes.axname(dcollect["time"].ax) == "time"
      @test dcollect["time"].ax.values == DateTime(2001,1,4):Day(1):DateTime(2001,1,13)
      @test dcollect["time"].offs == 2
      @test dcollect["d2"].ax isa RangeAxis
      @test YAXArrays.Cubes.Axes.axname(dcollect["d2"].ax) == "d2"
      @test dcollect["d2"].ax.values == 0.1:0.1:0.5
      @test dcollect["d2"].offs == 0
      @test dcollect["d3"].ax isa CategoricalAxis
      @test YAXArrays.Cubes.Axes.axname(dcollect["d3"].ax) == "d3"
      @test dcollect["d3"].ax.values == ["One","Two"]
      @test dcollect["d3"].offs == 0
      a1 = [0.1,0.2,0.3,0.4]
      a2 = [0.1,0.21,0.3,0.4]
      @test YAXArrays.Datasets.testrange(a1)== 0.1:0.1:0.4
      @test YAXArrays.Datasets.testrange(a2) isa Array
      @test YAXArrays.Datasets.testrange(a2) == [0.1,0.21,0.3,0.4]
    end
    @testset "open_dataset" begin
      ds = open_dataset("test.mock")
      @test size(ds.Var1) == (10,5,2)
      @test size(ds.Var2) == (10,5)
      @test all(in(keys(ds.axes)),(:time,:d2,:d3))
      ar = Cube(ds)
      @test ar isa YAXArray
      @test size(ar) == (10,5,2,2)
      @test YAXArrays.Cubes.Axes.axname.(ar.axes) == ["time","d2","d3","Variable"]
      @test ar.axes[4].values == ["Var1","Var3"]
    end
    @testset "Dataset creation" begin
      al = [RangeAxis("Time",Date(2001):Month(1):Date(2001,12,31)), CategoricalAxis("Variable",["A","B"]), RangeAxis("Xvals",1:10)]
      #Basic
      newds, newds2 = YAXArrays.Datasets.createdataset(MockDataset,al)
      @test YAXArrays.Cubes.axsym.(newds2.axes) == [:Time, :Xvals, :Variable]
      @test newds2.axes[1].values == Date(2001):Month(1):Date(2001,12,31)
      @test newds2.axes[3].values == ["A","B"]
      @test newds2.axes[2].values == 1:10
      @test newds2.data isa YAXArrays.Cubes.DiskArrayTools.DiskArrayStack
      # A bit more advanced
      fn = string(tempname(),".mock")
      newds, newds2 = YAXArrays.Datasets.createdataset(MockDataset,al,path = fn, persist = false,
      chunksize = (4,2,4), chunkoffset = (2,0,3), properties = Dict("att1"=>5), datasetaxis="A")
      @test size(newds.data) == (12,2,10)
      @test size(newds.data.parent) == (14,2,13)
      @test eltype(newds.data) <: Union{Float32,Missing}
      @test newds.properties["att1"] == 5
      @test isfile(fn)
      newds = nothing
      newds2 = nothing
      # Without missings
      fn = string(tempname(),"nomissings.mock")
      newds = YAXArrays.Datasets.createdataset(MockDataset,al,path = fn,datasetaxis="A")
    end
  end
end


text:
export impedanceM1, impedanceM2, impedanceM3, impedanceM4,
        impedanceM5, impedanceM6, impedanceM7, impedanceM8,
        impedanceM9, impedanceM9b

include("constant.jl")
#using .constant

"""

        impedanceM1(f,σ)

Delany and Bazley model

Compute characteristic ground impedance with 1 input parameter

# Arguments
- `σ`: Pa.s.m-2, flow resitivity (realistic values [5x10^3 3x10^7])

# Example:
```julia-repl
julia> impedanceM1(100, 200_000)
```
---
[1]. Denaly & Bazley, Acoustical properties of fibrous absorbent materials,
Appl. Acoustisc, 1970

[2]. Attenborough, Keith, Imran Bashir, and Shahram Taherzadeh. "Outdoor ground impedance models."
The Journal of the Acoustical Society of America 129.5 (2011): 2806-2819.

"""
function impedanceM1(f,σ)
    X = f*ϱ0/σ # in Ref. [2], this equation is typo error
    Zc = 1 + 0.0571*X^(-0.754) + 1im*0.087*X^(-0.732) # Eq. (1b) Ref. [2]
    return Zc
end

"""

        impedanceM2(f,σ)

        impedanceM2(f,σ,α)

Variable prositty models

Compute characteristic ground impedance with 2 input parameters

# Arguments

- `σ`: Pa.s.m-2, flow resitivity (realistic values [5x10^3 3x10^7])

- `α`: /m, rate of porosity variation (realistic values [1 10]), defaul α = 8 m-1

# Example:
```julia-repl
julia> impedanceM2(100,200_000)
julia> impedanceM2(100, 200_000,5)
```

---
[1]. Attenborough et al., Outdoor ground impedance models EuroNoise, 2015

"""
function impedanceM2(f,σ,α::Any=8)
    # default value for α = 8m-1 as used for institutional grass (Donato, 1977, JASA)
    Zc =    (1+1im)/sqrt(π*γa*ϱ0)*sqrt(σ/f) +
            1im*c0*α/(8*π*f)
    return Zc
end

"""

        impedanceM3(f,σ,Ω)

        impedanceM3(f,σ,Ω,svor,sent)

Wilson relaxation model

Compute characteristic ground impedance with 2 input parameters

# Arguments
- `σ`: Pa.s.m-2, flow resitivity (realistic values [5x10^3 3x10^7])

- `Ω`: prosity (realistic values [0.1 0.8]) (Table 10.1, ref. [1])

- `svor`, `sent`: shape factor, default values = 1

# Example:
```julia-repl
julia> impedanceM3(100,200_000,0.5)
```

---
[1]. Ostashev, Vladimir E., and D. Keith Wilson.
Acoustics in moving inhomogeneous media. CRC Press, 2015.

"""
function impedanceM3(f,σ,Ω, svor::Any=1,sent::Any=1)
    ω = 2π*f
    q = Ω^-0.25 # page 371
    τ0 = (2*ϱ0*q^2)/(σ*Ω) # Eq 10.53

    term1 = 2im/(ω*τ0) + 1 + 1/sqrt(1-1im*svor^2*ω*τ0)
    term2 = 1 + (γa-1)/sqrt(1 - 1im*sent^2*Npr*ω*τ0)
    Zc = (q/Ω)*(term1/term2)^0.5 # Eq 10.54

    return Zc
end

"""

        impedanceM4(f,σ,Ω)

        impedanceM4(f,σ,Ω,q)

Zwikker and Kosten model

Compute characteristic ground impedance with 3 input parameters

# Arguments

- `σ`: Pa.s.m-2, flow resitivity (realistic values [5x10^3 3x10^7])

- `Ω`: prosity (realistic values [0.1 0.8])

- `q`: tortuosity (realistic values [1.8 1.0]), default `q` = 1.0

# Example:
```julia-repl
julia> impedanceM4(100,200_000,0.5)
julia> impedanceM4(100,200_000,0.5,1.0)
```

---
[1]. Zwikker, Cornelis, and Cornelis Willem Kosten. Sound absorbing materials. Elsevier publishing company, 1949.

[2]. Attenborough, Keith, Imran Bashir, and Shahram Taherzadeh. "Outdoor ground impedance models."
The Journal of the Acoustical Society of America 129.5 (2011): 2806-2819.

"""
function impedanceM4(f,σ,Ω,q::Any=1.0)
    # the default value for q is for soft soil
    T = q*q
    ω = 2π*f

    Zc = (1/Ω) * ( T + 1im*Ω*σ/(ω*ϱ0) )^0.5 # Eq. 2a (ref. [2])
    return Zc
end

"""

        impedanceM5(f,σ,Ω)

        impedanceM5(f,σ,Ω,q)

Taraldsen model

Compute characteristic ground impedance with 3 input parameters

# Arguments

- `σ`: Pa.s.m-2, flow resitivity (realistic values [5x10^3 3x10^7])

- `Ω``: prosity (realistic values [0.1 0.8])

- `q`: tortuosity (realistic values [1.8 1.0]), default `q` = 1.0

# Example:
```julia-repl
julia> impedanceM5(100,200_000,0.5)
julia> impedanceM5(100,200_000,0.5,1.0)
```

---
[1]. Taraldsen, G. "The Delany-Bazley impedance model and Darcy's law."
Acta Acustica united with Acustica 91.1 (2005): 41-50.

"""
function impedanceM5(f,σ,Ω,q::Any=1.0)
    β✶ = sqrt(γa)*Ω/q
    f✶ = 1/(2π*ϱ0)*(Ω/q)*(σ/q) # Eq. 14 note √T = q
    β = β✶*(1+1im*f✶/f)^-0.5
    Zc = 1/β

    return Zc
end

"""

        impedanceM6(f,σ,Ω)

        impedanceM6(f,σ,Ω,q)

Hamet phenomenological models

Compute characteristic ground impedance with 3 input parameters

# Arguments

- `σ`: Pa.s.m-2, flow resitivity (realistic values [5x10^3 3x10^7])

- `Ω`: prosity (realistic values [0.1 0.8])

- `q`: tortuosity (realistic values [1.8 1.0]), default `q` = 1.0

# Example:
```julia-repl
julia> impedanceM6(100,200_000,0.5)
julia> impedanceM6(100,200_000,0.5,1.0)
```

---
[1]. Bérengier, M. C., et al. "Porous road pavements: Acoustical characterization and propagation effects."
The Journal of the Acoustical Society of America 101.1 (1997): 155-162.
"""
function impedanceM6(f,σ,Ω,q::Any=1.0)

    q² = q*q
    fμ = Ω*σ/(2π*ϱ0*q²) # Eq. 19
    fθ = σ/(2π*ϱ0*Npr)  # Eq. 20
    Fμ = 1+1im*fμ/f     # Eq. 21
    Fθ = 1 + 1im*fθ/f   # Eq. 22
    Zc = (q/Ω)*Fμ^0.5*( γa - (γa-1)/Fθ )^-0.5   # Eq. 23, without ρ0c term

    return Zc
end

"""

        impedanceM7(f,σ,Ω)

        impedanceM7(f,σ,Ω,q)

Identical tortuous pores model

Compute characteristic ground impedance with 3 input parameters

# Arguments

- `σ`: Pa s m-2, flow resitivity (realistic values [5x10^3 3x10^7])

- `Ω`: prosity (realistic values [0.1 0.8])

- `q`: tortuosity (realistic values [1.8 1.0]), default `q` = 1.0

# Example:
```julia-repl
julia> impedanceM7(100,200_000,0.5)
julia> impedanceM7(100,200_000,0.5,1.0)
```

---
[1]. Attenborough, Keith, Imran Bashir, and Shahram Taherzadeh. "Outdoor ground impedance models."
The Journal of the Acoustical Society of America 129.5 (2011): 2806-2819.

"""
function impedanceM7(f,σ,Ω,q::Any=1.0)
    ω = 2π*f
    q² = q*q # or = T
    GS(x) = 1-tanh(x*sqrt(-1im)) / (x*sqrt(-1im)) # Eq. 6e

    λ = sqrt((3*ϱ0*ω*q²)/(Ω*σ)) # Eq. 6f
    ρ_λ = ϱ0/GS(λ) # Eq. 6c
    C_λ = 1/(γa*P0) * ( γa - (γa-1)*GS(λ*sqrt(Npr)) ) # Eq. 6d

    Zc = (ϱ0*c0)^-1 * ( (q²/Ω^2)*ρ_λ/C_λ )^0.5 # Eq. 6b

    return Zc
end

"""

        impedanceM8(f,σ,Ω)

        impedanceM8(f,σ,Ω,sf)

        impedanceM8(f,σ,Ω,sf,q)

Attenborough four parameter model

Compute characteristic ground impedance with 4 input parameters

# Arguments

- `σ`: Pa s m-2, flow resitivity (realistic values [5x10^3 3x10^7])

- `Ω`: prosity (realistic values [0.1 0.8])

- `q`: tortuosity (realistic values [1.8 1.0]), default `q` = 1.0

- `sf`: pore shape factor (realistic values [0.6 1]), default `sf` = 0.75

# Example:
```julia-repl
julia> impedanceM8(100,200_000,0.5)
julia> impedanceM8(100,200_000,0.5,0.75)
julia> impedanceM8(100,200_000,0.5,0.75,1.0)
```

---
[1]. Attenborough, Keith. "Acoustical impedance models for outdoor ground surfaces."
Journal of Sound and Vibration 99.4 (1985): 521-544.

"""
function impedanceM8(f,σ,Ω,sf::Any=0.75,q::Any=1.0)
    ω = 2π*f
    q² = q*q

    kb =    (γa*Ω)^0.5 *
            ( (4/3 - (γa-1)/γa*Npr)*q²/Ω + 1im*sf^2*σ/(ω*ϱ0) )^0.5 # Eq. 12

    Zc = ( 4*q²/(3*Ω) + 1im*sf^2*σ/(ω*ϱ0) ) / kb # Eq. 11

    return Zc
end

"""

        impedanceM9(f,s̄,Ω)

        impedanceM9(f,s̄,Ω,σs)

Horoshenkov porous media model

Compute characteristic ground impedance with 3 input parameters

# Arguments

- `s̄`: median pore size (realistic values [2 2000 μm]), this parameter linked
to flow resistivity

- `Ω`: prosity (equivalent to Ω of other models) realistics values [0.1 0.8]

- `σs`: standard deviation of pore size distribution (realistic values [0 0.5]),
default `σs` = 0.3

# Example:
```julia-repl
julia> impedanceM9(100,200*10^-6,0.5)
julia> impedanceM9(100,200*10^-6,0.5,0.3)
```

---
[1]. Horoshenkov, Kirill V., Alistair Hurrell, and Jean-Philippe Groby. "A three-parameter analytical model for the acoustical properties of porous media."
The Journal of the Acoustical Society of America 145.4 (2019): 2512-2517.
"""
function impedanceM9(f,s̄,Ω,σs::Any=0.3)
    𝜙 = Ω # for consistent with original notations (ref. [1])

    ω = 2π*f
    α∞ = exp( 4*(σs*log(2))^2 ) # tourtuosity See Eq. 12

    θρ1 = 1/3
    θρ2 = exp( -1/2*(σs*log(2))^2 ) / sqrt(2)
    θρ3 = θρ1/θρ2

    σ = 8*η0*α∞ / (s̄^2*𝜙) * exp( 6*(σs*log(2))^2 )    # Eq. 15
    σ′ = 8*η0*α∞ / (s̄^2*𝜙) * exp( -6*(σs*log(2))^2 )  # Eq. 18

    ϵρ = sqrt( -1im*ω*ϱ0*α∞/(𝜙*σ) )
    F̃ρ = (1 + θρ3*ϵρ + θρ1*ϵρ^2) / (1 + θρ3*ϵρ)    # Eq. 14

    θc1 = 1/3
    θc2 = exp( 3/2*(σs*log(2))^2 ) / sqrt(2)
    θc3 = θc1/θc2

    ϵc = sqrt( -1im*ω*ϱ0*Npr*α∞ / (𝜙*σ′) )
    F̃c = (1 + θc3*ϵc + θc1*ϵc^2) / (1+ θc3*ϵc)     # Eq. 17

    C̃ω = 𝜙/(γa*P0) * ( γa - (γa-1)/(1+ϵc^-2*F̃c) )    # Eq 16
    ρ̃ω = ϱ0*α∞/𝜙 * (1 + ϵρ^-2*F̃ρ)    # Eq. 13

    Zc = 1/(ϱ0*c0) * sqrt( ρ̃ω / C̃ω)   # Eq. 19 characteristic impedance

    return Zc
end

"""

        impedanceM9(f,σ,Ω)

        impedanceM9(f,σ,Ω,σs)

Horoshenkov porous media model (σ input)

Compute characteristic ground impedance with 3 input parameters

# Arguments

- `σ`: Pa s m-2, flow resitivity (realistic values [5x10^3 3x10^7])

- `Ω`: prosity (equivalent to Ω of other models) realistics values [0.1 0.8]

- `σs`: standard deviation of pore size distribution (realistic values [0 0.5]),
default `σs` = 0.3

# Example:
```julia-repl
julia> impedanceM9(100,200_000,0.5)
julia> impedanceM9(100,200_000,0.5,0.3)
```

---
[1]. Horoshenkov, Kirill V., Alistair Hurrell, and Jean-Philippe Groby. "A three-parameter analytical model for the acoustical properties of porous media."
The Journal of the Acoustical Society of America 145.4 (2019): 2512-2517.
"""
function impedanceM9b(f,σ,Ω,σs::Any=0.3)
    𝜙 = Ω # for consistent with original notations (ref. [1])

    ω = 2π*f
    α∞ = exp( 4*(σs*log(2))^2 ) # tourtuosity See Eq. 12

    θρ1 = 1/3
    θρ2 = exp( -1/2*(σs*log(2))^2 ) / sqrt(2)
    θρ3 = θρ1/θρ2

    # convert σ to s̄ using equation 15
    s̄ = ( 8*η0*α∞ / (σ*𝜙) * exp( 6*(σs*log(2))^2 ) )^0.5

    #σ = 8*η0*α∞ / (s̄^2*𝜙) * exp( 6*(σs*log(2))^2 )    # Eq. 15
    σ′ = 8*η0*α∞ / (s̄^2*𝜙) * exp( -6*(σs*log(2))^2 )  # Eq. 18

    ϵρ = sqrt( -1im*ω*ϱ0*α∞/(𝜙*σ) )
    F̃ρ = (1 + θρ3*ϵρ + θρ1*ϵρ^2) / (1 + θρ3*ϵρ)    # Eq. 14

    θc1 = 1/3
    θc2 = exp( 3/2*(σs*log(2))^2 ) / sqrt(2)
    θc3 = θc1/θc2

    ϵc = sqrt( -1im*ω*ϱ0*Npr*α∞ / (𝜙*σ′) )
    F̃c = (1 + θc3*ϵc + θc1*ϵc^2) / (1+ θc3*ϵc)     # Eq. 17

    C̃ω = 𝜙/(γa*P0) * ( γa - (γa-1)/(1+ϵc^-2*F̃c) )    # Eq 16
    ρ̃ω = ϱ0*α∞/𝜙 * (1 + ϵρ^-2*F̃ρ)    # Eq. 13

    Zc = 1/(ϱ0*c0) * sqrt( ρ̃ω / C̃ω)   # Eq. 19 characteristic impedance

    return Zc
end


text:
mutable struct Person
    name::String
    age::Int
    print::Function # 멤버함수 필드

    function Person(name::String, age::Int)
        this = new() # 객체 생성
        this.name = name
        this.age = age
        # 멤버함수 삽입
        this.print=function() Person__print(this) end
        this # 반환
    end
end

# 멤버함수 구현
function Person__print(p::Person)
    println("이름:$(p.name), 나이:$(p.age)")
end

p = Person("홍길동", 22) # 객체 생성
p.print()

max_stars_repo_name: myyrakle/Julia-Example

max_issues_repo_name: myyrakle/Julia-Example

max_forks_repo_name: myyrakle/Julia-Example


text:
using MshReader
using WriteVTK

# runs MshFileReader
nodesCoordMat, connectivity, physicalNames, elemPhysNums = MshFileReader( joinpath("examples","cube.msh") )

#
print( "\n-------------------------\nexporting matrices")
print( "\nnodes Coords Mat", nodesCoordMat)
print( "\nconnec", connectivity )
print( "\nelem phys nums", elemPhysNums )

#
print( "\n-------------------------\nexporting to VTK...\n")
cells = Vector{MeshCell}( undef, length( connectivity ) )
MEBIs = Vector{String}(   undef, length( connectivity ) )

for i in 1:length(connectivity)
    aux = connectivity[i]
    print("  i: ", i, " aux:", aux ,"\n" )
    if length( aux ) == 1
        cells[i] = MeshCell( VTKCellTypes.VTK_VERTEX, connectivity[i] )
    elseif length( aux ) == 3
        cells[i] = MeshCell( VTKCellTypes.VTK_TRIANGLE, connectivity[i] )
    elseif length( aux ) == 4
        cells[i] = MeshCell( VTKCellTypes.VTK_TETRA, connectivity[i] )
    end

    if elemPhysNums[i] == 0
        MEBIs[i] = "null"
    else
        MEBIs[i] = physicalNames[ elemPhysNums[i] ]
    end
end

vtk_grid("cube.vtu", Array( nodesCoordMat' ) , cells, ascii=true) do vtk
    vtk["MEBI_params"] = MEBIs
end


text:
# von Mises distribution

immutable VonMises <: ContinuousUnivariateDistribution
    μ::Float64  # mean
    κ::Float64  # concentration

    function VonMises(μ::Real, κ::Real)
        κ > zero(κ) || error("kappa must be positive")
        new(float64(μ), float64(κ))
    end

    VonMises(κ::Real) = VonMises(0.0, float64(κ))
    VonMises() = VonMises(0.0, 1.0)
end

## Properties
circmean(d::VonMises) = d.μ
circmedian(d::VonMises) = d.μ
circmode(d::VonMises) = d.μ
circvar(d::VonMises) = 1.0 - besselix(1, d.κ) / besselix(0, d.κ)

function entropy(d::VonMises)
	I0κ = besselix(0.0, d.κ)
	log(twoπ * I0κ) - d.κ * (besselix(1, d.κ) / I0κ - 1.0)
end

## Functions
pdf(d::VonMises, x::Real) = exp(d.κ * (cos(x - d.μ) - 1.0)) / (twoπ * besselix(0, d.κ))
logpdf(d::VonMises, x::Real) = d.κ * (cos(x - d.μ) - 1.0) - log2π - log(besselix(0, d.κ))
cf(d::VonMises, t::Real) = besselix(abs(t), d.k) / besselix(0.0, d.κ) * exp(im * t * d.μ)
cdf(d::VonMises, x::Real) = cdf(d, x, d.μ - π)

function cdf(d::VonMises, x::Real, from::Real)
	const tol = 1.0e-20
	x = mod(x - from, twoπ)
	μ = mod(d.μ - from, twoπ)
	if μ == 0.0
		return vmcdfseries(d.κ, x, tol)
	elseif x <= μ
		upper = mod(x - μ, twoπ)
		if upper == 0.0
			upper = twoπ
		end
		return vmcdfseries(d.κ, upper, tol) - vmcdfseries(d.κ, mod(-μ, twoπ), tol)
	else
		return vmcdfseries(d.κ, x - μ, tol) - vmcdfseries(d.κ, μ, tol)
	end
end

function vmcdfseries(κ::Real, x::Real, tol::Real)
	j, s = 1, 0.0
	while true
		sj = besselix(j, κ) * sin(j * x) / j
		s += sj
		j += 1
		abs(sj) >= tol || break
	end
	x / twoπ + s / (π * besselix(0, κ))
end

## Sampling
sampler(d::VonMises) = VonMisesSampler(d.μ, d.κ)


text:
struct QuerySynthesisOCCOracle <: Oracle
    classifier::SVDD.OCClassifier
end

function QuerySynthesisOCCOracle(classifier_type, init_strategy, data_file::String, solver; classifier_params=Dict{Symbol, Any}(), adjust_K=true)
    data, labels = load_data(data_file)
    return QuerySynthesisOCCOracle(classifier_type, init_strategy, data, labels, solver, classifier_params=classifier_params, adjust_K=adjust_K)
end

function QuerySynthesisOCCOracle(classifier_type, init_strategy, data, labels, solver; classifier_params=Dict{Symbol, Any}(), adjust_K=true)
    pools = convert_labels_to_learning(labels)
    oracle = SVDD.instantiate(classifier_type, data, pools, classifier_params)
    SVDD.initialize!(oracle, init_strategy)
    SVDD.set_adjust_K!(oracle, true)
    SVDD.fit!(oracle, solver)
    return QuerySynthesisOCCOracle(oracle)
end

function QuerySynthesisOCCOracle(data, labels, params::Dict{Symbol, Any})
    # inline function declaration
    check_oracle_parameter(params, key) = key in keys(params) || throw(ArgumentError("Parameter '$key' missing for QuerySynthesisOCCOracle."))
    check_oracle_parameter(params, :classifier_type)
    check_oracle_parameter(params, :init_strategy)
    check_oracle_parameter(params, :solver)
    additional_params = Dict(k => params[k] for k in [:classifier_params, :adjust_K] if k in keys(params))
    return QuerySynthesisOCCOracle(params[:classifier_type], params[:init_strategy], data, labels, params[:solver], additional_params...)
end

function ask_oracle(oracle::QuerySynthesisOCCOracle, query_objects::Array{T, 2})::Vector{Symbol} where T <: Real
    return SVDD.classify.(SVDD.predict(oracle.classifier, query_objects))
end


text:
mutable struct ProblemData{T,X,U,D,O,FX,FU,FW,OX,OU,OXX,OUU,OUX}
    # current trajectory
    states::Vector{X}
    actions::Vector{U}

    # disturbance trajectory
    parameters::Vector{D}

    # nominal trajectory
    nominal_states::Vector{X}
    nominal_actions::Vector{U}

    # model data
    model::ModelData{T,FX,FU,FW}

    # objective data
    objective::ObjectiveData{O,OX,OU,OXX,OUU,OUX}

    # trajectory: z = (x1,..., xT, u1,..., uT-1) | Δz = (Δx1..., ΔxT, Δu1,..., ΔuT-1)
    trajectory::Vector{T}
end

function problem_data(dynamics, costs;
    parameters=[[zeros(d.num_parameter) for d in dynamics]..., zeros(0)])

    length(parameters) == length(dynamics) && (parameters = [parameters..., zeros(0)])
    @assert length(dynamics) + 1 == length(parameters)
    @assert length(dynamics) + 1 == length(costs)

	states = [[zeros(d.num_state) for d in dynamics]...,
            zeros(dynamics[end].num_next_state)]
    actions = [[zeros(d.num_action) for d in dynamics]..., zeros(0)]

    nominal_states = [[zeros(d.num_state) for d in dynamics]...,
            zeros(dynamics[end].num_next_state)]
    nominal_actions = [[zeros(d.num_action) for d in dynamics]..., zeros(0)]

    model = model_data(dynamics)
    objective = objective_data(dynamics, costs)

    trajectory = zeros(num_trajectory(dynamics))

    ProblemData(states, actions, parameters, nominal_states, nominal_actions, model, objective, trajectory)
end


text:
cd(@__DIR__); include("setups/distROT_distTSD_ratio.jl")
gr(dpi = 200)

## (a) P64
include("setups/path64.jl")
ρ1 = generate_ROT_TSD_ratio(nsim, 𝚽, ∇𝚽, 𝛌, Q)
plt = ROT_TSD_ratio_histogram(ρ1)
savefig(plt, "../figs/Path64_ROT_TSD.png")

## (b) P₇ x P₃
include("setups/grid7x3.jl")
ρ2 = generate_ROT_TSD_ratio(nsim, 𝚽, ∇𝚽, 𝛌, Q)
plt = ROT_TSD_ratio_histogram(ρ2)
savefig(plt, "../figs/Grid7x3_ROT_TSD.png")

## (c) Erdos Rényi
include("setups/er.jl")
ρ3 = generate_ROT_TSD_ratio(nsim, 𝚽, ∇𝚽, 𝛌, Q)
plt = ROT_TSD_ratio_histogram(ρ3)
savefig(plt, "../figs/Erdos_Renyi_ROT_TSD.png")

## (d) weighted RGC100
include("setups/rgc100.jl")
ρ4 = generate_ROT_TSD_ratio(nsim, 𝚽, ∇𝚽, 𝛌, Q; edge_length = edge_length)
plt = ROT_TSD_ratio_histogram(ρ4)
savefig(plt, "../figs/wRGC100_ROT_TSD.png")

## Table 4.1
display_basic_stats([ρ1, ρ2, ρ3, ρ4])


text:
#!/usr/bin/env julia
#
println(" Running with ", VERSION)
macro VarName(var)
    string(var)
end

function printsum(varname, a)
    println(" ", varname, " => ", summary(a), ": ", repr(a))
end

function printline(varname, a)
    println(" ", varname, " => ", a)
end

struct Person
	name::AbstractString
	male::Bool
	age::Float64
	children::Int
end

p = Person("Julia", false, 4, 0)
printsum(@VarName(p), p)
#> Person: Person("Julia",false,4.0,0)

people = Person[]
push!(people, Person("Steve", true, 42, 0))
push!(people, Person("Jade", false, 17, 3))
printsum(@VarName(people), people)
#> 2-element Array{Person,1}: [Person("Steve",true,42.0,0),Person("Jade",false,17.0,3)]

# types may also contains arrays and dicts
# constructor functions can be defined to easily create objects
struct Family
	name::AbstractString
	members::Array{AbstractString, 1}
	extended::Bool
	# constructor that takes one argument and generates a default
	# for the other two values
	Family(name::AbstractString) = new(name, AbstractString[], false)
	# constructor that takes two arguements and infers the third
	Family(name::AbstractString, members) = new(name, members, length(members) > 3)
end

fam1 = Family("blogs")
printline(@VarName(fam1), fam1)
#> Family("blogs",AbstractString[],false)
fam2 = Family("jones", ["anna", "bob", "charlie", "dick"])
printline(@VarName(fam2), fam2)
#> Family("jones",AbstractString["anna","bob","charlie","dick"],true)

max_stars_repo_name: wbeebe/julia

max_issues_repo_name: wbeebe/julia

max_forks_repo_name: wbeebe/julia


text:
"""
    greenseltype(::Type{DQMC}, m::Model)

Returns the type of the elements of the Green's function matrix. Defaults to
    `ComplexF64`.
"""
greenseltype(::Type{DQMC}, m::Model) = ComplexF64

"""
    hoppingeltype(::Type{DQMC}, m::Model)

Returns the type of the elements of the hopping matrix. Defaults to `Float64`.
"""
hoppingeltype(::Type{DQMC}, m::Model) = Float64

"""
    interaction_matrix_type(T::Type{DQMC}, m::Model)

Returns the (matrix) type of the interaction matrix. Defaults to
`Matrix{greenseltype(T, m)}`.
"""
interaction_matrix_type(T::Type{DQMC}, m::Model) = Matrix{greenseltype(T, m)}

"""
    hopping_matrix_type(T::Type{DQMC}, m::Model)

Returns the (matrix) type of the hopping matrix. Defaults to
`Matrix{hoppingeltype(T, m)}`.
"""
hopping_matrix_type(T::Type{DQMC}, m::Model) = Matrix{hoppingeltype(T, m)}

"""
    greens_matrix_type(T::Type{DQMC}, m::Model)

Returns the (matrix) type of the greens and most work matrices. Defaults to
`Matrix{greenseltype(T, m)}`.
"""
greens_matrix_type(T::Type{DQMC}, m::Model) = Matrix{greenseltype(T, m)}

"""
    init_interaction_matrix(m::Model)

Initializes the interaction matrix.
"""
function init_interaction_matrix(m::Model)
    N = length(lattice(m))
    flv = nflavors(m)
    zeros(greenseltype(DQMC, m), N*flv, N*flv)
end

"""
    energy_boson(mc::DQMC, m::Model, conf)

Calculate bosonic part (non-Green's function determinant part) of energy for
configuration `conf` for Model `m`.

This is required for global and parallel updates as well as boson energy
measurements.
"""
energy_boson(mc::DQMC, m::Model, conf) = throw(MethodError(energy_boson, (mc, m, conf)))


text:
printstyled("  gcdist\n", color=:light_green)
src = [46.8523, -121.7603]
src2 = [48.7767, -121.8144]
rec = [45.5135 -122.6801; 44.0442 -123.0925; 42.3265 -122.8756]

G0 = gcdist(src, rec)                           # vec, arr
G1 = gcdist(src[1], src[2], rec)                # lat, lon, arr
G2 = gcdist(src[1], src[2], rec[1,1], rec[1,2]) # s_lat, s_lon, r_lat, r_lon
G3 = gcdist(vcat(src',src2'), rec[1,:])             # arr, arr
G4 = gcdist(vcat(src',src2'), rec[1,:])             # arr, rec[1,:]

@test G0 == G1
@test G0[1,:] == G1[1,:] == G2[1,:]
@test G3 == G4


text:
# logreg.jl
#
# Code for logistic regression, but this version uses mini-batches.
#
# In this case, for a triple (a, b, x), where a, x are n-vectors and b
# is {-/+ 1} valued, the objective is
#
#  F(x, a, b) = log(1 + exp(-b * a' * x)).

using Distributions;

# (A, b, x) = GenerateLogisticData(m::Int64 = 200, n::Int64 = 20;
#                                  condition_number::Float64 = 1.0,
#                                  snr::Float64 = 1.0)
#
# Generates logistic regression data where the data matrix A is m-by-n,
# the label (target) vector b is length m, and the distribution of b is
#
#  p(b | a, x) = exp(b * a' * x) / (1 + exp(b * a' * x)),
#
# for b in {-1, 1}. The parameter snr governs the mean magnitude of
# the inner products A * x, that is, the margins. (High snr means less
# noise in the labels).
function GenerateLogisticData(m::Int64 = 200, n::Int64 = 20;
                              condition_number::Float64 = 1.0,
                              snr::Float64 = 1.0)
  A = randn(m, n);
  if (condition_number > 1.0)
    # Make A have condition number given
    condition_vec = Vector(LinRange(1.0, condition_number, n));
    condition_vec = condition_vec / norm(condition_vec);
    A = A .* repeat(condition_vec', outer = (m, 1));
  end
  x = randn(n);
  mean_inner_product = mean(abs.(A * x));
  # Make it so that the mean inner product between data and x is equal
  # to snr (the signal to noise ratio).
  x = snr * x / mean_inner_product;
  margins = A * x;
  probs = 1 ./ (1 .+ exp.(-margins));
  b = sign.(rand(m) .- 1 + probs);
  return (A, b, x);
end

# o = LogisticRegressionObj(A::Matrix{Float64}, b::Vector{Float64},
#                           x::Vector{Float64})
#
# Computes and returns logistic regression objective (mean log loss).
function LogisticRegressionObj(A::Matrix{Float64}, b::Vector{Float64},
                               x::Vector{Float64})
  margins = b .* (A * x);
  pos_inds = (margins .> 0);
  neg_inds = (margins .<= 0);
  obj = (sum(log.(1 .+ exp.(-margins[pos_inds])))
         + sum(log.(1 .+ exp.(margins[neg_inds])) - margins[neg_inds]));
  return obj / size(A, 1);
end

# g = LogisticRegressionGrad(A::Matrix{Float64}, b::Vector{Float64},
#                            x::Vector{Float64})
#
# Computes and returns the gradient of the logistic regression loss.
function LogisticRegressionGrad(A::Matrix{Float64}, b::Vector{Float64},
                                x::Vector{Float64})
  margins = b .* (A * x);
  probs = b ./ (1 .+ exp.(margins));
  grad = -A' * probs;
  return grad / size(A, 1);
end

# x = SGUpdateLogistic(A::SubArray{Float64, 2}, b::SubArray{Float64, 1},
#                        x_init::Vector{Float64}, stepsize::Float64)
#
# Computes and applies a minibatch stochastic gradient update to the
# vector x_init for the logistic regression loss. The entire minibatch is
# in the SubArrays (A, b), which can be treated like matrices.
function SGUpdateLogistic(A::SubArray{Float64, 2}, b::SubArray{Float64, 1},
                          x_init::Vector{Float64}, stepsize::Float64)
  # Replace this code to return the correct update.
  return x_init;
end

# x = TruncatedUpdateLogistic(A::SubArray{Float64, 2},
#                             b::SubArray{Float64, 1},
#                             x_init::Vector{Float64}, stepsize::Float64)
#
# Performs a mini-batch update using the truncated model for logistic
# regression. In particular, letting a_i be the rows of the matrix A and
# b_i the associated labels, uses that the logistic loss has lower-bound 0
# to make the approximation
#
#  F_approx(x) = max{ F_0 + g' * (x - x_init), 0 }
#
# where F_0 = (1/m) sum_{i=1}^m log(1 + exp(-b_i * a_i' * x)) and
#
#   g = (d / dx) (1/m) sum_{i = 1}^m log(1 + exp(-b_i * a_i' * x))
#
# is the gradient on the mini-batch.
function TruncatedUpdateLogistic(A::SubArray{Float64, 2},
                                 b::SubArray{Float64, 1},
                                 x_init::Vector{Float64}, stepsize::Float64)
  # Replace this code to return the correct update.
  return x_init;
end

function TruncatedUpdateLogistic(a::Vector{Float64},
                                 b::Float64,
                                 x_init::Vector{Float64}, stepsize::Float64)
  # Replace this code to return the correct update.
  return x_init;
end

function SGUpdateLogistic(a::Vector{Float64}, b::Float64,
                          x_init::Vector{Float64}, stepsize::Float64)
  # Replace this code to return the correct update.
  return x_init;
end

max_stars_repo_path: Machine Learning Summer School 2019 (London, UK)/tutorials/optimization/logreg.jl


text:
"""
    SegmentedPNorm{V <: AbstractVector{<:AbstractFloat}} <: AbstractAggregation

[`AbstractAggregation`](@ref) implementing segmented p-norm aggregation:

``
f(\\{x_1, \\ldots, x_k\\}; p, c) = \\left(\\frac{1}{k} \\sum_{i = 1}^{k} \\vert x_i - c \\vert ^ {p} \\right)^{\\frac{1}{p}}
``

Stores a vector of parameters `ψ` that are filled into the resulting matrix in case an empty bag is encountered,
and vectors of parameters `p` and `c` used during computation.

See also: [`AbstractAggregation`](@ref), [`AggregationStack`](@ref), [`pnorm_aggregation`](@ref),
    [`SegmentedMax`](@ref), [`SegmentedMean`](@ref), [`SegmentedSum`](@ref), [`SegmentedLSE`](@ref).
"""
struct SegmentedPNorm{V <: AbstractVector{<:AbstractFloat}} <: AbstractAggregation
    ψ::V
    ρ::V
    c::V
end

Flux.@functor SegmentedPNorm

SegmentedPNorm(T::Type, d::Int) = SegmentedPNorm(zeros(T, d), randn(T, d), randn(T, d))
SegmentedPNorm(d::Int) = SegmentedPNorm(Float32, d)

Flux.@forward SegmentedPNorm.ψ Base.getindex, Base.length, Base.size, Base.firstindex, Base.lastindex,
        Base.first, Base.last, Base.iterate, Base.eltype

Base.vcat(as::SegmentedPNorm...) = reduce(vcat, as |> collect)
function Base.reduce(::typeof(vcat), as::Vector{<:SegmentedPNorm})
    SegmentedPNorm(reduce(vcat, [a.ψ for a in as]),
                   reduce(vcat, [a.ρ for a in as]),
                   reduce(vcat, [a.c for a in as]))
end

p_map(ρ::T) where T = one(T) + softplus(ρ)
p_map(ρ::AbstractArray) = p_map.(ρ)
inv_p_map(p::T) where T = relu(p - one(T)) + log1p(-exp(-abs(p - one(T))))
inv_p_map(ρ::AbstractArray) = inv_p_map.(ρ)

function (a::SegmentedPNorm)(x::Missing, bags::AbstractBags,
                                w::Optional{AbstractVecOrMat}=nothing)
    _check_agg(a, x)
    segmented_pnorm_forw(x, a.ψ, nothing, bags, w)
end
function (a::SegmentedPNorm)(x::AbstractMatrix{T}, bags::AbstractBags,
                                w::Optional{AbstractVecOrMat{T}}=nothing) where T
    _check_agg(a, x)
    segmented_pnorm_forw(x .- a.c, a.ψ, p_map(a.ρ), bags, w)
end

function _pnorm_precomp(x::AbstractMatrix, bags)
    M = ones(eltype(x), size(x, 1), length(bags))
    @inbounds for (bi, b) in enumerate(bags)
        if !isempty(b)
            for j in b
                for i in 1:size(x, 1)
                    M[i, bi] = max(M[i, bi], abs(x[i, j]))
                end
            end
        end
    end
    M
end

function _segmented_pnorm_norm(a::AbstractMatrix, ψ::AbstractVector, p::AbstractVector,
                               bags::AbstractBags, w, M)
    isnothing(w) || @assert all(w .> 0)
    t = promote_type(eltype(a), eltype(ψ))
    y = zeros(t, size(a, 1), length(bags))
    M = _pnorm_precomp(a, bags)
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(ψ)
                y[i, bi] = ψ[i]
            end
        else
            ws = _bagnorm(w, b)
            for j in b
                for i in 1:size(a, 1)
                    y[i, bi] += _weight(w, i, j, t) * abs(a[i, j]/M[i, bi])^p[i]
                end
            end
            for i in 1:size(a, 1)
                y[i, bi] = M[i, bi] * (y[i, bi] / _weightsum(ws, i)) ^ (one(t)/p[i])
            end
        end
    end
    y
end

segmented_pnorm_forw(::Missing, ψ::AbstractVector, p, bags::AbstractBags, w) = repeat(ψ, 1, length(bags))
function segmented_pnorm_forw(a::Maybe{AbstractMatrix}, ψ::AbstractVector, p::AbstractVector, bags::AbstractBags, w)
    M = _pnorm_precomp(a, bags)
   _segmented_pnorm_norm(a, ψ, p, bags, w, M)
end

function segmented_pnorm_back(Δ, y, a, ψ, p, bags, w, M)
    da = similar(a)
    dp = zero(p)
    dps1 = zero(p)
    dps2 = zero(p)
    dψ = zero(ψ)
    dw = isnothing(w) ? nothing : zero(w)
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(ψ)
                dψ[i] += Δ[i, bi]
            end
        else
            ws = _bagnorm(w, b)
            dps1 .= zero(eltype(p))
            dps2 .= zero(eltype(p))
            for j in b
                for i in 1:size(a, 1)
                    ab = abs(a[i, j])
                    da[i, j] = Δ[i, bi] * _weight(w, i, j, eltype(p)) * sign(a[i, j]) / _weightsum(ws, i)
                    da[i, j] *= (ab / y[i, bi]) ^ (p[i] - one(eltype(p)))
                    ww = _weight(w, i, j, eltype(p)) * (ab / M[i, bi]) ^ p[i]
                    dps1[i] += ww * log(ab)
                    dps2[i] += ww
                end
            end
            for i in 1:size(a, 1)
                t = y[i, bi] / p[i]
                t *= dps1[i] / dps2[i] - (p[i] * log(M[i, bi]) + log(dps2[i]) - log(_weightsum(ws, i))) / p[i]
                dp[i] += Δ[i, bi] * t
            end
        end
    end
    da, dψ, dp, NoTangent(), @not_implemented("Not implemented yet!")
end

function segmented_pnorm_back(Δ, y, ψ, bags)
    dψ = zero(ψ)
    @inbounds for (bi, b) in enumerate(bags)
        for i in eachindex(ψ)
            dψ[i] += Δ[i, bi]
        end
    end
    ZeroTangent(), dψ, ZeroTangent(), NoTangent(), @not_implemented("Not implemented yet!")
end

∇dw_segmented_pnorm!(dw::ZeroTangent, Δ, a, y, w::Nothing, i, j, bi) = error("Not implemented yet!")
∇dw_segmented_pnorm!(dw::AbstractVector, Δ, a, y, w::AbstractVector, i, j, bi) = error("Not implemented yet!")
∇dw_segmented_pnorm!(dw::AbstractMatrix, Δ, a, y, w::AbstractMatrix, i, j, bi) = error("Not implemented yet!")

function ChainRulesCore.rrule(::typeof(segmented_pnorm_forw), a::AbstractMatrix, ψ, p, bags, w)
    M = _pnorm_precomp(a, bags)
    y = _segmented_pnorm_norm(a, ψ, p, bags, w, M)
    grad = Δ -> (NoTangent(), segmented_pnorm_back(Δ, y, a, ψ, p, bags, w, M)...)
    y, grad
end

function ChainRulesCore.rrule(::typeof(segmented_pnorm_forw), a::Missing, ψ, p, bags, w)
    y = segmented_pnorm_forw(a, ψ, p, bags, w)
    grad = Δ -> (NoTangent(), segmented_pnorm_back(Δ, y, ψ, bags)...)
    y, grad
end


text:
using Benchmarkplotting
using Gadfly
using Statistics

cases = [
    :range100   => 1:100,
    :range10000 => 1:10000,
    :array100   => [(1:100)...],
    :array10000 => [(1:10000)...],
]

implementations = [
    :map       =>
        function (x)
            map(e -> e + 1, x)
        end,
    :list_comp =>
        function (x)
            [e + 1 for e in x]
        end,
    :handwritten_loop =>
        function(x)
            i = 1
            arr = []
            while (iter = iterate(x, i); iter !== nothing)
                e, i = iter
                push!(arr, e + 1)
            end
            arr
        end,
    :handwritten_loop_prealloc =>
        function(x)
            i = 1
            arr = fill(0, length(x))
            while (iter = iterate(x, i); iter !== nothing)
                e, i = iter
                push!(arr, e + 1)
            end
            arr
    end,
    :for_iter =>
        function(x)
            arr = []
            for e in x
                push!(arr, e + 1)
            end
            arr
        end
]

theme = Theme(
    guide_title_position = :left,
    colorkey_swatch_shape = :circle,
    minor_label_font = "Consolas",
    major_label_font = "Consolas",
    point_size=6px
)

criterion(benchmark_result) = (time = mean(benchmark_result.times), )

df = bcompare(criterion, cases, implementations, quiet=false)
res = report(:time, df, Scale.y_log10, theme)
@info res[2]
draw(SVG("example.svg", 9inch, 6inch), res[1])


text:
skipcache(::Any) = false
skipcache(::Type{<:ConstantBlockHeader}) = true
skipcache(::Type{<:RunInfoBlockHeader}) = true
skipcache(::Type{<:CPUSplitBlockHeader}) = true
skipcache(::Type{<:PlainMeshBlockHeader}) = true

function cached_read(f, block, cache)
    skipcache(typeof(block)) && return read(f, block)

    raw_data = get!(cache, (f.name, block)) do
        read!(f, block)
    end

    apply_normalization!(raw_data, block)

    return add_units(raw_data, block)
end