text:
"""
    iscoprime(a,b)

Check if the pair of polynomials `a` and `b` is coprime, that is, it has no common factors.

# Examples

```julia
julia> a = Polynomial([1,2,3],:s);
julia> b = Polynomial([1,2],:s);
julia> iscoprime(a,b)
true

julia> a = Polynomial([1,2,3],:s);
julia> b = a*Polynomial([1,2],:s);
julia> iscoprime(a,b)
false
```
"""
function iscoprime(a::Polynomial,b::Polynomial)
    da = degree(a)
    db = degree(b)
    dx = db-1
    S = sylvestermatrix(a,b)
    return rank(S)==(da+db)
end

"""
    isschurstable(a)

Check if the polynomial `a` is Schur stable, that is, it has all its roots inside the unit disk.

# Examples

```julia
julia> a = fromroots([-1/2,0, 2/3,-1/3im, 1/3im]);
julia> isschurstable(a)
true

julia> a = fromroots([2,0, 2/3,-3im, 3im]);
julia> !isschurstable(a)
false
```
"""
function isschurstable(a::Polynomial)
    r = roots(a)
    all(abs.(r) .<= 1)
end

"""
    ishurwitzstable(a::Polynomial)

Check if the `a` is Schur stable, that is, it has all its roots inside the left complex half-plane.

# Examples

```julia
julia> a = fromroots([-1, -2, -3+4im, -3-4im]);
julia> ishurwitzstable(a)
true

julia> a = fromroots([-1, -2, 3+4im, 3-4im]);
julia> !ishurwitzstable(a)
false
```
"""
function ishurwitzstable(a::Polynomial)
    r = roots(a)
    return all(real(r) .<= 0)
end

"""
    isconjsymmetric(a::Polynomial)

Check if the scalar univariate polynomial `a` is conjugate symmetric, that is, if `a=ã`.

For a scalar univariate polynomial `a` with real coefficients, check if `a(s)=a(-s)`. If the coefficients are complex, the condition then extends to `a(s)=ā(-s)`, where `ā` is obtained from `a` by complex-conjugating the coefficients. A necessary and sufficient condition for conjugate symmetry is that the coefficients with odd powers of the variable are zero and the coefficients with even powers are real.

# Examples
```julia
julia> a = Polynomial([1,0,2,0,3],:s)
Polynomial(1 + 2*s^2 + 3*s^4)

julia> isconjsymmetric(a)
true
```
"""
function isconjsymmetric(a::Polynomial)
    return isequal(a,cconj(a))
end

"""
    isconjsymmetric(a::LaurentPolynomial)

Check if the scalar univariate Laurent polynomial `a` is conjugate symmetric, that is, if `a=ã`.

For a scalar univariate Laurent polynomial `a` with real coefficients, check if `a(z)=a(1/z)`. If the coefficients are complex, the condition then extends to `a(z)=ā(1/z)`, where `ā` is obtained from `a` by complex-conjugating the coefficients. A necessary and sufficient condition is that the coefficients are real and their order in the coefficient vector can be flipped with no impact.

# Examples
```julia
julia> a = LaurentPolynomial([3,2,1,2,3],-2:2,:z)
LaurentPolynomial(3*z⁻² + 2*z⁻¹ + 1 + 2*z + 3*z²)

julia> isconjsymmetric(a)
true
```
"""
function isconjsymmetric(a::LaurentPolynomial)
    return isequal(a,dconj(a))
end

ext: jl


text:
using Documenter, Julia_docs_test

makedocs(modules=[Julia_docs_test],
        doctest=true)

deploydocs(deps   = Deps.pip("mkdocs", "python-markdown-math"),
          repo = "github.com/nbavafa/Julia_docs_test.git",
          julia = "0.6",
	  osname = "osx
       )


text:
function my_lin_interp(grid,vals)
    function func(x::Real)
        lower_index = searchsortedlast(grid,x)
        if lower_index == 0
            interpolated_value  = ((vals[2]　-　vals[1])　/　(grid[2]　-　grid[1]))　*　(x　-　grid[1])　+　vals[1]
        elseif lower_index == length(grid)
            interpolated_value = ((vals[end]　-　vals[end　-　1])　/　(grid[end]　-　grid[end　-　1]))　*　(x　-　grid[end　-　1])　+　vals[end　-　1]
        else
            interpolated_value =  (vals[lower_index+1] - vals[lower_index]) * (x - grid[lower_index]) /
            (grid[lower_index+1] - grid[lower_index])  + vals[lower_index]
        end
        return interpolated_value
    end
    function func{T<:Real}(x::AbstractVector{T})
        n = length(x)
        out = Array(Float64, n)
        for i in 1 : n
            out[i] = func(x[i])
        end
        return out
    end
    return func
end


text:
module VimCompletion

export serverstart

include("complete.jl")
include("server.jl")

end # module


text:
using Rimu
using Rimu: fciqmc_step!
using Test

@testset "Allocations" begin
    # The purpose of these tests is to find type instabilities that might appear as the
    # Julia compiler changes. If allocations suddenly increase by a lot, there is a good
    # chance that dynamic dispatch is happening somewhere in the code.
    b1 = near_uniform(BoseFS{10,10})
    b2 = near_uniform(BoseFS{50,50})
    b3 = near_uniform(BoseFS{100,100})

    f1 = near_uniform(FermiFS{9,10})
    f2 = near_uniform(FermiFS{24,50})
    f3 = near_uniform(FermiFS{49,100})

    for H in (
        HubbardReal1D(b1),
        HubbardReal1D(b2),
        HubbardReal1D(b3),

        HubbardMom1D(b1),
        HubbardMom1D(b2),
        HubbardMom1D(b3),

        ExtendedHubbardReal1D(b1),
        ExtendedHubbardReal1D(b2),
        ExtendedHubbardReal1D(b3),

        BoseHubbardReal1D2C(BoseFS2C(b1, b1)),
        BoseHubbardReal1D2C(BoseFS2C(b2, b2)),
        BoseHubbardReal1D2C(BoseFS2C(b3, b3)),

        BoseHubbardMom1D2C(BoseFS2C(b1, b1)),
        BoseHubbardMom1D2C(BoseFS2C(b2, b2)),
        BoseHubbardMom1D2C(BoseFS2C(b3, b3)),

        HubbardRealSpace(b1),
        HubbardRealSpace(b2),
        HubbardRealSpace(b3),

        HubbardRealSpace(f1),
        HubbardRealSpace(f2),
        HubbardRealSpace(f3),

        HubbardRealSpace(CompositeFS(f1, f1)),
        HubbardRealSpace(CompositeFS(f2, f2)),
        HubbardRealSpace(CompositeFS(f3, f3)),
    )
        add = starting_address(H)
        for dv_type in (DVec, InitiatorDVec)
            for style in (
                IsDeterministic(),
                IsStochasticInteger(),
                IsStochasticWithThreshold(),
                IsDynamicSemistochastic(),
            )
                hamname = string(
                    nameof(typeof(H)), "(", num_modes(add), ")/", nameof(typeof(style))
                )
                @testset "Allocations for $(hamname)" begin
                    dτ = if num_modes(add) == 10
                        1e-4
                    elseif num_modes(add) == 50
                        1e-4
                    else
                        1e-6
                    end

                    dv = dv_type(add => 1.0, style=IsDynamicSemistochastic())
                    sizehint!(dv, 500_000)

                    # Warmup for lomc!
                   params = RunTillLastStep(shift=float(diagonal_element(H, add)), dτ=dτ)
                    _, st = lomc!(
                        H, dv; params, threading=false, maxlength=10_000, laststep=1
                    )

                    r = only(st.replicas)
                    p = r.params
                    t = Rimu.NoThreading()

                    # Warmup for step!
                    fciqmc_step!(t, r.w, H, r.v, p.shift, dτ)
                    fciqmc_step!(t, r.w, H, r.v, p.shift, dτ)
                    fciqmc_step!(t, r.w, H, r.v, p.shift, dτ)
                    fciqmc_step!(t, r.w, H, r.v, p.shift, dτ)

                    allocs_step = @allocated fciqmc_step!(t, r.w, H, r.v, p.shift, dτ)
                    @test allocs_step ≤ 512

                    dv = dv_type(add => 1.0, style=IsDynamicSemistochastic())
                    allocs_full = @allocated lomc!(
                        H, dv; dτ, laststep=200, threading=false, maxlength=10_000
                    )
                    @test allocs_full ≤ 1e8 # 100MiB

                    # Print out the results to make it easier to find problems.
                    print(rpad(hamname, 50))
                    print(": per step ", allocs_step, ", full ", allocs_full/(1024^2), "M\n")
                end
            end
        end
    end
end


text:
module Tune

using CpuId
using Hwloc
using BenchmarkTools
using MemoizedMethods
using DiskCaches
using Scratch
using SHA

@memoize arch_id() = (sha2_256(string(cpuinfo()))...,)

export @memoize DiskCache @get_scratch! arch_id

end


text:
using BenchmarkTools

function answer1(input = "input0")
    numtotal = 0
    for line in readlines(input)
        line = split.(line,"")
        last_opened = 5
        opened = []
        for character in line
            if character == "("
                last_opened = 1
                append!(opened,1)
            elseif character == "["
                last_opened = 2
                append!(opened,2)
            elseif character == "{"
                last_opened = 3
                append!(opened,3)
            elseif character == "<"
                last_opened = 4
                append!(opened,4)
            elseif character == ")"
                if opened[end] != 1
                    numtotal +=3
                    break
                else
                    pop!(opened)
                end
            elseif character == "]"
                if opened[end] != 2
                    numtotal +=57
                    break
                else
                    pop!(opened)
                end
            elseif character == "}"
                if opened[end] != 3
                    numtotal +=1197
                    break
                else
                    pop!(opened)
                end
            elseif character == ">"
                if opened[end] != 4
                    numtotal +=25137
                    break
                else
                    pop!(opened)
                end
            end
        end
    end
    numtotal
end
@show answer1("input")

using Statistics

function answer2(input = "input0")
    numtotal = 0
    autocomplete_score = []
    for line in readlines(input)
        line = split.(line,"")
        last_opened = 5
        opened = []
        line_autocomplete_score = 0
        corrupted = false
        for character in line
            if character == "("
                last_opened = 1
                append!(opened,1)
            elseif character == "["
                last_opened = 2
                append!(opened,2)
            elseif character == "{"
                last_opened = 3
                append!(opened,3)
            elseif character == "<"
                last_opened = 4
                append!(opened,4)
            elseif character == ")"
                if opened[end] != 1
                    numtotal +=3
                    corrupted = true
                    break
                else
                    pop!(opened)
                end
            elseif character == "]"
                if opened[end] != 2
                    numtotal +=57
                    corrupted = true
                    break
                else
                    pop!(opened)
                end
            elseif character == "}"
                if opened[end] != 3
                    numtotal +=1197
                    corrupted = true
                    break
                else
                    pop!(opened)
                end
            elseif character == ">"
                if opened[end] != 4
                    numtotal +=25137
                    corrupted = true
                    break
                else
                    pop!(opened)
                end
            end
        end
        if !corrupted
            for missing_character in reverse(opened)
                # @show missing_character
                line_autocomplete_score = line_autocomplete_score*5 + missing_character
            end
        append!(autocomplete_score, line_autocomplete_score)
        end
    end
    Int(median(autocomplete_score))
end

@show answer2("input")


text:
using OrdinaryDiffEq, Distributions,
      DiffEqUncertainty, Test, Quadrature, Cubature, Cuba,
      FiniteDiff, Zygote, ForwardDiff, DiffEqGPU, DiffEqSensitivity

quadalgs = [HCubatureJL(), CubatureJLh(), CubatureJLp(), CubaSUAVE(), CubaDivonne(), CubaCuhre()]
batchmode = [EnsembleSerial(), EnsembleThreads(), EnsembleCPUArray()]#, EnsembleGPUArray()]

function eom!(du,u,p,t)
  @inbounds begin
    du[1] = u[2]
    du[2] = -p[1]*u[1] - p[2]*u[2]
  end
  nothing
end

u0 = [1.0;1.0]
tspan = (0.0,3.0)
p = [1.0, 2.0]
prob = ODEProblem(eom!,u0,tspan,p)

A = [0.0 1.0; -p[1] -p[2]]
u0s_dist = [Uniform(1,10), Uniform(2,6)]

@testset "Koopman Expectation, nout = 1" begin
  g(sol) = sol[1,end]
  analytical = (exp(A*tspan[end])*mean.(u0s_dist))[1]

  for alg ∈ quadalgs
    @info "$alg"
    @test expectation(g, prob, u0s_dist, p, Koopman(), Tsit5(); quadalg = alg)[1] ≈ analytical rtol=1e-2
  end
end

@testset "Koopman Expectation, nout > 1" begin
  g(sol) = sol[:,end]
  analytical = (exp(A*tspan[end])*mean.(u0s_dist))

  for alg ∈ quadalgs
    @info "$alg"
    @test expectation(g, prob, u0s_dist, p, Koopman(), Tsit5(); quadalg = alg, nout =2) ≈ analytical rtol=1e-2
  end
end

@testset "Koopman Expectation, nout = 1, batch" begin
  g(sol) = sol[1,end]
  analytical = (exp(A*tspan[end])*mean.(u0s_dist))[1]

  for bmode ∈ batchmode
    for alg ∈ quadalgs
      if alg isa HCubatureJL
        continue
      end
      @info "nout = 1, batch mode = $bmode, $alg"
      @test expectation(g, prob, u0s_dist, p, Koopman(), Tsit5(), bmode; quadalg = alg, batch = 15)[1] ≈ analytical rtol=1e-2
    end
  end
end

@testset "Koopman Expectation, nout > 1, batch" begin
  g(sol) = sol[:,end]
  analytical = (exp(A*tspan[end])*mean.(u0s_dist))

  for bmode ∈ batchmode
    for alg ∈ quadalgs
      if alg isa HCubatureJL
        continue
      end
      @info "nout = 2, batch mode = $bmode, $alg"
      res = expectation(g, prob, u0s_dist, p, Koopman(), Tsit5(), bmode;
                        quadalg = alg, batch = 15, nout=2)
      @test res ≈ analytical rtol=1e-2
    end
  end
end

############## Koopman AD ###############

@testset "Koopman Expectation AD" begin
  g(sol) = sol[1,end]
  loss(p, alg) = expectation(g, prob, u0s_dist, p, Koopman(), Tsit5(); quadalg = alg)[1]
  dp1 = FiniteDiff.finite_difference_gradient(p->loss(p, HCubatureJL()),p)
  for alg ∈ quadalgs
    @info "$alg, ForwardDiff"
    @test ForwardDiff.gradient(p->loss(p,alg),p) ≈ dp1 rtol=1e-2
    @info "$alg, Zygote"
    @test Zygote.gradient(p->loss(p,alg),p)[1] ≈ dp1 rtol=1e-2
  end
end

@testset "Koopman Expectation AD, batch" begin
  g(sol) = sol[1,end]
  loss(p, alg, bmode) = expectation(g, prob, u0s_dist, p, Koopman(), Tsit5(), bmode; quadalg = alg, batch = 10)[1]
  dp1 = FiniteDiff.finite_difference_gradient(p->loss(p, CubatureJLh(), EnsembleThreads()),p)
  for bmode ∈ batchmode
    for alg ∈ quadalgs
      if alg isa HCubatureJL #no batch support
        continue
      end
      @info "$bmode, $alg, ForwardDiff"
      @test_broken ForwardDiff.gradient(p->loss(p,alg,bmode),p) ≈ dp1 rtol=1e-2
      @info "$bmode, $alg, Zygote"
      @test_broken  Zygote.gradient(p->loss(p,alg,bmode),p)[1] ≈ dp1 rtol=1e-2
    end
  end
end

########## Central Moment ##############

function eom!(du,u,p,t)
  @inbounds begin
    du[1] = p[1]*u[1]
  end
  nothing
end

u0 = [1.0]
tspan = (0.0,3.0)
p = [-.3]
prob = ODEProblem(eom!,u0,tspan,p)

u0s_dist = [Uniform(1,10)]

@testset "Koopman Central Moment" begin
  g(sol) = sol[1,end]
  analytical = [0.0, exp(2*p[1]*tspan[end])*var(u0s_dist[1]), 0.0]

  for alg ∈ quadalgs
    if alg isa CubaDivonne || alg isa CubaCuhre  #requires 2D spatial integration
      continue
    end
    r = centralmoment(3, g, prob, u0s_dist, p, Koopman(), Tsit5();
                      ireltol = 1e-8, iabstol = 1e-8, quadalg = alg)
    if alg isa CubaSUAVE
      @test_broken r ≈ analytical rtol=1e-2
    else
      @test r ≈ analytical rtol=1e-2
    end
  end
end

@testset "Koopman Central Moment, batch" begin
  g(sol) = sol[1,end]
  analytical = [0.0, exp(2*p[1]*tspan[end])*var(u0s_dist[1]), 0.0]

  for bmode ∈ batchmode
    for alg ∈ quadalgs
      if alg isa CubaDivonne || alg isa CubaCuhre || alg isa HCubatureJL #requires 2D spatial integration
        continue
      end
      @info "batch mode = $bmode, alg = $alg"
      r = centralmoment(3, g, prob, u0s_dist, p, Koopman(), Tsit5(), bmode;
                        ireltol = 1e-8, iabstol = 1e-8, quadalg = alg, batch=15)
      if alg isa CubaSUAVE
        @test_broken r ≈ analytical rtol=1e-2
      else
        @test r ≈ analytical rtol=1e-2
      end
    end
  end
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/relax.jl
# Defines routines used construct the relaxed subproblem.
#############################################################################

"""
$(FUNCTIONNAME)

Applies the safe cut checks detailed in Khajavirad, 2018 [Khajavirad, Aida,
and Nikolaos V. Sahinidis. "A hybrid LP/NLP paradigm for global optimization
relaxations." Mathematical Programming Computation 10.3 (2018): 383-421] to
ensure that only numerically safe affine relaxations are added. Checks that
i) ``|b| <= safe b`, ii) `safe_l <= abs(ai) <= safe u`, and iii) violates
`safe_l <= abs(ai/aj) <= safe_u`.
"""
function is_safe_cut!(m::GlobalOptimizer{N,T,S}, f::SAF{T}) where {N,T<:AbstractFloat,S}

    # violates |b| <= safe_b
    (abs(f.constant) > m.cut_safe_b) && return false

    safe_l = m.cut_safe_l
    safe_u = m.cut_safe_u
    term_count = length(f.terms)
    for i = 1:term_count
        ai = (@inbounds f.terms[i]).coefficient
        if !iszero(ai)
            ai_abs = abs(ai)
            # violates safe_l <= abs(ai) <= safe_u
            !(safe_l <= abs(ai) <= safe_u) && return false
            for j = i:term_count
                aj = (@inbounds f.terms[j]).coefficient
                if !iszero(aj)
                    # violates safe_l <= abs(ai/aj) <= safe_u
                    !(safe_l <= abs(ai/aj) <= safe_u) && return false
                end
            end
        end
    end

    return true
end

"""
$(FUNCTIONNAME)

Relaxs the constraint by adding an affine constraint to the model.
"""
function relax! end

"""
$(FUNCTIONNAME)

Default routine for relaxing quadratic constraint `func < 0.0` on node `n`.
Takes affine bounds of convex part at point `x0` and secant line bounds on
concave parts.
"""
function affine_relax_quadratic!(func::SQF{T}, d::Dict{Int,T}, saf::SAF{T},
                                 n::NodeBB, sol_to_branch_map::Vector{Int},
                                 x::Vector{T}) where T

    quadratic_constant = func.constant

    # Affine terms only contribute coefficients, so the respective
    # values do not contribute to the cut. Since all quadratic terms
    # are considered to be branch variables we exclude any potential
    # need to retrieve variable bounds from locations other than
    # the node.
    for term in func.quadratic_terms

        a = term.coefficient
        idx1 = term.variable_index_1.value
        idx2 = term.variable_index_2.value
        x1 = x[sol_idx1]                    # TODO FIX THIS
        xL1 = _lower_bound(FullVar, m, idx1)
        xU1 = _upper_bound(FullVar, m, idx1)

        if idx1 == idx2
            if a > zero(T)
                d[idx1] += a*x1
                quadratic_constant -= 0.5*a*x1*x1
            else
                if !isinf(xL1) && !isinf(xU1)
                    d[idx1] += 0.5*a*(xL1 + xU1)
                    quadratic_constant -= 0.5*a*xL1*xU1
                else
                    return false
                end
            end
        else
            x2 = x[sol_idx2]
            xL2 = _lower_bound(FullVar, m, idx2)
            xU2 = _upper_bound(FullVar, m, idx2)
            if a > 0.0
                if (!isinf(xL1) && !isinf(xL2)) &&
                   ((xU1 - xL1)*x2 + (xU2 - xL2)*x1 <= xU1*xU2 - xL1*xL2)
                    d[idx1] += a*xL2
                    d[idx2] += a*xL1
                    quadratic_constant -= a*xL1*xL2
                elseif !isinf(xU1) && !isinf(xU2)
                    d[idx1] += a*xU2
                    d[idx2] += a*xU1
                    quadratic_constant -= a*xU1*xU2
                else
                    return false
                end
            else
                if (!isinf(xU1) && !isinf(xL2)) &&
                   ((xU1 - xL1)*x2 - (xU2 - xL2)*x1 <= xU1*xL2 - xL1*xU2)
                    d[idx1] += a*xL2
                    d[idx2] += a*xU1
                    quadratic_constant -= a*xU1*xL2
                elseif !isinf(xL1) && !isinf(xU2)
                    d[idx1] += a*xU2
                    d[idx2] += a*xL1
                    quadratic_constant -= a*xL1*xU2
                else
                    return false
                end
            end
        end
    end
    foreach(x -> (d[x.variable_index.value] = x.coefficient;), func.affine_terms)

    i = 1
    for (k, v) in d
        saf.terms[i] = SAT(v, VI(k))
        d[k] = 0.0
        i += 1
    end
    saf.constant = quadratic_constant

    return true
end

"""
$(TYPEDSIGNATURES)
"""
function relax!(m::GlobalOptimizer{N,T,S}, f::BufferedQuadraticIneq{T}, safe::Bool) where {N,T<:AbstractFloat,S<:ExtensionType}

    finite_cut_generated = affine_relax_quadratic!(f.func, f.buffer, f.saf, m._current_node, m._sol_to_branch_map, m._current_xref)
    if finite_cut_generated && (!safe || is_safe_cut!(m, f.saf))
        lt = LT(-f.saf.constant + m.abs_constraint_feas_toleranc)
        f.saf.constant = zero(T)
        ci = MOI.add_constraint(m.relaxed_optimizer, f.saf, lt)
        push!(m._buffered_quadratic_ineq_ci, ci)
    end

    return nothing
end

"""
$(TYPEDSIGNATURES)
"""
function relax!(m::GlobalOptimizer{N,T,S}, f::BufferedQuadraticEq{T}, safe::Bool) where {N,T<:AbstractFloat,S<:ExtensionType}

    finite_cut_generated = affine_relax_quadratic!(f.func, f.buffer, f.saf, m._current_node, m._sol_to_branch_map, m._current_xref)
    if finite_cut_generated && (!safe || is_safe_cut!(m, f.saf))
        lt = LT(-f.saf.constant + m.abs_constraint_feas_tolerance)
        f.saf.constant = zero(T)
        ci = MOI.add_constraint(m.relaxed_optimizer, f.saf, lt)
        push!(m._buffered_quadratic_eq_ci, ci)
    end

    finite_cut_generated = affine_relax_quadratic!(f.minus_func, f.buffer, f.saf, m._current_node, m._sol_to_branch_map, m._current_xref)
    if finite_cut_generated && (!safe || is_safe_cut!(m, f.saf))
        lt = LT(-f.saf.constant + m.abs_constraint_feas_tolerance)
        f.saf.constant = zero(T)
        ci = MOI.add_constraint(m.relaxed_optimizer, f.saf, lt)
        push!(m._buffered_quadratic_eq_ci, ci)
    end

    return nothing
end

"""
$(FUNCTIONNAME)
"""
function affine_relax_nonlinear!(f::BufferedNonlinearFunction{MC{N,T}}, evaluator::Evaluator,
                                 use_cvx::Bool, new_pass::Bool) where {N,T<:RelaxTag}

    new_pass && forward_pass!(evaluator, f)
    x = evaluator.x
    finite_cut = true

    if f.expr.isnumber[1]
        f.saf.constant = f.expr.numberstorage[1]
        map!(i -> SAT(0.0, VI(i)), f.saf.terms, f.expr.grad_sparsity)
    else
        setvalue = f.expr.setstorage[1]
        finite_cut &= !(isempty(setvalue) || isnan(setvalue))
        if finite_cut
            value = f.expr.setstorage[1]
            f.saf.constant = use_cvx ? value.cv : -value.cc
            for i = 1:N
                vi = f.expr.grad_sparsity[i]
                c = use_cvx ? value.cv_grad[i] : -value.cc_grad[i]
                f.saf.terms[i] = SAT(c, VI(vi))
                f.saf.constant = sub_round(f.saf.constant , mul_round(coef, x[vi], RoundUp), RoundDown)
            end
            bnd_used =  use_cvx ? -f.upper_bound : f.lower_bound
            f.saf.constant = add_round(f.saf.constant, bnd_used, RoundDown)
        end
    end
    return finite_cut
end

"""
$(TYPEDSIGNATURES)
"""
function check_set_affine_nl!(m::GlobalOptimizer{N,T,S}, f::BufferedNonlinearFunction{MC{R,Q}},
                              finite_cut_generated::Bool, safe::Bool) where {N,T<:AbstractFloat,S,R,Q<:RelaxTag}
    if finite_cut_generated && (!safe || is_safe_cut!(m, f.saf))
        lt = LT{T}(-f.saf.constant + m.absolute_constraint_feas_tolerance)
        f.saf.constant = zero(T)
        ci = MOI.add_constraint(m.relaxed_optimizer, f.saf, lt)
        push!(m._buffered_nonlinear_ci, ci)
    end
    return nothing
end

"""
$(TYPEDSIGNATURES)
"""
function relax!(m::GlobalOptimizer{N,T,S}, f::BufferedNonlinearFunction{MC{R,Q}}, safe::Bool) where {N,T<:AbstractFloat,S,R,Q<:RelaxTag}
    evaluator = m._working_problem._relaxed_evaluator

    finite_cut_generated = affine_relax_nonlinear!(f, evaluator, true, true)
    check_set_affine_nl!(m, f, finite_cut_generated, safe)

    finite_cut_generated = affine_relax_nonlinear!(f, evaluator, false, false)
    check_set_affine_nl!(m, f, finite_cut_generated, safe)

    return nothing
end

"""
$(FUNCTIONNAME)

Adds linear objective cut constraint to the `x.relaxed_optimizer`.
"""
function objective_cut!(m::GlobalOptimizer{N,T,S}, saf::Bool) where {N,T<:AbstractFloat,S}
    u = m._global_upper_bound
    if u < Inf
        wp = m._working_problem
        if check_safe && is_safe_cut!(m, wp._objective)
            val = u - wp._objective.constant + m.abs_constraint_feas_tol
            MOI.set(m.mip_solver, MOI.ConstraintSet(), m._obj_ci, LT(val))
        end
    end
    return nothing
end

"""
$(TYPEDSIGNATURES)

A routine that adds relaxations for all nonlinear constraints and quadratic constraints
corresponding to the current node to the relaxed problem. This adds an objective cut
(if specified by `objective_cut_on`) and then sets the `_new_eval_constraint` flag
to false indicating that an initial evaluation of the constraints has occurred. If
the `objective_cut_on` flag is `true` then the `_new_eval_objective` flag is also
set to `false` indicating that the objective expression was evaluated.
"""
function relax_all_constraints!(t::ExtensionType, m::GlobalOptimizer{N,T,S}, q::Int) where {N,T<:AbstractFloat,S}

    wp = _working_problem(m)
    check_safe = (q === 1) ? false : m.cut_safe_on
    wp._relaxed_evaluator.is_first_eval = m._new_eval_constraint

    foreach(x -> relax!(m, x, check_safe), wp._sqf_leq)
    foreach(x -> relax!(m, x, check_safe), wp._sqf_eq)
    foreach(x -> relax!(m, x, check_safe), wp._nonlinear_constr)

    m._new_eval_constraint = false
    objective_cut!(m, check_safe)

    return nothing
end
function relax_constraints!(t::ExtensionType, m::GlobalOptimizer{N,T,S}, q::Int) where {N,T<:AbstractFloat,S}
    relax_all_constraints!(t, m, q)
end
function relax_constraints!(m::GlobalOptimizer{N,T,S}, q::Int) where {N,T<:AbstractFloat,S}
    relax_constraints!(m.ext_type, m, q)
end

"""
$(FUNCTIONNAME)

Deletes all nonlinear constraints added to the relaxed optimizer.
"""
function delete_nl_constraints!(m::GlobalOptimizer)
    MOI.delete.(m.relaxed_optimizer, m._buffered_quadratic_ineq_ci)
    MOI.delete.(m.relaxed_optimizer, m._buffered_quadratic_eq_ci)
    MOI.delete.(m.relaxed_optimizer, m._buffered_nonlinear_ci)
    empty!(m._buffered_quadratic_ineq_ci)
    empty!(m._buffered_quadratic_eq_ci)
    empty!(m._buffered_nonlinear_ci)
    return
end

"""
$(FUNCTIONNAME)

"""
function set_first_relax_point!(m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}
    m._working_problem._relaxed_evaluator.is_first_eval = true
    m._new_eval_constraint = true
    m._new_eval_objective = true
    @__dot__ m._current_xref = _mid(BranchVar, m, 1:N)
    unsafe_check_fill!(isnan, m._current_xref, zero(T), length(m._current_xref))
    return nothing
end


text:
include("HelperFunctions.jl")

export gen_news

"""
    This function calculates changes in news by using a given DFM results  structure. It inputs two datasets, DFM parameters, target time index, and target variable index. The function then produces Nowcast updates and decomposes the changes into news.
     parameters:
        data_old : DataFrame
            old/vintage dataframe, including a date column. Must have same columns as data_new
        data_new : DataFrame
            new data dataframe, including a date column
        output_dfm : Dict
            output of the estimate_dfm function
        target_variable : Symbol
            target variable column name
        target_period : Date
            target period dates
    returns: Dict
       y_old => old nowcast
       y_new => new nowcast
       singlenews => news for each data series
       actual => observed series release values
       forecast => forecasted series values
       weight => news weight
       t_miss => time index for data releases
       v_miss => series index for data releases
       innov => difference between observed and predicted series values ("innovation")
"""
function news_dfm(;old_y, new_y, output_dfm, target_variable, target_period)
    # initializing data
    old_date_col = date_col_name(old_y)
    new_date_col = date_col_name(new_y)
    old_dates = old_y[!, old_date_col]
    new_dates = new_y[!, new_date_col]
    data_old = old_y[!, Not(old_date_col)]
    data_new = new_y[!, Not(new_date_col)]
    # making sure old data has same number of rows/dates as new data
    old_y = leftjoin(DataFrame(date=new_y[!, date_col_name(new_y)]), old_y, on=:date=>date_col_name(old_y))
    data_old = old_y[!, Not(:date)]
    target_index = findall(x-> x == target_variable, Symbol.(names(data_new)))
    target_old = target_new = zeros(1,length(target_index))
    target_period_index = findall(x-> x==target_period, new_dates)[1]

    r = size(output_dfm[:C])[2]
    N = size(data_new)[2]
    singlenews = repeat([0.0], N) |> x-> reshape(x, 1, length(x))

    y_old = y_new = repeat([0.0],1,length(target_index))

    # forecasting
    if !ismissing(data_new[target_period_index, target_index][1]) # no forecast if value for target already exists at target date
        results_old = kalman_filter_constparams(data_old; output_dfm=output_dfm, lag=0)
        # Loop for each target variable
        for i in 1:length(target_index)
            # (Observed value) - (predicted value)
            singlenews[:,target_index[i]] .= data_new[target_period_index, target_index[i]] - results_old[:X_smooth][target_period_index, target_index[i]]
            # Set predicted and observed y values
            y_old[1,i] = results_old[:X_smooth][target_period_index, target_index[i]]
            y_new[1,i] = data_new[target_period_index, target_index[i]]
        end
        # Forecast-related output set to empty
        actual = missing; forecast = missing; weight = missing; t_miss = missing; v_miss = missing; innov = missing; row_miss = missing; col_miss = missing;
    else # yes forecast case, steps A and B
        # Initialize series mean/standard deviation
        means = output_dfm[:means]
        stds = output_dfm[:stds]
        # Calculate indicators for missing values (1 if missing, 0 otherwise)
        miss_old = ismissing.(data_old)
        miss_new = ismissing.(data_new)
        # Indicator for missing--combine above information to single matrix where:
        # (i) -1: Value is in the old data, but missing in new data
        # (ii) 1: Value is in the new data, but missing in old data
        # (iii) 0: Values are either both missing or both available in the datasets
        i_miss = miss_old .- miss_new
        row_miss = [i[1] for i in findall(x->x==1, Array(i_miss))] # Time/variable indices where case (ii) is true
        col_miss = [i[2] for i in findall(x->x==1, Array(i_miss))]

        # Forecast subcase (A): No new information
        if length(col_miss) == 0
            # Fill in missing variables using a Kalman filter
            results_old = kalman_filter_constparams(data_old; output_dfm=output_dfm, lag=0)
            results_new = kalman_filter_constparams(data_new; output_dfm=output_dfm, lag=0)
            # Set predicted and observed y values. New y value is set to old
            y_old = results_old[:X_smooth][target_period_index, target_index]
            y_new = y_old
            # No news, so nothing returned for news-related output
            groupnews = missing; singlenews = missing; gain = missing; gainSer = missing;
            actual = missing; forecast = missing; weight = missing; t_miss = missing; v_miss = missing; innov = missing
        # Forecast subcase (b): new information
        else
            # Difference between forecast time and new data time
            lag = target_period_index .- row_miss
            # Gives biggest time interval between forecast and new data
            k = maximum([abs.(lag); (maximum(lag) - minimum(lag))])
            C =- output_dfm[:C] # Measurement matrix
            R = transpose(output_dfm[:R]) # Covariance for measurement matrix residuals
            # Number of new events
            n_news = length(lag)
            # Smooth old dataset
            results_old = kalman_filter_constparams(data_old; output_dfm=output_dfm, lag=k)
            Plag = results_old[:Plag]
            # Smooth new dataset
            results_new = kalman_filter_constparams(data_new; output_dfm=output_dfm, lag=0)
            # Subset for target variable and forecast time
            y_old = results_old[:X_smooth][target_period_index, target_index]
            y_new = results_new[:X_smooth][target_period_index, target_index]
            Vs = results_old[:Vsmooth][:,:,end-1]
            P1 = missing  # Initialize projection onto updates

            # Cycle through total number of updates
            for i in 1:n_news
                h = abs(target_period_index - row_miss[i])
                m = maximum([row_miss[i], target_period_index])
                # If location of update is later than the forecasting date
                if row_miss[i] > target_period_index
                    Pp = Plag[h + 1][:,:,m]  # P[1:r, h*r+1:h*r+r, m]'
                else
                    Pp = transpose(Plag[h + 1][:,:,m])  # P[1:r, h*r+1:h*r+r, m]
                end
                P1 = !ismissing(P1) ? hcat(P1, Pp * C[col_miss[i], 1:r]) : Pp * C[col_miss[i], 1:r]
            end

            innov = Array{Union{Missing, Float64},2}(undef, length(row_miss), 1)
            for i in 1:length(row_miss)
                # Standardize predicted and observed values
                X_new_norm = (data_new[row_miss[i], col_miss[i]] - means[col_miss[i]]) / stds[col_miss[i]]
                X_sm_norm = (results_old[:X_smooth][row_miss[i], col_miss[i]] - means[col_miss[i]]) / stds[col_miss[i]]
                # Innovation: gives [observed] data - [predicted data]
                innov[i] = X_new_norm - X_sm_norm
            end
            ins = size(innov)[1]
            P2 = missing
            p2 = missing
            WW = repeat([0.0], N, N)

            # Gives non-standardized series weights
            for i in 1:length(lag)
                for j in 1:length(lag)
                    h = abs(lag[i] - lag[j])
                    m = maximum([row_miss[i], row_miss[j]])
                    if row_miss[j] > row_miss[i]
                        Pp = Plag[h + 1][:,:,m]
                    else
                        Pp = transpose(Plag[h + 1][:,:,m])
                    end
                    if (col_miss[i] == col_miss[j]) & (row_miss[i] != row_miss[j])
                        WW[col_miss[i], col_miss[j]] = 0
                    else
                        WW[col_miss[i], col_miss[j]] = R[col_miss[i], col_miss[j]]
                    end
                    p2 = !ismissing(p2) ? hcat(p2, transpose(C[col_miss[i], 1:r]) * Pp * C[col_miss[j], 1:r] + WW[col_miss[i], col_miss[j]]) : transpose(C[col_miss[i], 1:r]) * Pp * C[col_miss[j], 1:r] + WW[col_miss[i], col_miss[j]]
                end
                P2 = vcat(P2, p2)
                p2 = missing
            end
            P2 = P2[.![sum(ismissing.(Array(row))) == length(row) for row in eachrow(P2)],:] # drop all missing rows

            totnews = repeat([0.0], 1, length(target_index))
            temp = Array{Union{Missing, Float64},3}(undef, 1, n_news, length(target_index))
            gain = Array{Union{Missing, Float64},3}(undef, 1, n_news, length(target_index))

            for i in 1:length(target_index) # loop on v_news
                # Convert to real units (unstadardized data)
                totnews[1,i] = (stds[target_index[i]] .* transpose(C[target_index[i], 1:r]) * P1 * inv(P2) * innov)[1]
                temp[1,:,i] = stds[target_index[i]] * transpose(C[target_index[i], 1:r]) * P1 * inv(P2) .* transpose(innov)
                gain[:,:,i] = stds[target_index[i]] * transpose(C[target_index[i], 1:r]) * P1 * inv(P2)
            end

            # Initialize output objects
            singlenews = Array{Union{Missing, Float64},3}(undef, maximum(row_miss) - minimum(row_miss) + 1, N, length(target_index))
            actual = Array{Union{Missing, Float64},2}(undef, N, 1) # Actual forecasted values
            forecast =  Array{Union{Missing, Float64},2}(undef, N, 1) # Forecasted values
            weight = Array{Union{Missing, Float64},3}(undef, N, 1, length(target_index))

            # Fill in output values
            for i in 1:length(innov)
                actual[col_miss[i],1] = data_new[row_miss[i], col_miss[i]]
                forecast[col_miss[i],1] = results_old[:X_smooth][row_miss[i], col_miss[i]]
                for j in 1:length(target_index)
                    singlenews[row_miss[i] - minimum(row_miss) + 1, col_miss[i], j] = temp[1,i,j]
                    weight[col_miss[i],1,j] = (gain[:,i,j] / stds[col_miss[i]])[1]
                end
            end
            singlenews = [sum(skipmissing(col)) for col in eachcol(singlenews[:,:,1])]
            col_miss = unique(col_miss)
        end
    end
    return Dict(
        :y_old => y_old,
        :y_new => y_new,
        :singlenews => singlenews,
        :actual => actual,
        :forecast => forecast,
        :weight => weight,
        :row_miss => row_miss,
        :col_miss => col_miss,
        :innov => innov
    )
end

"""
    This function calculates the DFM nowcast estimates at a current data vintage and compare it with the results from a previous data vintage (pass the same dataset for old_y and new_y to generate an artifical 1-month lagged dataset). The changes are classified as data revisions and data releases, and their impact on the nowcast estimate is presented separately.
     parameters:
        old_y : DataFrame
            DataFrame of older data, including date column
        new_y : DataFrame of newer data
            DataFrame of newest data, including date column
        output_dfm : Dict
            Output/parameters of the estimate_dfm function. No model estimation happens here, this compares new data revisions with pre-estimated parameters.
        target_variable : Symbol
            the target variable column name
        target_period : Dates.Date
            the date for which a forecast is desired
    returns: Dict
        :news_table : DataFrame
           series => the name of the variable
           forecast => the forecast value for the desired time period
           actual => the most recent actually observed value of the variable
           weight => weight of each data release
           impact_releases => impact of data releases on nowcast
           impact_total => impact of data_release + data_revisions on nowcast
           data_release => was this series released between the old and newer datasets
       :y_old : Float
            previous prediction for y
       :y_hat : Float
            prediction for y from new data
"""
function gen_news(;old_y::DataFrame, new_y::DataFrame, output_dfm::Dict, target_variable::Symbol, target_period::Dates.Date)
    # artificial lagged dataset if no old data given
    if isequal(old_y, new_y)
        old_y = create_lag(new_y, 1)
    end
    # making sure old data has same number of rows/dates as new data
    old_y = leftjoin(DataFrame(date=new_y[!, date_col_name(new_y)]), old_y, on=:date=>date_col_name(old_y))

    # add 12 months to each dataset to allow for forecasting
    months_ahead = 12
    last_date = new_y[end, :date]
    row = DataFrame(eachrow(new_y)[1:months_ahead])
    row[!, Not(date_col_name(new_y))] .= missing
    old_y = vcat(old_y, row)
    new_y = vcat(new_y, row)
    old_y[end-months_ahead+1:end, :date] = new_y[end-months_ahead+1:end, :date] = [last_date + Month(i) for i in 1:months_ahead]

    # Update nowcast for target variable 'series' (i) at horizon 'target' (t)
    # Relate nowcast update into news from data releases:
    #   a. Compute the impact from data revisions
    #   b. Compute the impact from new data releases
    data_rev = copy(new_y)
    old_missing = ismissing.(old_y)
    for i in 1:ncol(old_y)
        if eltype(data_rev[!, i]) != Dates.Date
            data_rev[old_missing[!, i], i] .= missing
        end
    end

    # Compute impact from data revisions
    results_old = news_dfm(;old_y=old_y, new_y=data_rev, output_dfm=output_dfm, target_variable=target_variable, target_period=target_period)
    y_old = results_old[:y_old]; old_forecast = results_old[:forecast]

    # Compute impact from data releases
    results_new = news_dfm(;old_y=data_rev, new_y=new_y, output_dfm=output_dfm, target_variable=target_variable, target_period=target_period)
    y_rev = results_new[:y_old]; y_new = results_new[:y_new];
    actual = results_new[:actual]; forecast = results_new[:forecast]; weight = results_new[:weight]

    impact_revisions = y_rev - y_old # Impact from revisions
    news = actual .- forecast # News from releases
    impact_releases = weight[:,:,1] .* news # Impact of releases

    news_table = DataFrame(
        forecast=forecast[:,1],
        actual=actual[:,1],
        weight=weight[:,1],
        impact_releases=impact_releases[:,1],
        impact_total=impact_releases[:,1] .+ impact_revisions
    )
    news_table = DataFrame(Array(news_table) * diagm([100,100,1,100, 100])) |> x-> rename!(x, names(news_table))
    news_table[!, :series] = names(new_y) |> x-> x[.!occursin.(string.(x), String(date_col_name(new_y)))] .|> string
    news_table = news_table[!, [:series; Symbol.(names(news_table)[1:end-1])]]
    data_released = ismissing.(old_y[old_y[!, date_col_name(old_y)] .== target_period, Not(date_col_name(old_y))]) .& # data newly released between two datasets
        .!ismissing.(new_y[new_y[!, date_col_name(new_y)] .== target_period, Not(date_col_name(new_y))]) |> Array |> x-> reshape(x, size(x)[2])
    news_table[!, :data_release] = data_released

    return Dict(
        :news_table => news_table,
        :y_old => y_old,
        :y_new => y_new
    )
end


text:
include("../../src/preprocessing/scaler.jl")
using Test, DataFrames

test_data = DataFrame(
    w = ["one", "two", "three", "four"],
    x = [1, 2, 3, missing],
    y = [15, 45, 90, 15],
    z = [-1, -20, 20, -20],
)

@testset "min_max_scaler" begin
    # all columns
    @test isequal(min_max_scaler(test_data), DataFrame(
        w = ["one", "two", "three", "four"],
        x = [0.0, 0.5, 1.0, missing],
        y = [0.0, 0.4, 1.0, 0.0],
        z = [0.475, 0.0, 1.0, 0.0]
    ))
    # subset of columns
    @test isequal(min_max_scaler(test_data, [:x, :y]), DataFrame(
        w = ["one", "two", "three", "four"],
        x = [0.0, 0.5, 1.0, missing],
        y = [0.0, 0.4, 1.0, 0.0],
        z = [-1, -20, 20, -20]
    ))
    # subset of columns as strings
    @test isequal(min_max_scaler(test_data, ["x", "y"]), DataFrame(
        w = ["one", "two", "three", "four"],
        x = [0.0, 0.5, 1.0, missing],
        y = [0.0, 0.4, 1.0, 0.0],
        z = [-1, -20, 20, -20]
    ))
    # subset of columns + custom scaler
    @test isequal(min_max_scaler(test_data, [:x, :y], scale_min=-100, scale_max=100), DataFrame(
        w = ["one", "two", "three", "four"],
        x = [-100, 0.0, 100.0, missing],
        y = [-100.0, -20.0, 100.0, -100.0],
        z = [-1, -20, 20, -20]
    ))
end

@testset "standard_scaler" begin
    # all columns
    @test isequal(standard_scaler(test_data).w, ["one", "two", "three", "four"])
    @test isequal(standard_scaler(test_data).x, [-1.0, 0.0, 1.0, missing])
    @test isequal(round.(standard_scaler(test_data).y, digits=4), [-0.7406, 0.1058, 1.3754, -0.7406])
    @test isequal(round.(standard_scaler(test_data).z, digits=4), [0.2229, -0.7736, 1.3242, -0.7736])

    # subset of columns
    @test isequal(standard_scaler(test_data, [:x, :y]).w, ["one", "two", "three", "four"])
    @test isequal(standard_scaler(test_data, [:x, :y]).x, [-1.0, 0.0, 1.0, missing])
    @test isequal(round.(standard_scaler(test_data, [:x, :y]).y, digits=4), [-0.7406, 0.1058, 1.3754, -0.7406])
    @test isequal(round.(standard_scaler(test_data, [:x, :y]).z, digits=4), [-1, -20, 20, -20])

    # subset of columns as strings
    @test isequal(standard_scaler(test_data, ["x", "y"]).w, ["one", "two", "three", "four"])
    @test isequal(standard_scaler(test_data, ["x", "y"]).x, [-1.0, 0.0, 1.0, missing])
    @test isequal(round.(standard_scaler(test_data, ["x", "y"]).y, digits=4), [-0.7406, 0.1058, 1.3754, -0.7406])
    @test isequal(round.(standard_scaler(test_data, ["x", "y"]).z, digits=4), [-1, -20, 20, -20])
end


text:
import .FFTW: fft, rfft, plan_fft, plan_rfft

# It's better not to define fft on Colorant arrays, because keeping
# track of the color dimension and the fft-dims is prone to omissions
# or problems due to later operations. So we put the bookkeeping on
# the user, but we try to give helpful suggestions.

function throw_ffterror(f, x, dims=1:ndims(x))
    newdims = plus(dims, channelview_dims_offset(x))
    error("$f not defined for eltype $(eltype(x)). Use channelview, and likely $newdims for the dims in the fft.")
end

for f in (:fft, :rfft)
    pf = Symbol("plan_", f)
    @eval begin
        $f(x::AbstractArray{C}) where {C<:Colorant} = throw_ffterror($f, x)
        $f(x::AbstractArray{C}, dims) where {C<:Colorant} = throw_ffterror($f, x, dims)
        $pf(x::AbstractArray{C}; kws...) where {C<:Colorant} = throw_ffterror($pf, x)
    end
end


text:
include("Include.jl")

function plot_mRNA_C1_data(time_array, simulation_data_array,experimental_data_dictionary)

    (NR,NC) = size(time_array)
    for index = 1:20:NC
        plot(time_array[:,index],1000.0*simulation_data_array[:,index],color="#7AD7F0",alpha=0.80,lw=0.5)
    end

    # Plot the region -
    # calculate the mean and std -
    μ = mean(1000.0*simulation_data_array,dims=2)
    σ = std(1000.0*simulation_data_array,dims=2)
    LB = μ .- (2.58/sqrt(1))*σ
    UB = μ .+ (2.58/sqrt(1))*σ
    fill_between(time_array[:,1], vec(UB), vec(LB), color="#d9d9d9")

    # show std error region -
    LB = μ .- (2.58/sqrt(NC))*σ
    UB = μ .+ (2.58/sqrt(NC))*σ
    fill_between(time_array[:,1], vec(UB), vec(LB), color="#fe9929")

    # Plot mean -
    plot(time_array[:,1],μ,"--",color="#fff7bc",lw=1.5)

    # plot the experimemtal data -
    TEXP = experimental_data_dictionary["c1_ssrA_mRNA_data_table"][!,:Time_in_hr]
    DATA = experimental_data_dictionary["c1_ssrA_mRNA_data_table"][!,:Concentration_in_nM]
    STD = (2.58/sqrt(5))*experimental_data_dictionary["c1_ssrA_mRNA_data_table"][!,:stdev]
    yerr_array = transpose([STD STD])
    errorbar(TEXP, DATA,yerr=yerr_array,fmt="o",mfc="#252525",mec="#252525",color="#252525", lw=1)

    # labels -
    xlabel("Time (hr)",fontsize=16)
    ylabel("Concentration (nM)",fontsize=16)

    # set the face color -
    fig = gcf()
    ax = gca()
    ax.set_facecolor("#F5FCFF")
    fig.set_facecolor("#ffffff")

    # write -
    savefig("./plots/mRNA-C1ssrA-Ensemble.pdf", facecolor=fig.get_facecolor(), edgecolor="none")
end

function main(path_to_simulation_dir::String, path_to_plot_file::String)

    state_index = 5

    # load the experimemtal data -
    exp_dd = load_experimental_data_dictionary("./data")

    # how many files?
    file_name_array = searchdir(path_to_simulation_dir, ".dat")
    number_of_trials = length(file_name_array)

    # initialize -> hardcode the dimension for now
    data_array = zeros(1601,number_of_trials)
    time_array = zeros(1601,number_of_trials)

    # read the simulation dir -
    for (file_index,file_name) in enumerate(file_name_array)

        # load -
        file_path = "$(path_to_simulation_dir)/$(file_name)"
        sim_data_array = readdlm(file_path)

        # what size?
        (NR,NC) = size(sim_data_array)
        for step_index = 1:NR
            data_array[step_index,file_index] = sim_data_array[step_index,(state_index+1)]
            time_array[step_index,file_index] = sim_data_array[step_index,1]
        end
    end

    # plot -
    plot_mRNA_C1_data(time_array, data_array, exp_dd)
end

path_to_simulation_dir = "./simulation_files"
path_to_plot_file = "./plots"
main(path_to_simulation_dir, path_to_plot_file)

max_stars_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_issues_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_forks_repo_name: varnerlab/Biophysical-TXTL-Model-Code


text:
using Zygote, Test, LinearAlgebra

@test gradient(x -> real(abs(x)*exp(im*angle(x))), 10+20im)[1] ≈ 1
@test gradient(x -> imag(real(x)+0.3im), 0.3)[1] ≈ 0
@test gradient(x -> imag(conj(x)+0.3im), 0.3)[1] ≈ -1im
@test gradient(x -> abs((imag(x)+0.3)), 0.3)[1] == 1im

@test gradient(a -> real((a*conj(a))), 0.3im)[1] == 0.6im
@test gradient(a -> real((a.*conj(a))), 0.3im)[1] == 0.6im
@test gradient(a -> real(([a].*conj([a])))[], 0.3im)[1] == 0.6im
@test gradient(a -> real(([a].*conj.([a])))[], 0.3im)[1] == 0.6im
@test gradient(a -> real.(([a].*conj.([a])))[], 0.3im)[1] == 0.6im

@test gradient(x -> norm((im*x) ./ (im)), 2)[1] == 1
@test gradient(x -> norm((im) ./ (im*x)), 2)[1] == -1/4
@test gradient(x -> real(det(x)), [1 2im; 3im 4])[1] ≈ [4 3im; 2im 1]
@test gradient(x -> real(logdet(x)), [1 2im; 3im 4])[1] ≈ [4 3im; 2im 1]/10
@test gradient(x -> real(logabsdet(x)[1]), [1 2im; 3im 4])[1] ≈ [4 3im; 2im 1]/10

# https://github.com/FluxML/Zygote.jl/issues/705
@test gradient(x -> imag(sum(exp, x)), [1,2,3])[1] ≈ real(im .* exp.(1:3))
@test gradient(x -> imag(sum(exp, x)), [1+0im,2,3])[1] ≈ im .* exp.(1:3)

fs_C_to_R = (real,
             imag,
             abs,
             abs2,
             z -> abs(z)*cos(im*angle(z)),
             z->abs(cos(exp(z))),
             z->3*real(z)^3-2*imag(z)^5
             )
@testset "C->R" begin
    for f in fs_C_to_R
        for z in (1.0+2.0im, -2.0+pi*im)
            grad_zygote = gradient(real∘f, z)[1]
            ε = 1e-8
            grad_fd = (f(z+ε)-f(z))/ε + im*(f(z+ε*im)-f(z))/ε
            @test abs(grad_zygote - grad_fd) < sqrt(ε)
        end
    end
end

fs_C_to_C_holomorphic = (cos,
                         exp,
                         log,
                         z->z^2,
                         z->(real(z)+im*imag(z))^2,
                         z->real(z)^2 - imag(z)^2 +2im*(real(z)*imag(z)),
                         z->exp(cos(log(z))),
                         z->abs(z)*exp(im*angle(z)),
                         )
@testset "C->C holomorphic" begin
    for f in fs_C_to_C_holomorphic
        for z in (1.0+2.0im, -2.0+pi*im)
            grad_zygote = gradient(real∘f, z)[1]
            ε = 1e-8
            grad_fd_r = (f(z+ε)-f(z))/ε
            grad_fd_i = (f(z+ε*im)-f(z))/(ε*im)
            @assert abs(grad_fd_r - grad_fd_i) < sqrt(ε) # check the function is indeed holomorphic
            @test abs(grad_zygote - conj(grad_fd_r)) < sqrt(ε)
        end
    end
end

fs_C_to_C_non_holomorphic = (conj,
                             z->abs(z)+0im,
                             z->im*abs(z),
                             z->abs2(z)+0im,
                             z->im*abs2(z),
                             z->z'z,
                             z->conj(z)*z^2,
                             )
@testset "C->C non-holomorphic" begin
    for f in (fs_C_to_C_holomorphic...,fs_C_to_C_holomorphic...)
        for z in (1.0+2.0im, -2.0+pi*im)
            grad_zygote = gradient(real∘f, z)[1]
            ε = 1e-8
            grad_fd = real(f(z+ε)-f(z))/ε + im*real(f(z+ε*im)-f(z))/ε
            @test abs(grad_zygote - grad_fd) < sqrt(ε)
        end
    end
end


text:
########################################
## File Name: forward_backward_simulation.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/14
## Description: Forward Backward Simulation for SACBP
########################################

abstract type SACSimulationModel end

struct SimulatePosRangeLocalization2D <: SACSimulationModel
    continuousTransModel::TransModel_Pos
    discreteTransModel::UKFPosRange
    costModel::CostPosRangeLocalization
    dtc::Float64
    dto::Float64
    dtexec::Float64
    Q::Matrix{Float64}
    Cu::Matrix{Float64}
end
function SimulatePosRangeLocalization2D(continuousTransModel::TransModel_Pos,
                                        discreteTransModel::UKFPosRange,
                                        dtc::Float64,dto::Float64,dtexec::Float64,
                                        Q::Matrix{Float64},Cu::Matrix{Float64})
    Int64(dto/dtc);  # Inexact Error Check.
    return SimulatePosRangeLocalization2D(continuousTransModel,discreteTransModel,CostPosRangeLocalization(),dtc,dto,dtexec,Q,Cu)
end
function simulateForward(simModel::SimulatePosRangeLocalization2D,
                         s_init::AugState{<:PhysPos,<:VecBelMvNormal},
                         UArray::Vector{<:PosControl},
                         rng::MersenneTwister)
    SArray = Vector{typeof(s_init)}(undef, length(UArray)+1);
    SArray_before_update = typeof(s_init)[];
    YArray = Vector{Float64}[];
    SArray[1] = s_init;
    o_c_ratio = Int64(simModel.dto/simModel.dtc);
    kk = 0;
    for t = 1:length(UArray)
        kk += 1;
        p = trans(simModel.continuousTransModel,SArray[t].p,UArray[t],simModel.dtc);
        bVec = similar(s_init.b);
        for ii = 1:length(bVec)
            bVec[ii] = BelMvNormal(p.t,SArray[t].b[ii].params)
        end
        SArray[t+1] = AugState(p.t,p,bVec);
        if kk == o_c_ratio
            push!(SArray_before_update,SArray[t]);
            yVec,SArray[t+1] = trans(simModel.discreteTransModel,SArray[t+1],simModel.dto,simModel.Q,rng);
            push!(YArray,yVec);
            kk = 0;
        end
    end
    YArray,SArray,SArray_before_update
end
function simulateForward(simModel::SimulatePosRangeLocalization2D,
                         nominalPolicy::GradientMultiTargetLocalizationPolicy,
                         s_init::AugState{<:PhysPos,<:VecBelMvNormal},
                         u_param_min::Vector{<:Real},
                         u_param_max::Vector{<:Real},
                         UArray::Vector{<:PosControl},
                         rng::MersenneTwister)
    UArrayClosed = similar(UArray);
    SArray = Vector{typeof(s_init)}(undef, length(UArray)+1);
    SArray_before_update = typeof(s_init)[];
    YArray = Vector{Float64}[];
    SArray[1] = s_init;
    UArrayClosed[1] = control_nominal(nominalPolicy,simModel.discreteTransModel,s_init,u_param_min,u_param_max,simModel.dto,simModel.Q);
    o_c_ratio = Int64(simModel.dto/simModel.dtc);
    kk = 0;
    for t = 1:length(UArray)
        kk += 1;
        p = trans(simModel.continuousTransModel,SArray[t].p,UArrayClosed[t],simModel.dtc);
        bVec = similar(s_init.b);
        for ii = 1:length(bVec)
            bVec[ii] = BelMvNormal(p.t,SArray[t].b[ii].params)
        end
        SArray[t+1] = AugState(p.t,p,bVec);
        if t != length(UArray)
            UArrayClosed[t+1] = PosControl(UArray[t+1].t,vec(UArrayClosed[t]));
        end
        if kk == o_c_ratio
            push!(SArray_before_update,SArray[t]);
            yVec,SArray[t+1] = trans(simModel.discreteTransModel,SArray[t+1],simModel.dto,simModel.Q,rng);
            if t != length(UArray)
                UArrayClosed[t+1] = control_nominal(nominalPolicy,simModel.discreteTransModel,SArray[t+1],u_param_min,u_param_max,simModel.dto,simModel.Q);
            end
            push!(YArray,yVec);
            kk = 0;
        end
    end
    YArray,SArray,SArray_before_update,UArrayClosed
end
function simulateBackward(simModel::SimulatePosRangeLocalization2D,
                          SArray::Vector{<:AugState{<:PhysPos,<:VecBelMvNormal}},
                          SArray_before_update::Vector{<:AugState{<:PhysPos,<:VecBelMvNormal}},
                          UArray::Vector{<:PosControl},
                          YArray::Vector{Vector{Float64}})
    CoSArray = similar(SArray);
    coSP_end = PhysPos(SArray[end].t,termCost_grad_p(simModel.costModel,SArray[end]));
    coSB_end = similar(SArray[end].b);
    for ii = 1:length(coSB_end)
        coSB_end[ii] = BelMvNormal(SArray[end].t,termCost_grad_b(simModel.costModel,SArray[end])[ii]);
    end
    CoSArray[end] = AugState(coSP_end,coSB_end);
    o_c_ratio = Int64(simModel.dto/simModel.dtc);
    kk = o_c_ratio;
    ll = length(SArray_before_update);
    for t = Iterators.reverse(1:length(UArray))
        if kk == o_c_ratio
            CoSArray[t] = cotrans(simModel.discreteTransModel,simModel.costModel,CoSArray[t+1],SArray_before_update[ll],YArray[ll],simModel.dto,simModel.Q);
            ll -= 1;
            kk = 1;
        else
            CoSArray[t] = cotrans(simModel.continuousTransModel,simModel.costModel,CoSArray[t+1],SArray[t+1],UArray[t+1],simModel.dtc);
            kk += 1;
        end
    end
    return CoSArray
end

function evaluateCost(simModel::SimulatePosRangeLocalization2D,
                      s_init::AugState{<:PhysPos,<:VecBelMvNormal},
                      UArray::Vector{<:PosControl},
                      rng::MersenneTwister)
    YArray,SArray,SArray_before_update = simulateForward(simModel,s_init,UArray,rng);
    C = termCost(simModel.costModel,SArray[end]);
    for t = 1:length(SArray)-1;
        C += instCost(simModel.costModel,SArray[t],UArray[t],simModel.Cu)*simModel.dtc
    end
    return C
end

struct SimulateManipulate2D  <: SACSimulationModel
    continuousTransModel::CPredictManipulate2D
    discreteTransModel::DUpdateManipulate2D
    costModel::CostManipulate2D
    dtc::Float64       # Discrete Time Step for Continuous Simulation
    dto::Float64       # Discrete Interval for Observations.
    dtexec::Float64
    Q::Matrix{Float64} # Symmetric PSD Matrix. Transition Covariance.
    R::Matrix{Float64} # Symmetric PSD Matrix. Observation Covariance.
    Cs::Matrix{Float64} # Symmetric PSD Matrix. State Cost Coefficient Matrix.
    Cu::Matrix{Float64} # Symmetric PD Matrix. Control Cost Coefficient Matrix
end
function SimulateManipulate2D(dtc::Real,dto::Real,dtexec::Float64,
                              Q::Matrix{Float64},R::Matrix{Float64},Cs::Matrix{Float64},Cu::Matrix{Float64})
    cTransModel = CPredictManipulate2D();
    dTransModel = DUpdateManipulate2D();
    costModel = CostManipulate2D();
    Int64(dto/dtc) # Inexact Error Check.
    return SimulateManipulate2D(cTransModel,dTransModel,costModel,dtc,dto,dtexec,Q,R,Cs,Cu)
end
function simulateForward(simModel::SimulateManipulate2D,s_init::BelMvNormal,UArray::Vector{<:MControl2D},rng::MersenneTwister)
    SArray = Vector{typeof(s_init)}(undef, length(UArray)+1);
    SArray_before_update = typeof(s_init)[];
    YArray = Vector{Float64}[];
    SArray[1] = s_init;
    o_c_ratio = Int64(simModel.dto/simModel.dtc);
    kk = 0;
    for t = 1:length(UArray)
        #d = Distributions.MvNormal(SArray[t]);
        kk += 1;
        #println("In simulateForward(): t = $(SArray[t].t), $(SArray[t].params[SArray[t].dim+1:end])")
        SArray[t+1] = trans(simModel.continuousTransModel,SArray[t],UArray[t],simModel.dtc,simModel.Q);
        try d = Distributions.MvNormal(SArray[t+1])
        catch
            println("In simulateForward() 1: t = $(SArray[t+1].t), $(SArray[t+1].params[SArray[t+1].dim+1:end])")
        end
        if kk == o_c_ratio
            push!(SArray_before_update,SArray[t]);
            y,SArray[t+1] = trans(simModel.discreteTransModel,SArray[t+1],MControl2D(SArray[t+1].t,vec(UArray[t])),simModel.R,rng);
            try d = Distributions.MvNormal(SArray[t+1])
            catch
                println("In simulateForward() 2: t = $(SArray[t+1].t), $(SArray[t+1].params[SArray[t+1].dim+1:end])")
            end
            push!(YArray,y);
            kk = 0;
        end
    end
    return YArray,SArray,SArray_before_update
end
function simulateForward(simModel::SimulateManipulate2D,
                         nominalPolicy::ManipulatePositionControlPolicy,
                         s_init::BelMvNormal,
                         u_param_min::Vector{<:Real},
                         u_param_max::Vector{<:Real},
                         posGain,
                         rotGain,
                         UArray::Vector{<:MControl2D},
                         rng::MersenneTwister)
    SArray = Vector{typeof(s_init)}(undef, length(UArray)+1);
    SArray_before_update = typeof(s_init)[];
    YArray = Vector{Float64}[];
    UArrayClosed = similar(UArray);
    SArray[1] = s_init;
    UArrayClosed[1] = control_nominal(nominalPolicy,s_init,u_param_min,u_param_max,posGain,rotGain);
    o_c_ratio = Int64(simModel.dto/simModel.dtc);
    kk = 0;
    for t = 1:length(UArray)
        #d = Distributions.MvNormal(SArray[t]);
        kk += 1;
        #println("In simulateForward(): t = $(SArray[t].t), $(SArray[t].params[SArray[t].dim+1:end])")
        SArray[t+1] = trans(simModel.continuousTransModel,SArray[t],UArrayClosed[t],simModel.dtc,simModel.Q);
        if t != length(UArray)
            UArrayClosed[t+1] = MControl2D(UArray[t+1].t,vec(UArrayClosed[t]));
        end
        try d = Distributions.MvNormal(SArray[t+1])
        catch
            println("In simulateForward() 1: t = $(SArray[t+1].t), $(SArray[t+1].params[SArray[t+1].dim+1:end])")
        end
        if kk == o_c_ratio
            push!(SArray_before_update,SArray[t]);
            y,SArray[t+1] = trans(simModel.discreteTransModel,SArray[t+1],MControl2D(SArray[t+1].t,vec(UArrayClosed[t])),simModel.R,rng);
            if t != length(UArray)
                UArrayClosed[t+1] = control_nominal(nominalPolicy,SArray[t+1],u_param_min,u_param_max,posGain,rotGain);
            end
            try d = Distributions.MvNormal(SArray[t+1])
            catch
                println("In simulateForward() 2: t = $(SArray[t+1].t), $(SArray[t+1].params[SArray[t+1].dim+1:end])")
            end
            push!(YArray,y);
            kk = 0;
        end
    end
    return YArray,SArray,SArray_before_update,UArrayClosed
end
function simulateBackward(simModel::SimulateManipulate2D,SArray::VecBelMvNormal,SArray_before_update::VecBelMvNormal,
                          UArray::Vector{<:MControl2D},YArray::Vector{Vector{Float64}})
    CoSArray = similar(SArray);
    CoSArray[end] = BelMvNormal(SArray[end].t,termCost_grad(simModel.costModel,SArray[end],simModel.Cs));
    o_c_ratio = Int64(simModel.dto/simModel.dtc);
    kk = o_c_ratio;
    ll = length(SArray_before_update);
    for t = Iterators.reverse(1:length(UArray))
        if kk == o_c_ratio
            CoSArray[t] = cotrans(simModel.discreteTransModel,simModel.costModel,CoSArray[t+1],SArray_before_update[ll],UArray[t],YArray[ll],simModel.R);
            ll -= 1;
            kk = 1;
        else
            CoSArray[t] = cotrans(simModel.continuousTransModel,simModel.costModel,CoSArray[t+1],SArray[t+1],UArray[t+1],simModel.dtc,simModel.Q,simModel.Cs);
            kk += 1;
        end
    end
    return CoSArray
end

function evaluateCost(simModel::SimulateManipulate2D,
                      s_init::BelMvNormal,
                      UArray::Vector{<:MControl2D},
                      rng::MersenneTwister)
    YArray,SArray,SArray_before_update = simulateForward(simModel,s_init,UArray,rng);
    C = termCost(simModel.costModel,SArray[end],simModel.Cs);
    for t = 1:length(SArray)-1;
        C += instCost(simModel.costModel,SArray[t],UArray[t],simModel.Cs,simModel.Cu)*simModel.dtc
    end
    return C
end


text:
using ANBIMA
using Dates
using BusinessDays
using Test

dt = Date(2011,1,3)
@info "Reading IMAs on $dt"
data = ANBIMA.read_IMA(dt)
for (k,v) in data.elements
    println("  $k: $(v.index)")
end

dt = BusinessDays.advancebdays(:BRSettlement, Dates.today(), -2)
@info "Reading IDAs on $dt"
data = ANBIMA.read_IDA(dt)
for (k,v) in data.elements
    println("  $k: $(v.index)")
end
dt_old = Date(2020,9,1)
@test_throws AssertionError ANBIMA.read_IDA(dt_old)

@info "Reading ETTJs on $dt"
data = ANBIMA.read_ettj(dt)
vertices = [1, 252,504, 10000]
for (k,v) in data
	taxas = [t => 100*ANBIMA.zerorate(v, t) for t in vertices]
	println("  $k")
	for x in taxas
		t = first(x)
		r = last(x)
		println("    $t: $r")
	end
end
@test_throws AssertionError ANBIMA.read_ettj(dt_old)

@info "Reading credit spread curves on $dt"
data = ANBIMA.read_credit_spread_data(dt)
curves = ANBIMA.calibrate_credit_spread_curves(dt)
for (k,v) in curves
	println("$k => $v")
end


text:
#
# diff-tests.jl -
#
# Tests for finite difference operators.
#
module TestingLazyAlgebraDiff

using Test
using MayOptimize
using LazyAlgebra
using LazyAlgebra.Foundations
using LazyAlgebra.FiniteDifferences: limits, ArrayAxis, ArrayAxes,
    differentiation_order, dimensions_of_interest, optimization_level

include("common.jl")

@inline colons(n::Integer) = ntuple(x -> Colon(), max(Int(n), 0))

# Slice along the last dimension.
slice(A::AbstractArray{T,N}, d::Int) where {T,N} =
    view(A, colons(N - 1)..., d)

# Reference implementation of finite differences.

gram_diff_ref(order::Integer, A::AbstractArray; dim::Int=1) =
    diff_ref(order, diff_ref(order, A; dim=dim, adj=false); dim=dim, adj=true)

diff_ref(order::Integer, A::AbstractArray; kwds...) =
    (order == 1 ? diff1_ref(A; kwds...) :
     order == 2 ? diff2_ref(A; kwds...) :
     error("invalid derivative order ", order))

diff1_ref(A::AbstractArray; kwds...) = diff1_ref!(similar(A), A; kwds...)

diff2_ref(A::AbstractArray; kwds...) = diff2_ref!(similar(A), A; kwds...)

function diff1_ref!(dst::AbstractArray{<:Any,N},
                    src::AbstractArray{<:Any,N};
                    dim::Int=1,
                    adj::Bool=false) where {N}
    inds = axes(src)
    @assert axes(dst) == inds
    @assert 1 ≤ dim ≤ N
    I = colons(dim - 1)
    J = inds[dim]
    K = colons(N - dim)
    diff1_ref!(dst, src, I, J, K, adj)
end

function diff2_ref!(dst::AbstractArray{<:Any,N},
                    src::AbstractArray{<:Any,N};
                    dim::Int=1,
                    adj::Bool=false) where {N}
    inds = axes(src)
    @assert axes(dst) == inds
    @assert 1 ≤ dim ≤ N
    I = colons(dim - 1)
    J = inds[dim]
    K = colons(N - dim)
    diff2_ref!(dst, src, I, J, K)
end

#
# Code to apply 1st order finite difference operator D or its adjoint D' (given
# below) along a given dimensions of an array.
#
#     D = [ -1   1   0   0
#            0  -1   1   0
#            0   0  -1   1
#            0   0   0   0];
#
#     D' = [ -1   0   0   0
#             1  -1   0   0
#             0   1  -1   0
#             0   0   1   0];
#
function diff1_ref!(dst::AbstractArray{<:Any,N},
                    src::AbstractArray{<:Any,N},
                    I::Tuple{Vararg{Colon}},
                    J::ArrayAxis,
                    K::Tuple{Vararg{Colon}},
                    adj::Bool) where {N}
    # Forward 1st order finite differences assuming flat boundary conditions
    # along a given dimension.
    len = length(J)
    if len > 1
        j_first = first(J)
        j_last = last(J)
        if adj
            map!(-,
                 view(dst, I..., j_first, K...),
                 view(src, I..., j_first, K...))
            if len > 2
                map!(-,
                     view(dst, I..., (j_first+1):(j_last-1), K...),
                     view(src, I..., (j_first  ):(j_last-2), K...),
                     view(src, I..., (j_first+1):(j_last-1), K...))
            end
            copyto!(view(dst, I..., j_last,   K...),
                    view(src, I..., j_last-1, K...))
        else
            map!(-,
                 view(dst, I..., (j_first  ):(j_last-1), K...),
                 view(src, I..., (j_first+1):(j_last  ), K...),
                 view(src, I..., (j_first  ):(j_last-1), K...))
            fill!(view(dst, I..., j_last, K...), zero(eltype(dst)))
        end
    else
        fill!(dst, zero(eltype(dst)))
    end
    return dst
end

#
# 2nd order finite differences D2 with flat boundary conditions is given by:
#
#     D2 = [-1   1   0   0
#            1  -2   1   0
#            0   1  -2   1
#            0   0   1  -1]
#
# Note that this operator is self-adjoint, so there is no `adj` argument.
#
function diff2_ref!(dst::AbstractArray{<:Any,N},
                    src::AbstractArray{<:Any,N},
                    I::Tuple{Vararg{Colon}},
                    J::ArrayAxis,
                    K::Tuple{Vararg{Colon}}) where {N}
    # Forward 1st order finite differences assuming flat boundary conditions
    # along a given dimension.
    len = length(J)
    if len ≥ 2
        j_first = first(J)
        j_last = last(J)
        map!((b, c) -> c - b,
             view(dst, I..., j_first,   K...),
             view(src, I..., j_first,   K...),
             view(src, I..., j_first+1, K...))
        if len > 2
            map!((a, b, c) -> (c - b) - (b - a),
                 view(dst, I..., (j_first+1):(j_last-1), K...),
                 view(src, I..., (j_first  ):(j_last-2), K...),
                 view(src, I..., (j_first+1):(j_last-1), K...),
                 view(src, I..., (j_first+2):(j_last  ), K...))
        end
        map!((a, b) -> a - b,
             view(dst, I..., j_last,   K...),
             view(src, I..., j_last-1, K...),
             view(src, I..., j_last,   K...))
    else
        fill!(dst, zero(eltype(dst)))
    end
    return dst
end

@testset "Finite differences" begin
    # First test the correctness of the result compared to the above reference
    # implementation and use the saclar product to check the adjoint.  Use
    # integer values for exact computations, over a reduced range to avoid
    # overflows and have exact floating-point representation.
    vmin = -(vmax = Float64(7_000))
    vals = vmin:vmax
    @testset "Differentiation order = $order" for order in 1:2
        # Dimensions to test: 1, 2 and any ≥ 3 for 1st order derivatives, 1:4
        # and any ≥ 5 for 2nd order derivatives.
        dimlist = (order == 1 ? (1, 2, 5) : ((1:4)..., 6))
        sizes = Tuple{Vararg{Int}}[]
        for dim1 in dimlist
            push!(sizes, (dim1,))
        end
        for dim1 in dimlist, dim2 in dimlist
            push!(sizes, (dim1,dim2))
        end
        for dim1 in dimlist, dim2 in dimlist, dim3 in dimlist
            push!(sizes, (dim1,dim2,dim3))
        end
        for dims in sizes
            # Random values (an copies to check that input variables do not
            # change).
            x = rand(vals, dims)
            y = rand(vals, (size(x)..., ndims(x)))
            xsav = copy(x)
            ysav = copy(y)
            atol = zero(eltype(x))
            rtol = 16*eps(eltype(x))
            # Apply along all dimensions specified as a colon.
            D_all = Diff(order,:,Debug)
            @test differentiation_order(D_all) === order
            @test dimensions_of_interest(D_all) === Colon
            @test optimization_level(D_all) === Debug
            D_all_x = D_all*x;   @test x == xsav
            Dt_all_y = D_all'*y; @test y == ysav
            @test vdot(x, Dt_all_y) == vdot(y, D_all_x)
            test_api(Direct, D_all, x, y; atol=atol, rtol=rtol)
            test_api(Adjoint, D_all, x, y; atol=atol, rtol=rtol)
            DtD_all = D_all'*D_all
            @test DtD_all === gram(D_all)
            @test isa(DtD_all, Gram{typeof(D_all)})
            @test differentiation_order(DtD_all) === order
            @test dimensions_of_interest(DtD_all) === dimensions_of_interest(D_all)
            @test optimization_level(DtD_all) === optimization_level(D_all)
            @test DtD_all*x == D_all'*D_all_x; @test x == xsav
            test_api(Direct, DtD_all, x, x; atol=atol, rtol=rtol)
            # Apply along all dimensions specified as a list.
            D_lst = Diff(order,((1:ndims(x))...,),Debug)
            @test differentiation_order(D_lst) === order
            @test dimensions_of_interest(D_lst) === ((1:ndims(x))...,)
            @test optimization_level(D_lst) === Debug
            D_lst_x = D_lst*x;   @test x == xsav
            Dt_lst_y = D_lst'*y; @test y == ysav
            @test vdot(x, Dt_lst_y) == vdot(y, D_lst_x)
            test_api(Direct, D_lst, x, y; atol=atol, rtol=rtol)
            test_api(Adjoint, D_lst, x, y; atol=atol, rtol=rtol)
            DtD_lst = D_lst'*D_lst
            @test DtD_lst === gram(D_lst)
            @test isa(DtD_lst, Gram{typeof(D_lst)})
            @test differentiation_order(DtD_lst) === order
            @test dimensions_of_interest(DtD_lst) === dimensions_of_interest(D_lst)
            @test optimization_level(DtD_lst) === optimization_level(D_lst)
            @test DtD_lst*x == D_lst'*D_lst_x
            test_api(Direct, DtD_lst, x, x; atol=atol, rtol=rtol)
            if ndims(x) > 1
                # Apply along all dimensions in reverse order specified as a
                # range and to avoid all dimensions at a time version.
                D_rev = Diff(order,ndims(x):-1:1,Debug)
                @test differentiation_order(D_rev) === order
                @test dimensions_of_interest(D_rev) === ((ndims(x):-1:1)...,)
                @test optimization_level(D_rev) === Debug
                D_rev_x = D_rev*x;   @test x == xsav
                Dt_rev_y = D_rev'*y; @test y == ysav
                @test vdot(x, Dt_rev_y) == vdot(y, D_rev_x)
                test_api(Direct, D_rev, x, y; atol=atol, rtol=rtol)
                test_api(Adjoint, D_rev, x, y; atol=atol, rtol=rtol)
                DtD_rev = D_rev'*D_rev
                @test DtD_rev === gram(D_rev)
                @test isa(DtD_rev, Gram{typeof(D_rev)})
                @test differentiation_order(DtD_rev) === order
                @test dimensions_of_interest(DtD_rev) === dimensions_of_interest(D_rev)
                @test optimization_level(DtD_rev) === optimization_level(D_rev)
                @test DtD_rev*x == D_rev'*D_rev_x
                test_api(Direct, DtD_rev, x, x; atol=atol, rtol=rtol)
            end
            for d in 1:ndims(x)
                # Apply along a single dimension.
                D_one = Diff(order,d,Debug)
                @test differentiation_order(D_one) === order
                @test dimensions_of_interest(D_one) === d
                @test optimization_level(D_one) === Debug
                z = slice(y, d)
                D_ref_x = diff_ref(order, x, dim=d, adj=false)
                Dt_ref_z = diff_ref(order, z, dim=d, adj=true)
                @test D_one*x == D_ref_x;   @test x == xsav
                @test D_one'*z == Dt_ref_z; @test y == ysav # z is a view on y
                @test slice(D_all_x, d) == D_ref_x
                @test slice(D_lst_x, d) == D_ref_x
                if ndims(x) > 1
                    r = ndims(x) - d + 1
                    @test slice(D_rev_x, r) == D_ref_x
                end
                test_api(Direct, D_one, x, z; atol=atol, rtol=rtol)
                test_api(Adjoint, D_one, x, z; atol=atol, rtol=rtol)
                DtD_one = D_one'*D_one
                @test DtD_one === gram(D_one)
                @test isa(DtD_one, Gram{typeof(D_one)})
                @test differentiation_order(DtD_one) === order
                @test dimensions_of_interest(DtD_one) === dimensions_of_interest(D_one)
                @test optimization_level(DtD_one) === optimization_level(D_one)
                DtD_one_x = DtD_one*x; @test x == xsav
                @test DtD_one_x == gram_diff_ref(order, x, dim=d)
                @test DtD_one_x == D_one'*(D_one*x)
                test_api(Direct, DtD_one, x, x; atol=atol, rtol=rtol)
            end
        end
    end
end

nothing

end # module


text:
using SparseArrays

include("PushRelabelMaxflow.jl")
# include("GurobiMaxflow.jl")
include("pwl_approx.jl")
function TwoNodeEdge(w)
    # a = weight from each node to the source
    # b = weight from one node to the otehr
    #c = weight from each node to the sink
    a = w[1]/2
    b = w[2] - w[1]/2 - w[3]/2
    c = w[3]/2

    return a,b,c
end

function SparseCard(Edges,EdgesW,n,epsilon;refined=true,flowsolver="pushrelabel",gurobioutput=0)

    """
    Wrapper for SparseCard, for solving a Card-DSFM problem.

        refined = true gives better minimal approximate reductions
                (same density of graph, better reductions)
        flowsolver = choose from push-relabel or solving LP with Gurobi
    """

    # Reduce to a graph s-t cut problem
    tic = time()
    A, svec, tvec, worst_eps = CardDSFM_reduction(Edges,EdgesW,n,epsilon,refined)
    reducetime = time()-tic
    N = size(A,1)

    # Solve with a max-flow solver
    if flowsolver == "pushrelabel"

        tic = time()
        F = maxflow(A,svec,tvec,0)
        flowtime = time()-tic
        flowval = F.flowvalue

        # get source set, requires removing source node s = 1
        S_in_A = source_nodes_min(F,1e-10)[2:end].-1

        # ignore auxiliary nodes
        S = intersect(1:n,S_in_A)
        eS = zeros(n)
        eS[S] .= 1
        cutS = eval_fS(Edges,EdgesW,eS)

    else

        Svec, flowval,f,flowtime = Gubori_MaxFlow(A,svec,tvec,gurobioutput)
        eS = Svec[1:n]

    end
    cutval = eval_fS(Edges,EdgesW,eS)

    # The value of the flow in the approximate graph approximates the
    approxcut = flowval
    return eS, cutval, approxcut, flowtime, reducetime, worst_eps
end

function SymmetricCard_reduction(Edges,EdgesW,n,epsilon,returnIJV=true)
    """
    Reduce a hypergraph with symmetric cardinality-based submodular
    splitting functions to a directed graph preserving cuts.

    Edges[i] gives nodes in the ground set of function i in the sum
    EdgesW[i] gives the corresponding CB function penalties
        EdgesW[i][j] gives function evaluation at input (j-1)

    NOTE:

    N = number of objects in the ground set
    epsilon[i] = approximation guarantee for reducing a function/hyperedge of size i
    """
    if length(epsilon) == 1
        K = maximum(length(e) for e in Edges)
        epsilon = epsilon*ones(K)
    end
    N = n
    I = Vector{Int64}()
    J = Vector{Int64}()
    W = Vector{Float64}()
    for i = 1:length(Edges)
        e = Edges[i]
        w = EdgesW[i]
        k = length(e)
        @assert(is_increasing(w))
        @assert(check_cb_submodular(w))
        # @assert(length(w) == floor(Int64,(k)/2)+1)
        # @assert(w[1] == 0)
        if w[1] != 0
            println("For symmetric splitting functions, include the zero: w[1] = w_0 = 0.")
            return
        end
        if length(w) != floor(Int64,(k)/2)+1
            println("Give one penalty for each possible small side of the cut.")
            wrong = length(w)
            right = floor(Int64,(k)/2)+1
            println("length(w) should be $right, but it's $wrong")
            return
        end

        if k == 2
            # just an edge
            push!(I,e[1])
            push!(J,e[2])
            push!(W,w[2])
            push!(I,e[2])
            push!(J,e[1])
            push!(W,w[2])

        elseif k == 3
            # just a triangle
            for ii = 1:3
                for jj = 1:3
                    if ii != jj
                        push!(I,e[ii])
                        push!(J,e[jj])
                        push!(W,w[2]/2)
                    end
                end
            end
        else
            a, b = SymmetricSCB_to_Gadget(w,epsilon[k])

            @assert(minimum(a) > 0)
            @assert(minimum(b) > 0)

            L = length(a)
            for j = 1:L
                # add two new auxiliary nodes for this gadget
                push!(I,N+1)
                push!(J,N+2)
                push!(W,a[j]*b[j])
                for v in e
                    push!(I,v)
                    push!(J,N+1)
                    push!(W,a[j])
                    push!(I,N+2)
                    push!(J,v)
                    push!(W,a[j])
                end
                N += 2
            end
        end
    end

    A = sparse(I,J,W,N,N)
    if returnIJV
        return A,I,J,W
    else
        return A
    end
end

function CardDSFM_reduction(Edges,EdgesW,n,epsilon,refined=true)
    """
    Reduce an instance of Card-DSFM to a graph where the s-t cut
    structure approximately models the cardinality-based decomposable
    submodular function.

    A minimum s-t cut in the graph will then solve the problem.

    Think of the problem as a hypergraph cut problem.

    Edges[i] gives nodes in the ground set of function i in the sum
    EdgesW[i] gives the corresponding CB function penalties
        EdgesW[i][j] gives function evaluation at input (j-1)
    N = number of objects in the ground set
    epsilon[i] = approximation guarantee for reducing a function/hyperedge of size i
    refined = true takes a tiny bit long, but gives better gadget reductions.
        (it finds a better approximation, AND a minimum number of gadgets)
    """
    N = n
    I = Vector{Int64}()
    J = Vector{Int64}()
    W = Vector{Float64}()
    svec = zeros(n)
    tvec = zeros(n)
    if refined
        worst_eps = 0
    else
        worst_eps = maximum(epsilon)
    end
    for i = 1:length(Edges)
        e = Edges[i]
        w = EdgesW[i]
        k = length(e)
        @assert(check_cb_submodular(w))
        if k == 2
            z0,uv,zk = TwoNodeEdge(w)
            svec[e[1]] += z0
            svec[e[2]] += z0
            tvec[e[1]] += zk
            tvec[e[2]] += zk
            push!(I,e[1])
            push!(J,e[2])
            push!(W,uv)
            push!(I,e[2])
            push!(J,e[1])
            push!(W,uv)
        else
            if refined
                # don't just find minimum number of CB-gadgets,
                # find best approx with minimum number of gadgets
                z0, zk, a, b, cgf, eps_val = Refined_SCB_to_CGF(w,epsilon[k])
                if eps_val > worst_eps
                    worst_eps = eps_val
                end
            else
                z0, zk, a, b, cgf = SubCardFun_to_CGF_weights(w,epsilon[k],false)
            end
            @assert(min(z0,zk) > -1e-14)
            z0 = max(z0,0)
            zk = max(zk,0)
            L = length(a)
            if z0 < 0
                @show z0
            end
            if zk < 0
                @show zk
                @show w[end]
            end
            if length(a) > 0 && minimum(a) < 0
                @show minimum(a)
            end
            if length(b) > 0 && minimum(b) < 0
                @show minimum(b)
            end
            for v in e
                svec[v] += z0
                tvec[v] += zk
            end
            for j = 1:L
                for v in e
                    push!(I,v)
                    push!(J,N+j)
                    push!(W,a[j]*(k-b[j]))
                    push!(I,N+j)
                    push!(J,v)
                    push!(W,a[j]*b[j])
                end
            end
            N += L
        end
    end

    A = sparse(I,J,W,N,N)
    svec = [svec;zeros(N-n)]
    tvec = [tvec;zeros(N-n)]
    return A, svec, tvec, worst_eps
end

function eval_fS(Edges::Vector{Vector{Int64}},EdgesW::Vector{Vector{Float64}},eS::Vector)
    """
    Evaluate the cardinality-based decomposable submodular function f at
    set S, where the function is given by ground sets in Edges, each
    with the cardinality-based penalties given in EdgesW.
    """
    cutval = 0
    for i = 1:length(Edges)
        e = Edges[i]
        w = EdgesW[i]
        # t = length(intersect(e,S))
        t = round(Int64,sum(eS[e]))
        cutval += w[t+1]
    end

    return cutval
end

function eval_min_dircut(A,svec,tvec,S,T)
    """
    A is the adjacency matrix between non terminal nodes.
    svec and tvec give terminal edges.
    Fix S to be on the source side, and T to be on the sink side.
    Find the minimum S-T cut with these nodes anchored in this way.
    """
    svec = copy(svec)
    tvec = copy(tvec)
    svec[S] .= 2*sum(A)
    tvec[T] .= 2*sum(A)
    F = maxflow(A,svec,tvec,0)
    S = F.source_nodes[2:end].-1
    cutval = F.cutvalue
    return S, cutval
end

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
# -----------------------------------------------------------------
# scitype function (generic)

"""
    scitype(X)

The scientific type (interpretation) of `X`, as distinct from its
machine type, as specified by the active convention.

### Examples from the MLJ convention

```
julia> using ScientificTypes # or `using MLJ`
julia> scitype(3.14)
Continuous

julia> scitype([1, 2, 3, missing])
AbstractArray{Union{Missing, Count},1}

julia> scitype((5, "beige"))
Tuple{Count, Textual}

julia> using CategoricalArrays
julia> X = (gender = categorical(['M', 'M', 'F', 'M', 'F']),
            ndevices = [1, 3, 2, 3, 2])
julia> scitype(X)
Table{Union{AbstractArray{Count,1}, AbstractArray{Multiclass{2},1}}}
```

The specific behavior of `scitype` is governed by the active
convention, as returned by `ScientificTypesBase.convention()`. The
[ScientificTypes.jl
documentation](https://alan-turing-institute.github.io/ScientificTypes.jl/dev/)
details the convention demonstrated above.

"""
scitype(X;    kw...) = scitype(X, convention();     kw...)
scitype(X, C; kw...) = scitype(X, C, Val(trait(X)); kw...)

scitype(X, C, ::Val{:other}; kw...) = Unknown
scitype(::Missing;           kw...) = Missing
scitype(::Nothing;           kw...) = Nothing

scitype(t::Tuple, ::Convention; kw...) = Tuple{scitype.(t; kw...)...}

# -----------------------------------------------------------------
# convenience methods for scitype over unions

"""
    scitype_union(A)

Return the type union, over all elements `x` generated by the iterable `A`,
of `scitype(x)`. See also [`scitype`](@ref).
"""
function scitype_union(A)
    isempty(A) && return scitype(eltype(A))
    reduce((a,b)->Union{a,b}, (scitype(el) for el in A))
end

# -----------------------------------------------------------------
# Scitype for arrays

"""
    Scitype(T, C)

Method for implementers of a convention `C` to enable speed-up of
scitype evaluations for arrays.

In general, one cannot infer the scitype of an object of type
`AbstractArray{T, N}` from the machine type `T` alone.

Nevertheless, for some *restricted* machine types `U`, the statement
`type(X) == AbstractArray{T, N}` for some `T<:U` already allows one
deduce that `scitype(X) = AbstractArray{S,N}`, where `S` is determined
by `U` alone. This is the case in the *MLJ* convention, for example,
if `U = Integer`, in which case `S = Count`.

Such shortcuts are specified as follows:

```
Scitype(::Type{<:U}, ::C) = S
```

which incurs a considerable speed-up in the computation of `scitype(X)`.
There is also a speed-up for the case that `T <: Union{U, Missing}`.

For example, in the *MLJ* convention, one has

```
Scitype(::Type{<:Integer}, ::MLJ) = Count
```
"""
function Scitype end

Scitype(::Type, ::Convention) = Unknown

# to distinguish between Any type and Union{T,Missing} for some more
# specialised `T`, we define the Any case explicitly
Scitype(::Type{Any}, ::Convention) = Unknown

# for the case Union{T,Missing} we return Union{S,Missing} with S
# the scientific type corresponding to T
Scitype(::Type{Union{T,Missing}}, C::Convention) where T =
    Union{Missing,Scitype(T, C)}

# for the case Missing, we return Missing
Scitype(::Type{Missing}, C::Convention) = Missing

Scitype(::Type{Nothing}, C::Convention) = Nothing

# Broadcasting over arrays

scitype(A::Arr{T}, C::Convention, ::Val{:other}; kw...) where T =
    arr_scitype(A, C, Scitype(T, C); kw...)

"""
    arr_scitype(A, C, S; tight)

Return the scitype associated with an  array `A` of type `{T,N}` assuming an
explicit `Scitype` correspondance exist mapping `T` to `S`.
If `tight=true` and `T>:Missing` then the function checks whether there are
"true missing values", otherwise it constructs a "tight copy" of the array
without a `Union{Missing,S}` type.
"""
function arr_scitype(A::Arr{T,N}, C::Convention, S::Type;
                     tight::Bool=false) where {T,N}
    # no explicit scitype available
    S === Unknown && return Arr{scitype_union(A),N}
    # otherwise return `Arr{S,N}` or `Arr{Union{Missing,S},N}`
    if T >: Missing
        if tight
            has_missings = findfirst(ismissing, A) !== nothing
            !has_missings && return Arr{nonmissing(S),N}
        end
        return Arr{Union{S,Missing},N}
    end
    return Arr{S,N}
end

"""
    elscitype(A)

Return the element scientific type of an abstract array `A`. By definition, if
`scitype(A) = AbstractArray{S,N}`, then `elscitype(A) = S`.
"""
elscitype(X::Arr; kw...) = scitype(X; kw...) |> _get_elst

_get_elst(st::Type{Arr{T,N}}) where {T,N} = T


text:
using Documenter
using AnyMOD

makedocs(sitename="AnyMOD.jl",
    authors = "Leonard Goeke",
    pages = [
        "Introduction" => "index.md",
        "Model Elements" => Any["Sets and Mappings" => "sets.md",
                                "Parts" => "parts.md",
                                "Parameter" => Any["Parameter overview" => "parameter_overview.md","Parameter list" => "parameter_list.md"],
                                "Variables" => "variables.md",
                                "Constraints" => "constraints.md"],
        "Tools" => Any["Reporting" => Any["Plots" => "plots.md","Data files" => "data.md"],
                        "Performance and stability" => "performance.md"],
        "Annex" => Any["Functions and types" => "functions.md",
                        "Common mistakes" => "common_mistakes.md"]
        ],
    )

deploydocs(
    repo = "github.com/leonardgoeke/AnyMOD.jl.git",
)