text:
# `target` will typically a `Function` or `Type` but could also be a function-like object
struct Patch{T}
    target::T
    alternate::Function
end

macro patch(expr::Expr)
    def = splitdef(expr)

    if haskey(def, :name) && haskey(def, :body)
        target = def[:name]
    elseif !haskey(def, :name)
        throw(ArgumentError("Function definition must be a named function"))
    else
        throw(ArgumentError("Function definition must not be an empty function"))
    end

    # Include the target function name in the patch to make stack traces easier to read.
    # If the provided target uses a fully-qualified reference we'll just extract the name
    # of the function (e.g `Base.foo` -> `foo`).
    target_name = if Meta.isexpr(target, :.)
        target.args[2].value
    else
        target
    end

    def[:name] = gensym(string(target_name, "_patch"))
    alternate = combinedef(def)

    # We need to evaluate the alternate function in the context of the `@patch` macro in
    # order to support closures.
    return esc(:(Mocking.Patch($target, $alternate)))
end

struct PatchEnv
    mapping::Dict{Any,Vector{Function}}
    debug::Bool
end

function PatchEnv(patches, debug::Bool=false)
    pe = PatchEnv(debug)
    apply!(pe, patches)
    return pe
end

PatchEnv(debug::Bool=false) = PatchEnv(Dict{Any,Vector{Function}}(), debug)

function Base.:(==)(pe1::PatchEnv, pe2::PatchEnv)
    return pe1.mapping == pe2.mapping && pe1.debug == pe2.debug
end

"""
    merge(pe1::PatchEnv, pe2::PatchEnv) -> PatchEnv

Merge the two `PatchEnv` instances.

This is done in such a way that the following always holds:

```
patches_1 = Patch[...]
patches_2 = Patch[...]
patches = vcat(patches_1, patches_2)

pe1 = PatchEnv(patches_1)
pe2 = PatchEnv(patches_2)
pe = PatchEnv(patches)

@assert pe == merge(pe1, pe2)
```

The `debug` flag will be set to true if either `pe1` or `pe2` have it set to true.
"""
function Base.merge(pe1::PatchEnv, pe2::PatchEnv)
    mapping = mergewith(vcat, pe1.mapping, pe2.mapping)
    return PatchEnv(mapping, pe1.debug || pe2.debug)
end

function apply!(pe::PatchEnv, p::Patch)
    alternate_funcs = get!(Vector{Function}, pe.mapping, p.target)
    push!(alternate_funcs, p.alternate)
    return pe
end

function apply!(pe::PatchEnv, patches)
    for p in patches
        apply!(pe, p)
    end
    return pe
end

"""
    apply(body::Function, patches; debug::Bool=false)
    apply(body::Function, pe::PatchEnv)

Convenience function to run `body` in the context of the given `patches`.

This is intended to be used with do-block notation, e.g.:

```
patch = @patch ...
apply(patch) do
    ...
end
```

## Nesting

Note that calls to apply will nest the patches that are applied. If multiple patches
are made to the same method, the innermost patch takes precedence.

The following two examples are equivalent:

```
patch_2 = @patch ...
apply([patch, patch_2]) do
    ...
end
```

```
apply(patch) do
    apply(patch_2) do
        ...
    end
end
```
"""
function apply(body::Function, pe::PatchEnv)
    prev_pe = get_active_env()
    set_active_env(merge(prev_pe, pe))
    try
        return body()
    finally
        set_active_env(prev_pe)
    end
end

function apply(body::Function, patches; debug::Bool=false)
    return apply(body, PatchEnv(patches, debug))
end

const PATCH_ENV = Ref{PatchEnv}(PatchEnv())
set_active_env(pe::PatchEnv) = (PATCH_ENV[] = pe)
get_active_env() = PATCH_ENV[]


text:
###############################################################################
#
# # ListofStates tutorial
#
# We illustrate here the use of ListofStates in dealing with a warm start
# procedure.
#
# ListofStates can also prove the user history over the iteration process.
#
# We compare the resolution of a convex unconstrained problem with 3 variants:
# - a steepest descent method
# - an inverse-BFGS method
# - a mix with 5 steps of steepest descent and then switching to BFGS with
#the history (using the strength of the ListofStates).
#
###############################################################################
using Stopping, NLPModels, LinearAlgebra, Test, Printf

import Stopping.armijo
function armijo(xk, dk, fk, slope, f)
  t = 1.0
  fk_new = f(xk + dk)
  while f(xk + t * dk) > fk + 1.0e-4 * t * slope
    t /= 1.5
    fk_new = f(xk + t * dk)
  end
  return t, fk_new
end

#Newton's method for optimization:
function steepest_descent(stp :: NLPStopping)

  xk = stp.current_state.x
  fk, gk = objgrad(stp.pb, xk)

  OK = update_and_start!(stp, fx = fk, gx = gk)

  @printf "%2s %9s %7s %7s %7s\n" "k" "fk" "||∇f(x)||" "t" "λ"
  @printf "%2d %7.1e %7.1e\n" stp.meta.nb_of_stop fk norm(stp.current_state.current_score)
  while !OK
    dk = - gk
    slope = dot(dk, gk)
    t, fk = armijo(xk, dk, fk, slope, x->obj(stp.pb, x))
    xk += t * dk
    fk, gk = objgrad(stp.pb, xk)

    OK = update_and_stop!(stp, x = xk, fx = fk, gx = gk)

    @printf "%2d %9.2e %7.1e %7.1e %7.1e\n" stp.meta.nb_of_stop fk norm(stp.current_state.current_score) t slope
  end
  return stp
end

function bfgs_quasi_newton_armijo(stp :: NLPStopping; Hk = nothing)

  xk = stp.current_state.x
  fk, gk = objgrad(stp.pb, xk)
  gm = gk

  dk, t = similar(gk), 1.
  if isnothing(Hk)
    Hk = I #start from identity matrix
  end

  OK = update_and_start!(stp, fx = fk, gx = gk)

  @printf "%2s %7s %7s %7s %7s\n" "k" "fk" "||∇f(x)||" "t" "cos"
  @printf "%2d %7.1e %7.1e\n" stp.meta.nb_of_stop fk norm(stp.current_state.current_score)

  while !OK
    if stp.meta.nb_of_stop != 0
      sk = t * dk
      yk = gk - gm
      ρk = 1/dot(yk, sk)
      #we need yk'*sk > 0 for instance yk'*sk ≥ 1.0e-2 * sk' * Hk * sk
      Hk = ρk ≤ 0.0 ? Hk : (I - ρk * sk * yk') * Hk * (I - ρk * yk * sk') + ρk * sk * sk'
      if norm(sk) ≤ 1e-14
        break
      end
      #H2 = Hk + sk * sk' * (dot(sk,yk) + yk'*Hk*yk )*ρk^2 - ρk*(Hk * yk * sk' + sk * yk'*Hk)
    end
    dk = - Hk * gk
    slope = dot(dk, gk) # ≤ -1.0e-4 * norm(dk) * gnorm
    t, fk = armijo(xk, dk, fk, slope, x->obj(stp.pb, x))

    xk = xk + t * dk
    gm = copy(gk)
    gk = grad(stp.pb, xk)

    OK = update_and_stop!(stp, x = xk, fx = fk, gx = gk)
    @printf "%2d %7.1e %7.1e %7.1e %7.1e\n" stp.meta.nb_of_stop fk norm(stp.current_state.current_score) t slope
  end
  stp.stopping_user_struct = Dict( :Hk => Hk)
  return stp
end
using Test

############ PROBLEM TEST #############################################
fH(x) = (x[2]+x[1].^2-11).^2+(x[1]+x[2].^2-7).^2
nlp = ADNLPModel(fH, [10., 20.])
stp = NLPStopping(nlp, optimality_check = unconstrained_check,
                 atol = 1e-6, rtol = 0.0, max_iter = 100)

reinit!(stp, rstate = true, x = nlp.meta.x0)
steepest_descent(stp)

@test status(stp) == :Optimal
@test stp.listofstates == VoidListofStates()

@show elapsed_time(stp)
@show nlp.counters

reinit!(stp, rstate = true, x = nlp.meta.x0, rcounters = true)
bfgs_quasi_newton_armijo(stp)

@test status(stp) == :Optimal
@test stp.listofstates == VoidListofStates()

@show elapsed_time(stp)
@show nlp.counters

NLPModels.reset!(nlp)
stp_warm = NLPStopping(nlp, optimality_check = unconstrained_check,
                      atol = 1e-6, rtol = 0.0, max_iter = 5,
                      n_listofstates = 5) #shortcut for list = ListofStates(5, Val{NLPAtX{Float64,Array{Float64,1},Array{Float64,2}}}()))
steepest_descent(stp_warm)
@test status(stp_warm) == :IterationLimit
@test length(stp_warm.listofstates) == 5

Hwarm = I
for i=2:5
  sk = stp_warm.listofstates.list[i][1].x - stp_warm.listofstates.list[i-1][1].x
  yk = stp_warm.listofstates.list[i][1].gx - stp_warm.listofstates.list[i-1][1].gx
  ρk = 1/dot(yk, sk)
  if ρk > 0.0
    global Hwarm = (I - ρk * sk * yk') * Hwarm * (I - ρk * yk * sk') + ρk * sk * sk'
  end
end

reinit!(stp_warm)
stp_warm.meta.max_iter = 100
bfgs_quasi_newton_armijo(stp_warm, Hk = Hwarm)
status(stp_warm)

@show elapsed_time(stp_warm)
@show nlp.counters

#=
 k        fk ||∇f(x)||       t       λ
 0 1.7e+05 3.2e+04
 1  2.73e+04 8.6e+03 1.0e-03 -1.1e+09
 2  1.80e+03 1.1e+03 2.3e-03 -7.3e+07
 3  1.24e+03 7.9e+02 1.2e-02 -1.3e+06
 4  6.37e+01 2.4e+01 1.2e-02 -6.3e+05
 5  1.34e+01 5.8e+01 2.0e-01 -8.3e+02
 6  5.87e+00 2.5e+01 1.3e-01 -3.5e+03
 7  2.88e+00 2.4e+01 2.6e-02 -6.7e+02
 8  2.42e+00 1.8e+01 1.7e-02 -6.1e+02
 9  6.58e-01 1.2e+01 1.2e-02 -6.1e+02
10  1.64e-01 5.3e+00 1.2e-02 -1.7e+02
11  4.96e-02 3.2e+00 1.2e-02 -4.4e+01
12  1.44e-02 1.6e+00 1.2e-02 -1.3e+01
13  4.35e-03 9.2e-01 1.2e-02 -3.9e+00
14  1.29e-03 5.0e-01 1.2e-02 -1.2e+00
15  3.87e-04 2.7e-01 1.2e-02 -3.5e-01
16  1.15e-04 1.5e-01 1.2e-02 -1.0e-01
17  3.45e-05 8.2e-02 1.2e-02 -3.1e-02
18  1.03e-05 4.5e-02 1.2e-02 -9.2e-03
19  3.08e-06 2.4e-02 1.2e-02 -2.8e-03
20  9.21e-07 1.3e-02 1.2e-02 -8.2e-04
21  2.75e-07 7.3e-03 1.2e-02 -2.5e-04
22  8.23e-08 4.0e-03 1.2e-02 -7.4e-05
23  2.46e-08 2.2e-03 1.2e-02 -2.2e-05
24  7.35e-09 1.2e-03 1.2e-02 -6.6e-06
25  2.20e-09 6.5e-04 1.2e-02 -2.0e-06
26  6.57e-10 3.6e-04 1.2e-02 -5.9e-07
27  1.96e-10 1.9e-04 1.2e-02 -1.8e-07
28  5.87e-11 1.1e-04 1.2e-02 -5.3e-08
29  1.75e-11 5.8e-05 1.2e-02 -1.6e-08
30  5.24e-12 3.2e-05 1.2e-02 -4.7e-09
31  1.57e-12 1.7e-05 1.2e-02 -1.4e-09
32  4.68e-13 9.5e-06 1.2e-02 -4.2e-10
33  1.40e-13 5.2e-06 1.2e-02 -1.3e-10
34  4.18e-14 2.8e-06 1.2e-02 -3.7e-11
35  1.25e-14 1.6e-06 1.2e-02 -1.1e-11
36  3.74e-15 8.5e-07 1.2e-02 -3.3e-12
elapsed_time(stp) = 0.7508440017700195
nlp.counters =   Counters:
             obj: ████████████████████ 889               grad: █⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 37                cons: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
            jcon: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                jgrad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                  jac: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
           jprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0               jtprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                 hess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
           hprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0               jhprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0

 k      fk ||∇f(x)||       t     cos
 0 1.7e+05 3.2e+04
 1 2.7e+04 8.6e+03 1.0e-03 -1.1e+09
 2 1.8e+04 4.5e+03 1.2e-02 -1.8e+06
 3 2.5e+03 1.3e+03 1.0e+00 -7.1e+04
 4 1.2e+03 8.5e+02 1.0e+00 -1.7e+03
 5 3.2e+02 3.3e+02 1.0e+00 -1.4e+03
 6 9.8e+01 1.4e+02 1.0e+00 -3.2e+02
 7 2.7e+01 6.0e+01 1.0e+00 -1.1e+02
 8 6.4e+00 2.4e+01 1.0e+00 -3.0e+01
 9 9.9e-01 7.9e+00 1.0e+00 -8.2e+00
10 6.3e-02 1.9e+00 1.0e+00 -1.5e+00
11 8.7e-04 3.2e-01 1.0e+00 -1.1e-01
12 3.6e-05 7.9e-02 1.0e+00 -1.6e-03
13 1.4e-05 4.2e-02 1.0e+00 -2.9e-05
14 2.0e-07 3.4e-03 1.0e+00 -2.6e-05
15 4.1e-09 4.9e-04 1.0e+00 -3.6e-07
16 2.9e-12 2.5e-05 1.0e+00 -8.1e-09
17 2.5e-15 6.3e-07 1.0e+00 -5.6e-12
elapsed_time(stp) = 0.017869949340820312
nlp.counters =   Counters:
             obj: ████████████████████ 91                grad: ████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 18                cons: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
            jcon: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                jgrad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                  jac: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
           jprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0               jtprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                 hess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
           hprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0               jhprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0

 k        fk ||∇f(x)||       t       λ
 0 1.7e+05 3.2e+04
 1  2.73e+04 8.6e+03 1.0e-03 -1.1e+09
 2  1.80e+03 1.1e+03 2.3e-03 -7.3e+07
 3  1.24e+03 7.9e+02 1.2e-02 -1.3e+06
 4  6.37e+01 2.4e+01 1.2e-02 -6.3e+05
 5  1.34e+01 5.8e+01 2.0e-01 -8.3e+02
 6  5.87e+00 2.5e+01 1.3e-01 -3.5e+03
 k      fk ||∇f(x)||       t     cos
 0 5.9e+00 2.5e+01
 1 3.8e+00 2.7e+01 1.7e-02 -1.1e+03
 2 2.8e+00 2.4e+01 4.4e-01 -1.1e+01
 3 1.4e+00 1.2e+01 3.0e-01 -3.0e+01
 4 1.1e-02 1.3e+00 1.0e+00 -2.5e+00
 5 9.0e-05 9.2e-02 1.0e+00 -2.5e-02
 6 7.9e-08 3.9e-03 1.0e+00 -1.8e-04
 7 7.7e-10 3.8e-04 1.0e+00 -1.4e-07
 8 1.3e-19 4.2e-09 1.0e+00 -1.5e-09
elapsed_time(stp_warm) = 0.01520395278930664
nlp.counters =   Counters:
             obj: ████████████████████ 192               grad: ██⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 16                cons: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
            jcon: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                jgrad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                  jac: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
           jprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0               jtprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                 hess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
           hprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0               jhprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0

  Counters:
             obj: ████████████████████ 192               grad: ██⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 16                cons: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
            jcon: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                jgrad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                  jac: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
           jprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0               jtprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                 hess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
           hprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0               jhprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0

=#

ext: jl


text:
## Cartpole
import TrajectoryOptimization
traj_folder = joinpath(dirname(pathof(TrajectoryOptimization)),"..")
urdf_folder = joinpath(traj_folder, "dynamics/urdf")
urdf_cartpole = joinpath(urdf_folder, "cartpole.urdf")

cartpole_urdf = Model(urdf_cartpole,[1.;0.]) # underactuated, only control of slider

function cartpole_dynamics!(ẋ::AbstractVector{T}, x::AbstractVector{T}, u::AbstractVector{T}) where T
    mc = 1.0  # mass of the cart in kg (10)
    mp = 0.2    # mass of the pole (point mass at the end) in kg
    l = 0.5   # length of the pole in m
    g = 9.81  # gravity m/s^2

    q = x[1:2]
    qd = x[3:4]

    if isfinite(q[2])
        s = sin(q[2])
        c = cos(q[2])
    else
        s = Inf
        c = Inf
    end

    H = [mc+mp mp*l*c; mp*l*c mp*l^2]
    C = [0 -mp*qd[2]*l*s; 0 0]
    G = [0; mp*g*l*s]
    B = [1; 0]

    qdd = -H\(C*qd + G - B*u[1])

    ẋ[1:2] = qd
    ẋ[3:4] = qdd
    return nothing
end

n,m = 4,1
cartpole = Model(cartpole_dynamics!,n,m)

# unknown friction
function cartpole_dynamics_uncertain!(ẋ, x, u, w)
    mc = 1.0  # mass of the cart in kg (10)
    mp = 0.2    # mass of the pole (point mass at the end) in kg
    l = 0.5   # length of the pole in m
    g = 9.81  # gravity m/s^2
    mu = 0.1

    q = x[1:2]
    qd = x[3:4]

    if isfinite(q[2])
        s = sin(q[2])
        c = cos(q[2])
    else
        s = Inf
        c = Inf
    end

    H = [mc+mp mp*l*c; mp*l*c mp*l^2]
    C = [0 -mp*qd[2]*l*s; 0 0]
    G = [0; mp*g*l*s]
    B = [1; 0]

    qdd = -H\(C*qd + G - B*u[1] - B*w[1])

    ẋ[1:2] = qd
    ẋ[3:4] = qdd

    return nothing
end

n = 4; m = 1; r = 1
cartpole_uncertain = UncertainModel(cartpole_dynamics_uncertain!,n,m,r)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
###################################################
# Script to train a simple perceptron

###################################################
# Include libraries and core

using Pkg; Pkg.activate("./../core/"); Pkg.instantiate()
include("./../core/src/ratenets/core.jl")

###################################################
# Define task and network

dataset = "MNIST" # dataset to be used: "MNIST" or "CIFAR10"
hidden_sizes = [0] # number of neurons in hidden layer (0 for simple perceptron)
n_inits = 3 # number of initializations per net (for averaging)
iterations = 10^6 # iterations per initialization
learningrates = [1e-4] # learningrates per layer
nonlinearity = [relu!] #nonlinearities
nonlinearity_diff = [relu_diff!]

###################################################
# Load and preprocess data

smallimgs, labels, smallimgstest, labelstest = import_data(dataset)
subtractmean!(smallimgs) # substract pixel-wise mean
subtractmean!(smallimgstest)

###################################################
# Train network

print("train simple perceptron (without hidden layer)...\n")
for j in 1:n_inits
  net = Network([size(smallimgs)[1], 10])
  learn!(net, getsmallimg, getlabel, iterations, learningrates, 0.,
                     nonlinearity = nonlinearity, nonlinearity_diff = nonlinearity_diff)
  error_train = geterrors(net, smallimgs, labels, nonlinearity = nonlinearity)
  error_test = geterrors(net, smallimgstest, labelstest, nonlinearity = nonlinearity)
  print(string(100 * error_train, " % on training set\n"))
  print(string(100 * error_test, " % on test set\n"))
end


text:
using Revise
using Test: @testset, @test

include("Interface.jl")
include("ExampleFunctions.jl")

using .Interface: @ODE
using .StrangeAttractors: yuwang, TSUCS2, lorenz, TSUCS1

u = [[1., 2., 3.] [4., 5., 6.] [7., 8., 9.]]

lorenz(u)

@testset "Interface" begin
    @testset "Three-Dimensional Functions" begin

        @testset "Typing" begin
            u = [[1., 2., 3.] [4., 5., 6.] [7., 8., 9.]]

            @test lorenz(u) ≈ [[10., 10., 10.] [23., 83., 125.] [-6., 4., 32.]]
        end

        @testset "Time Functionality" begin
            u = [[1., 2., 3.] [4., 5., 6.] [7., 8., 9.]]

            time_ODE = Interface.@ODE quote
                dx = t
                dy = x + y
                dz = y + z
            end

            @test length(methods(time_ODE)) == 1

            @test time_ODE(0., u)[:,1] ≈ Vector{Float64}([0, 0, 0])
            @test time_ODE(2., u)[:,1] ≈ Vector{Float64}([2, 2, 2])
        end
    end

    @testset "Two-Dimensional Functions" begin

    end
end

lorenz([[1.,2.,3.] [3.,4.,5.] [6.,7., 8.]], σ = 10, ρ = 28, β = 8/3)
TSUCS2(0., [[1,2,3] [3,4,5] [6,7, 8]], α = 1, δ = 2, ς = 4, ζ = 8, β = 5, ϵ = 6)
yuwang(0., [[1,2,3] [3,4,5] [6,7, 8]], α = 1, β = 2, σ = 4, δ = 8)


text:
# Autogenerated wrapper script for HelloWorldC_jll for armv7l-linux-gnueabihf
export hello_world

JLLWrappers.@declare_executable_product(hello_world)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        hello_world,
        "bin/hello_world",
    )

end  # __init__()


text:
## Test nonlinear function
function nonlinear_function(y, u, v)
    return ((0.8 - 0.5*exp(-y[1]^2))*y[1] - (0.3 + 0.9*exp(-y[1]^2))*y[2] + u[1]
            + 0.2*u[2] + 0.1*u[1]*u[2] + 0.1*v[1] + 0.3*v[2])
end

@testset "Test simulation for arbitrary nonlinear function" begin
    u = NarmaxLasso.generate_random_input(10, 5)
    y = NarmaxLasso.simulate(u, nonlinear_function, [1, 2], [1, 2], [1, 2], σv=0.1)
    @test y ≈ [0.0, 0.0, 0.318563, 0.480975, 0.134545, -0.139843, 0.458256,
                   0.603922, 0.422271, -0.0208464] rtol=1e-2
    @test length(y) == length(u)
end

## Test one step ahead
m1 = NarmaxLasso.Monomial([1, 0, 0])
m2 = NarmaxLasso.Monomial([0, 0, 1])
m3 = NarmaxLasso.Monomial([0, 2, 0])

@testset "Test one-step-ahead" begin
    @test NarmaxLasso.one_step_ahead([1.0], [2.0], [3.0],
                                     [m1, m2, m3],
                                     [1.0, 2.0, 3.0]) ≈ 19.0
    @test NarmaxLasso.one_step_ahead([1.0], [2.0, 3.0],
                                     [m1, m2, m3],
                                     [1.0, 2.0, 3.0]) ≈ 19.0
end

# simulate
@testset "Test simulate" begin
    mdl1 = NarmaxLasso.NarmaxRegressors([m1, m2, m3], [1], [1], [1])
    mdl2 = NarmaxLasso.NarmaxRegressors([m1, m2, m3], [1], [1, 2])
    u = NarmaxLasso.generate_random_input(20, 20)
    y = NarmaxLasso.simulate(u, mdl1, [0.5, 1, 0])
    @test (1-0.5)*y[end]/(u[end]^2) ≈ 1.0 rtol=1e-2

    y = NarmaxLasso.simulate(u, mdl2, [0.5, 1, 0.2])
    @test (1-0.5)*y[end]/(0.2*u[end]^2+u[end]) ≈ 1.0 rtol=1e-2
end


text:
@testset "Norms" begin
    x = [2.0im, 3.0 - 1im, 5.0 + 2.0im]
    y = [-2.0im, 3.0 - 1im, 5.0 + 2.0im]
    Norm = HC.InfNorm()
    @test Norm isa HC.AbstractNorm
    @test Norm(x) ≈ abs(5.0 + 2.0im)
    @test Norm(x, x) ≈ 0.0
    @test Norm(x, y) ≈ abs(4im)

    WeightedNorm = HC.WeightedNorm(Norm, x)
    @test HC.weights(WeightedNorm) == ones(3)
    WeightedNorm .= 2.0
    @test HC.weights(WeightedNorm) == [2.0, 2, 2]
    WeightedNorm[1] = 4.0
    @test WeightedNorm isa HC.AbstractNorm
    @test WeightedNorm(x) ≈ 0.5 * abs(5 + 2.0im)
    @test WeightedNorm(x, x) ≈ 0.0
    @test WeightedNorm(x, y) ≈ 0.25 * abs(4im)

    Norm = HC.EuclideanNorm()
    @test Norm isa HC.AbstractNorm

    @test Norm(x) ≈ LinearAlgebra.norm(x, 2)
    @test Norm(x, x) ≈ 0.0
    @test Norm(x, y) ≈ abs(4im)

    WeightedNorm = HC.WeightedNorm([4.0, 2.0, 2.0], Norm)
    @test WeightedNorm isa HC.AbstractNorm
    @test WeightedNorm(x) ≈ sqrt(abs2(x[1]) / 16 + abs2(x[2]) / 4 + abs2(x[3]) / 4)
    @test WeightedNorm(x, x) ≈ 0.0
    @test WeightedNorm(x, y) ≈ sqrt(abs2(4im) / 16)
end


text:
export BFFPSV18

# ===============================================================
# Bogomolov, Forets, Frehse, Podelski, Schilling, Viry. HSCC 2018
# ===============================================================

# dummy functions for option :lazy_inputs_interval
lazy_inputs_interval_always = (k -> true)
lazy_inputs_interval_never = (k -> false)

function ispartition(partition::AbstractVector{<:AbstractVector{Int}})
    current = 1
    for block in partition
        for i in block
            if i != current
                return false
            end
            current += 1
        end
    end
    return true
end

function options_BFFPSV18()
    return OptionSpec[
        # general options
        OptionSpec(:discretization, "forward", domain=String,
            info="model for bloating/continuous time analysis"),
        OptionSpec(:algorithm, "explicit", domain=String, domain_check=(
            v  ->  v in ["explicit", "wrap"]), info="algorithm backend"),
        OptionSpec(:δ, 1e-2, domain=Float64, aliases=[:sampling_time],
            domain_check=(v  ->  v > 0.), info="time step"),
        OptionSpec(:vars, Int[], domain=AbstractVector{Int}, domain_check=(
            v  ->  length(v) > 0 && all(e -> e > 0, v)),
            info="variables of interest; default: all variables"),
        OptionSpec(:partition, [Int[]],
            domain=AbstractVector{<:AbstractVector{Int}}, domain_check=
            ispartition,
            info="block partition; a block is represented by a vector " *
                 "containing its indices"),

        # discretization options
        OptionSpec(:sih_method, "concrete", domain=String,
            info="method to compute the symmetric interval hull in discretization"),
        OptionSpec(:exp_method, "base", domain=String,
            info="method to compute the matrix exponential"),
        OptionSpec(:assume_sparse, false, domain=Bool,
            info="use an analysis for sparse discretized matrices?"),

        # reachability options
        OptionSpec(:lazy_inputs_interval, lazy_inputs_interval_always,
            domain=Union{Int, Function},
            domain_check=(v  ->  !(v isa Int) || v >= -1),
            info="length of interval in which the inputs are handled as a " *
                 "lazy set (``-1`` for 'never'); may generally also be a " *
                 "predicate over indices; the default corresponds to ``-1``"),

        # approximation options
        OptionSpec(:block_options, nothing, domain=Any,
            info="short hand to set ':block_options_init' and " *
                 "':block_options_iter'"),
        OptionSpec(:block_options_init, nothing, domain=Any,
            info="option for the approximation of the initial states " *
                 "(during decomposition)"),
        OptionSpec(:block_options_iter, nothing, domain=Any,
            info="option for the approximation of the states ``X_k``, ``k>0``"),

        # convenience options
        OptionSpec(:assume_homogeneous, false, domain=Bool,
            info="ignore dynamic inputs during the analysis?"),
        OptionSpec(:eager_checking, true, domain=Bool,
            info="terminate as soon as property violation was detected?"),
    ]
end

function normalization_BFFPSV18!(𝑂::TwoLayerOptions)
    # :lazy_inputs_interval option: convert integers to functions
    if haskey_specified(𝑂, :lazy_inputs_interval)
        v = 𝑂[:lazy_inputs_interval]
        if v isa Int
            if v == -1
                𝑂.specified[:lazy_inputs_interval] = lazy_inputs_interval_never
            elseif v == 0
                𝑂.specified[:lazy_inputs_interval] = lazy_inputs_interval_always
            else
                𝑂.specified[:lazy_inputs_interval] = (k -> k % v == 0)
            end
        end
    end

    # :block_options options: use convenience option for '_init' and '_iter'
    if !haskey_specified(𝑂, :block_options_init) &&
            haskey_specified(𝑂, :block_options)
        𝑂.specified[:block_options_init] = 𝑂[:block_options]
    end
    if !haskey_specified(𝑂, :block_options_iter) &&
            haskey_specified(𝑂, :block_options)
        𝑂.specified[:block_options_iter] = 𝑂[:block_options]
    end

    nothing
end

function validation_BFFPSV18(𝑂)
    # block_options
    for b_options in [:block_options, :block_options_init, :block_options_iter]
        if haskey_specified(𝑂, b_options)
            bo = 𝑂[b_options]
            if bo isa Type{<:LazySet} || bo isa Type{<:AbstractDirections}
                # uniform options
            elseif bo isa AbstractVector || bo isa Dict{Int, Any}
                # mapping
            elseif bo isa Real || bo isa Pair{<:UnionAll, <:Real}
                ε = bo isa Real ? bo : bo[2]
                if ε <= 0
                    throw(DomainError(ε, "the `$b_options` option must be " *
                                         "positive"))
                end
            elseif b_options == :block_options_iter && bo == nothing
                # no overapproximation
            else
                throw(DomainError(bo == nothing ? "nothing" : bo,
                    "the `$b_options` option does not accept the given input"))
            end
        end
    end

    nothing
end

"""
    BFFPSV18 <: ContinuousPost

Implementation of the reachability algorithm using decompositions from [1].

This algorithm applies to continuous affine ordinary differential equations of
the form

```
    x'(t) = Ax(t) + u(t),\\qquad x(0) ∈ X0, u(t) ∈ U ∀ t ∈ [0, T],    (1)
```
where `T` is the (finite) time horizon. The result of the algorithm is a
sequence of sets `Xk` whose union overapproximates the exact flowpipe of (1).

### Fields

- `options` -- an `Options` structure that holds the algorithm-specific options

The following options are available:

```julia
$(print_option_spec(options_BFFPSV18()))
```

See the `Notes` below for examples and an explanation of the options.

### Algorithm

We refer to [1] for technical details.

### Notes

The compositional approach offers a tradeoff between accuracy and performance,
since one can choose different partitions and different set representations.

These notes illustrate the options that can be passed to the solver. We begin
by explaining the default options used by the algorithm. Then we show how to
specify the computation of a subset of the flowpipe, which helps to improve
performance, specially for large systems.

We continue to show some lower-level options such as how to specify a custom
block decomposition (the set types used and the sizes of the blocks).

We conclude showing that there are two different *modes* of operation:
flowpipe computation and safety property checking, and examples of use.

#### Running example and default options

By default, `BFFPSV18` uses a uniform partition of dimension *one* over all
variables.

For example, consider the two-dimensional example from [2]:

```jldoctest BFFPSV18_example
julia> using Reachability, MathematicalSystems

julia> A = [-1. -4   0  0   0;
             4  -1   0  0   0;
             0   0  -3  1   0;
             0   0  -1  3   0;
             0   0   0  0  -2];

julia> X₀ = BallInf(ones(5), 0.1);

julia> problem = InitialValueProblem(LinearContinuousSystem(A), X₀);

julia> sol = solve(problem, Options(:T=>5.0), op=BFFPSV18(Options(:δ=>0.04)));
```
Let's check that the dimension of the first set computed is 5:

```jldoctest BFFPSV18_example
julia> X₁ = set(first(first(sol.flowpipes).reachsets));

julia> dim(X₁)
5

julia> sol.options[:partition]
5-element Array{Array{Int64,1},1}:
 [1]
 [2]
 [3]
 [4]
 [5]
```
Here, `[1]` corresponds to a block of size one for variable `5`, etc., until
variable `5`.
Let's check the set type used is `Interval`:

```jldoctest BFFPSV18_example
julia> all(Xi isa Interval for Xi in array(X₁))
true
```

#### Specifying a subset of all variables

Suppose that we are only interested in variables `1` and `2`. Then we can
specify it using the `:vars` option:

```jldoctest BFFPSV18_example
julia> sol = solve(problem, Options(:T=>5.0),
                   op=BFFPSV18(Options(:δ=>0.04, :vars=>[1, 2])));

julia> X₁ = set(first(first(sol.flowpipes).reachsets));

julia> dim(X₁)
2
```
The sets used are `Interval`s:

```jldoctest BFFPSV18_example
julia> all(Xi isa Interval for Xi in array(X₁))
true
```

#### Specifying the set type

We can as well specify using other set types for the overapproximation of the
flowpipe, such as hyperrectangles. It can be specified using the
`:overapproximation` option:

```jldoctest BFFPSV18_example
julia> sol = solve(problem, Options(:T=>5.0),
                   op=BFFPSV18(Options(:δ=>0.04, :vars=>[1, 2],
                                       :block_options=>Hyperrectangle)));

julia> X₁ = set(first(first(sol.flowpipes).reachsets));

julia> all(Xi isa Hyperrectangle for Xi in array(X₁))
true
```

To increase the accuracy, the `overapproximation` option offers other
possibilities.
For example, one can use template directions instead of a fixed set type.
Below we use octagonal directions:

```jldoctest BFFPSV18_example
julia> sol = solve(problem, Options(:T=>5.0),
                   op=BFFPSV18(Options(:δ=>0.04, :vars=>[1, 2],
                                       :block_options=>OctDirections)));

julia> X₁ = set(first(first(sol.flowpipes).reachsets));

julia> all(Xi isa HPolytope for Xi in array(X₁))
true
```

### References

[1] [Reach Set Approximation through Decomposition with Low-dimensional Sets
and High-dimensional Matrices](https://dl.acm.org/citation.cfm?id=3178128).
S. Bogomolov, M. Forets, G. Frehse, A. Podelski, C. Schilling, F. Viry.
HSCC '18 Proceedings of the 21st International Conference on Hybrid Systems:
Computation and Control (part of CPS Week).

[2] Althoff, Matthias. Reachability analysis and its application to the safety
assessment of autonomous cars. Diss. Technische Universität München, 2010.
"""
struct BFFPSV18 <: ContinuousPost
    options::TwoLayerOptions

    function BFFPSV18(𝑂::Options)
        normalized_𝑂 = validate_and_wrap_options(𝑂, options_BFFPSV18();
            validation=validation_BFFPSV18,
            normalization=normalization_BFFPSV18!)
        return new(normalized_𝑂)
    end
end

# convenience constructor from pairs of symbols
BFFPSV18(𝑂::Pair{Symbol,<:Any}...) = BFFPSV18(Options(Dict{Symbol,Any}(𝑂)))

# default options
BFFPSV18() = BFFPSV18(Options())

init(𝒫::BFFPSV18, 𝑆::AbstractSystem, 𝑂::Options) = init!(𝒫, 𝑆, copy(𝑂))

function init!(𝒫::BFFPSV18, 𝑆::AbstractSystem, 𝑂::Options)
    # state dimension for (purely continuous or purely discrete systems)
    𝑂copy = copy(𝑂)
    𝑂copy[:n] = statedim(𝑆)

    # solver-specific options (adds default values for unspecified options)
    𝑂validated = validate_solver_options_and_add_default_values!(𝑂copy)

    # :vars option; default: all variables
    if haskey_specified(𝒫.options, :vars)
        𝑂validated[:vars] = 𝒫.options[:vars]
    else
        𝑂validated[:vars] = 1:𝑂validated[:n]
    end

    # :partition option: use 1D blocks
    if haskey_specified(𝒫.options, :partition)
        𝑂validated[:partition] = 𝒫.options[:partition]
    else
        𝑂validated[:partition] = [[i] for i in 1:𝑂validated[:n]]
    end

    # :blocks option (internal only)
    # list of all interesting block indices in the partition
    𝑂validated[:blocks] = compute_blocks(𝑂validated[:vars], 𝑂validated[:partition])

    # :block_options_init & :block_options_iter options:
    # set default according to :partition
    if !haskey_specified(𝒫.options, :block_options_init)
        𝑂validated[:block_options_init] =
            compute_default_block_options(𝑂validated[:partition])
    end
    if !haskey_specified(𝒫.options, :block_options_iter)
        𝑂validated[:block_options_iter] =
            compute_default_block_options(𝑂validated[:partition])
    end

    if 𝑂validated[:project_reachset]
        𝑂validated[:output_function] = nothing
    else
        𝑂validated[:output_function] = 𝑂validated[:projection_matrix]
    end

    return 𝑂validated
end

"""
    post(𝒫::BFFPSV18, 𝑆::AbstractSystem, 𝑂::Options)

Calculate the reachable states of the given initial value problem using `BFFPSV18`.

### Input

- `𝒫` -- post operator of type `BFFPSV18`
- `𝑆` -- sytem, initial value problem for a continuous ODE
- `𝑂` -- algorithm-specific options
"""
function post(𝒫::BFFPSV18, 𝑆::AbstractSystem, 𝑂_input::Options)
    𝑂 = TwoLayerOptions(merge(𝑂_input, 𝒫.options.specified), 𝒫.options.defaults)

    if 𝑂[:mode] == "reach"
        info("Reachable States Computation...")
        @timing begin
            flowpipe = reach(𝑆, 𝑂)
            info("- Total")
        end

        # Projection
        if 𝑂[:project_reachset]
            info("Projection...")
            flowpipe_proj = @timing project(flowpipe, 𝑂)
        else
            flowpipe_proj = flowpipe
        end

        return ReachSolution(flowpipe_proj, 𝑂_input)

    elseif 𝑂[:mode] == "check"
        info("invariants are currently not supported in 'check' mode")

        # Input -> Output variable mapping in property
        property = inout_map_property(𝑂[:property], 𝑂[:partition], 𝑂[:blocks], 𝑂[:n])

        # =================
        # Property checking
        # =================
        info("Property Checking...")
        @timing begin
            answer = check_property(𝑆, property, 𝑂)
            info("- Total")
        end

        if answer == 0
            info("The property is satisfied!")
            return CheckSolution(true, -1, 𝑂_input)
        else
            info("The property may be violated at index $answer," *
                " (time point $(answer * 𝑂[:δ]))!")
            return CheckSolution(false, answer, 𝑂_input)
        end
    else
        error("unsupported mode $(𝑂[:mode])")
    end # mode
end

function compute_blocks(vars, partition)
    blocks = Vector{Int}()
    sizehint!(blocks, length(vars))
    next = 0
    var_idx = 1
    for (i, block) in enumerate(partition)
        next += length(block)
        if vars[var_idx] <= next
            push!(blocks, i)
            var_idx += 1
            while var_idx <= length(vars) && vars[var_idx] <= next
                var_idx += 1
            end
            if var_idx > length(vars)
                break
            end
        end
    end
    @assert var_idx == length(vars) + 1
    sizehint!(blocks, length(blocks))
    return blocks
end

function compute_default_block_options(partition)
    # use Interval for 1D blocks and Hyperrectangle otherwise
    block_options = Vector{Type{<:LazySet}}(undef, length(partition))
    for (i, block) in enumerate(partition)
        block_options[i] = length(block) == 1 ? Interval : Hyperrectangle
    end
    return block_options
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

function zernikelens()
    topsurface = leaf(AcceleratedParametricSurface(ZernikeSurface(9.5, zcoeff = [(1, 0.0), (4, -1.0), (7, 0.5)]), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0))
    botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)))
    barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64))))
    return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0))
end

function asphericlens()
    topsurface = leaf(AcceleratedParametricSurface(ZernikeSurface(9.5, aspherics = [(2, 0.0), (4, -0.001)]), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0))
    botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)))
    barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64))))
    return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0))
end

function coniclensZ()
    topsurface = leaf(AcceleratedParametricSurface(ZernikeSurface(9.5, radius = -15.0, conic = -5.0), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0))
    botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)))
    barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64))))
    return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0))
end

function coniclensQ()
    topsurface = leaf(AcceleratedParametricSurface(QTypeSurface(9.5, radius = -15.0, conic = -5.0), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0))
    botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)))
    barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64))))
    return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0))
end

function qtypelens()
    topsurface = leaf(AcceleratedParametricSurface(QTypeSurface(9.0, radius = -25.0, conic = 0.3, αcoeffs = [(1, 0, 0.3), (1, 1, 1.0)], βcoeffs = [(1, 0, -0.1), (2, 0, 0.4), (3, 0, -0.6)], normradius = 9.5), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0))
    botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)))
    barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64))))
    return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0))
end

function chebyshevlens()
    topsurface = leaf(AcceleratedParametricSurface(ChebyshevSurface(9.0, 9.0, [(1, 0, -0.081), (2, 0, -0.162), (0, 1, 0.243), (1, 1, 0.486), (2, 1, 0.081), (0, 2, -0.324), (1, 2, -0.405), (2, 2, -0.81)], radius = -25.0), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0))
    botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)))
    barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64))))
    return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0))
end

function gridsaglens()
    topsurface = leaf(AcceleratedParametricSurface(GridSagSurface{Float64}(joinpath(@__DIR__, "../../test/test.GRD"), radius = -25.0, conic = 0.4, interpolation = GridSagBicubic), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0))
    botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)))
    barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64))))
    return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0))
end


text:
function device_commitment(ps_m::CanonicalModel,
                        initial_conditions::Vector{InitialCondition},
                        time_range::UnitRange{Int64},
                        cons_name::Symbol,
                        var_names::Tuple{Symbol,Symbol,Symbol})

    ps_m.constraints[cons_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, ps_m.variables[var_names[1]].axes[1], time_range)

    for i in initial_conditions

        ps_m.constraints[cons_name][i.device.name, 1] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_names[3]][i.device.name, 1] == i.value + ps_m.variables[var_names[1]][i.device.name, 1] - ps_m.variables[var_names[2]][i.device.name, 1])

    end

    for t in time_range[2:end], i in initial_conditions

        ps_m.constraints[cons_name][i.device.name, t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_names[3]][i.device.name, t] == ps_m.variables[var_names[3]][i.device.name, t-1] + ps_m.variables[var_names[1]][i.device.name, t] - ps_m.variables[var_names[2]][i.device.name, t])

    end

    return

end


text:
# #########################################################################
# ##############################   GaN  ###################################
# #########################################################################

# mutable struct GaN <: MixedPotential
#     lj_Ga_Ga::LennardJones
#     lj_N_N::LennardJones
#     bm_Ga_N::BornMayer
#     c::Coulomb
# end

# function GaN()
#     return GaN(LennardJones(), LennardJones(), BornMayer(), Coulomb())
# end

# function get_trainable_params(g::GaN)
#     p = Parameter{(:Ga_ϵ, :Ga_σ)}(values(get_trainable_params(g.lj_Ga_Ga)))
#     p = p ⊕ Parameter{(:N_ϵ, :N_σ)}(values(get_trainable_params(g.lj_Ga_Ga)))
#     p = p ⊕ get_trainable_params(g.bm_Ga_N)
#     return p
# end

# function get_nontrainable_params(g::GaN)
#     p = get_nontrainable_params(g.c)
#     return p
# end

# ########################################################################
# include("energies.jl")
# include("forces.jl")
# include("virials.jl")
# include("gradients.jl")


text:
"""
	project(cmd0::String="", arg1=[], kwargs...)

Project data onto lines or great circles, generate tracks, or translate coordinates.

Full option list at [`project`](http://gmt.soest.hawaii.edu/doc/latest/project.html)

Parameters
----------

- **C** : **origin** : -- Str or list --    Flags = cx/cy

    Sets the origin of the projection, in Definition 1 or 2.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/project.html#c)

- **A** : **azim** : -- Number --    Flags = azimuth

    Defines the azimuth of the projection (Definition 1).
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/project.html#a)
- **E** : **end_point** : -- Str or List --    Flags = bx/by

    bx/by defines the end point of the projection path (Definition 2).
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/project.html#e)
- **F** : **out_flags** : -- wStr --    Flags = xyzpqrs

    Specify your desired output using any combination of xyzpqrs, in any order [Default is xyzpqrs].
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/project.html#f)
- **G** : **no_input** : -- Str or Number --    Flags = dist[/colat][+h]

    Generate mode. No input is read. Create (r, s, p) output points every dist units of p. See Q option.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/project.html#g)
- **L** : **length_control** : -- Str or list --    Flags = [w|l_min/l_max]

    Length controls. Project only those points whose p coordinate is within l_min < p < l_max.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/project.html#l)
- **N** : **flat_earth** : -- Bool or [] --

    Flat Earth. Make a Cartesian coordinate transformation in the plane. [Default uses spherical trigonometry.]
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/project.html#n)
- **Q** : **units** : -- Bool or [] --

    Map type units.
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/project.html#q)
- **S** : **sort** : -- Bool or [] --

    Sort the output into increasing p order. Useful when projecting random data into a sequential profile.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/project.html#s)
- **T** : **pole** : -- Str or list --    Flags = px/py

    px/py sets the position of the rotation pole of the projection. (Definition 3).
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/project.html#t)
- **W** : **width_control** : -- Str or list --    Flags = w_min/w_max

    Width controls. Project only those points whose q coordinate is within w_min < q < w_max.
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/project.html#w)

- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_V)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_s)
- $(GMT.opt_swap_xy)
"""
function project(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("project", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd = parse_V_params("", d)
	cmd, = parse_b(cmd, d)
	cmd, = parse_d(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_g(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_s(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'A', d, [:A :azim])
	cmd = add_opt(cmd, 'C', d, [:C :origin])
	cmd = add_opt(cmd, 'E', d, [:E :end_point])
	cmd = add_opt(cmd, 'F', d, [:F :out_flags])
	cmd = add_opt(cmd, 'G', d, [:G :no_input])
	cmd = add_opt(cmd, 'L', d, [:L :length_control])
	cmd = add_opt(cmd, 'N', d, [:N :flat_earth])
	cmd = add_opt(cmd, 'Q', d, [:Q :units])
	cmd = add_opt(cmd, 'S', d, [:S :sort])
	cmd = add_opt(cmd, 'T', d, [:T :pole])
	cmd = add_opt(cmd, 'W', d, [:W :width_control])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "project", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
project(arg1=[], cmd0::String=""; kw...) = project(cmd0, arg1; kw...)


text:
function dalphainit!(dalpha::Float64, params::Parameters)
  for i in 2:length(params.nodes)
    dz = params.nodes[i].z - params.nodes[i-1].z
    params.nodes[i].inclination = params.nodes[i-1].inclination + dalpha*params.dtor*dz/1200.0
    params.nodes[i].dip = params.nodes[i-1].dip - dalpha*params.dtor*dz/1200.0
    params.elements[i-1].inclination = (params.nodes[i-1].inclination + params.nodes[i].inclination)/2.0
  end
end


text:
# CUDA-specific runtime libraries

import Base.Sys: WORD_SIZE

## GPU runtime library

# reset the runtime cache from global scope, so that any change triggers recompilation
GPUCompiler.reset_runtime()

# load or build the runtime for the most likely compilation job given a compute capability
function precompile_runtime(caps=CUDA.llvm_compat(LLVM.version()).cap)
    dummy_source = FunctionSpec(()->return, Tuple{})
    params = CUDACompilerParams()
    Context() do ctx
        for cap in caps
            # NOTE: this often runs when we don't have a functioning set-up,
            #       so we don't use CUDACompilerTarget(...) which requires NVML
            target = PTXCompilerTarget(; cap=cap)
            job = CompilerJob(target, dummy_source, params)
            GPUCompiler.load_runtime(job; ctx)
        end
    end
    return
end

@eval @inline exception_flag() =
    Base.llvmcall(
        $("""@exception_flag = weak externally_initialized global i$(WORD_SIZE) 0
             define i64 @entry() #0 {
                 %ptr = load i$(WORD_SIZE), i$(WORD_SIZE)* @exception_flag, align 8
                 ret i$(WORD_SIZE) %ptr
             }
             attributes #0 = { alwaysinline }
          """, "entry"), Ptr{Cvoid}, Tuple{})

function signal_exception()
    ptr = exception_flag()
    if ptr !== C_NULL
        unsafe_store!(convert(Ptr{Int}, ptr), 1)
        threadfence_system()
    else
        @cuprintf("""
            WARNING: could not signal exception status to the host, execution will continue.
                     Please file a bug.
            """)
    end
    return
end

function report_exception(ex)
    @cuprintf("""
        ERROR: a %s was thrown during kernel execution.
               Run Julia on debug level 2 for device stack traces.
        """, ex)
    return
end

function report_oom(sz)
    @cuprintf("ERROR: Out of dynamic GPU memory (trying to allocate %i bytes)\n", sz)
    return
end

function report_exception_name(ex)
    @cuprintf("""
        ERROR: a %s was thrown during kernel execution.
        Stacktrace:
        """, ex)
    return
end

function report_exception_frame(idx, func, file, line)
    @cuprintf(" [%i] %s at %s:%i\n", idx, func, file, line)
    return
end

## CUDA device library

function load_libdevice(cap; ctx)
    path = libdevice()
    parse(LLVM.Module, read(path); ctx)
end

function link_libdevice!(mod::LLVM.Module, cap::VersionNumber, undefined_fns)
    ctx = LLVM.context(mod)

    # only link if there's undefined __nv_ functions
    if !any(fn->startswith(fn, "__nv_"), undefined_fns)
        return
    end
    lib::LLVM.Module = load_libdevice(cap; ctx)

    # override libdevice's triple and datalayout to avoid warnings
    triple!(lib, triple(mod))
    datalayout!(lib, datalayout(mod))

    GPUCompiler.link_library!(mod, lib)

    ModulePassManager() do pm
        push!(metadata(mod)["nvvm-reflect-ftz"],
              MDNode([ConstantInt(Int32(1); ctx)]; ctx))
        run!(pm, mod)
    end
end


text:
# GenSPGL
# spgInfo

export spgInfo

"""
GenSPGL output vars composite type. \n

Tip: Use print(ec::spgExitCondition) to display an exit status.
"""
mutable struct spgInfo

    tau::Number
    rNorm::Number
    gNorm::Number
    rErr::Number
    exit_status::spgExitCondition
    iter::Int64
    nProdA::Int64
    nProdAt::Int64
    nNewton::Int64
    timeProject::Float64
    timeMatProd::Float64
    options::spgOptions
    xNorm1::Vector{Number}
    rNorm2::Vector{Number}
    lambda::Vector{Number}

end


text:
function dmrg_x_solver(PH, t, psi0; kwargs...)
  H = contract(PH, ITensor(1.0))
  D, U = eigen(H; ishermitian=true)
  u = uniqueind(U, H)
  max_overlap, max_ind = findmax(abs, array(psi0 * U))
  U_max = U * dag(onehot(u => max_ind))
  return U_max, nothing
end

function dmrg_x(PH, psi0::MPS; reverse_step=false, kwargs...)
  t = Inf
  return tdvp(dmrg_x_solver, PH, t, psi0; reverse_step, kwargs...)
end


text:
# This file was generated, do not modify it. # hide
using Plots
Plots.plot([1,2,3])
Plots.savefig(joinpath(@OUTPUT, "grplot.svg")) # hide


text:
using StructJuMP

m = StructuredModel()

n = 10
@defVar(m, x[1:n] >= 0)
@defVar(m, y[1:n] >= 0)
for i in 1:n
    @addConstraint(m, x[i] - y[i] <= i)
end

p = 3
q = 2
for s = 1:p
    bl = StructuredModel(parent=m)
    @defVar(bl, z[1:n] >= 0)
    @addConstraint(bl, x + sum{z[i], i=1:n} == 1)
    for t = 1:q
        bll = StructuredModel(parent=bl)
        @defVar(bll, w[1:n] >= 0)
        @addConstraint(bll, sum{w[i], i=1:n; iseven(i)} == 1)
        par = parent(bll)
        @addConstraint(bll, sum{z[i], i=1:n} - sum{w[i], i=1:n} >= 0)
        @addConstraint(bll, sum{z[i], i=1:n} >= 0)
        @addConstraint(bll, 0 >= sum{w[i], i=1:n})
    end
end

@addConstraint(m, sum{z, bl in children(m), z=getVar(bl,:z)} == 1)

sum{indx; cond(indx)} (models[indx].y[1] + x[1]) = 0


text:
function aa_comp(aa)
    amino_acids = Dict('W'=>0,'C'=>0,'H'=>0,'M'=>0,'Y'=>0,'Q'=>0,'F'=>0,'N'=>0,'P'=>0,'T'=>0,'R'=>0,'I'=>0,'D'=>0,'G'=>0,'A'=>0,'K'=>0,'E'=>0,'L'=>0,'S'=>0, 'V'=>0)
    total = 0
    for letter in aa
        if haskey(amino_acids, letter)
            amino_acids[letter] += 1
            total += 1
        end
    end
    println(total)
    amino_acids = [(k,v) = (k,v/total) for (k,v) in amino_acids]
    answer = sort([(k,v) = (v,k) for (k,v) in amino_acids])
    return answer
end
aa_comp("MVGGKKKTKICDKVSHEEDRISQLPEPLISEILFHLSTKDSVRTSALSTKWRYLWQSVPG
LDLDPYASSNTNTIVSFVESFFDSHRDSWIRKLRLDLGYHHDKYDLMSWIDAATTRRIQH
LDVHCFHDNKIPLSIYTCTTLVHLRLRWAVLTNPEFVSLPCLKIMHFENVSYPNETTLQK
LISGSPVLEELILFSTMYPKGNVLQLRSDTLKRLDINEFIDVVIYAPLLQCLRAKMYSTK
NFQIISSGFPAKLDIDFVNTGGRYQKKKVIEDILIDISRVRDLVISSNTWKEFFLYSKSR
PLLQFRYISHLNARFYISDLEMLPTLLESCPKLESLILVMSSFNPS")

function aa_pairs()
    aa_string = "WCHMYQFNPTRIDGAKEVLS"
    aa = sort(split(aa_string,""))
    pairs = String[]
    for x = 1:length(aa)-1
        for i = x:length(aa)-1
            push!(pairs, aa[x]*aa[i+1])
        end
    end
    return pairs
end
aa_pairs()

function anti(DNA)
    anti = ""
    if length(DNA) == 0
        return nothing
    end
    complement = Dict([('A','T'),('T','A'),('C','G'),('G','C')])
    for i = 0:length(DNA)-1
        bp = complement[DNA[length(DNA)-i]]
        anti = anti*bp
    end
    return anti
end
anti("ACTGAAAAAAAAAAA")

using Random
function at_seq()
    sequence = ""
    bp = ['A','A','T','T','C','G']
    at = 0
    for i = 1:30
        next_ = rand(bp)
        sequence = sequence*next_
        if next_ == 'A'||next_ == 'T'
            at+=1
        end
    end
    at_fraction = at/30
    length_ = length(sequence)
    return at_fraction, length_, sequence
end
at_seq()

function birthday(people)
    days = Dict([ ("January", 1:31),
    ("February", 1:28),
    ("March", 1:31),
    ("April", 1:30),
    ("May", 1:31),
    ("June", 1:30),
    ("July", 1:31),
    ("August", 1:31),
    ("September", 1:30),
    ("October", 1:31),
    ("November", 1:30),
    ("December", 1:31)
    ])
    repeat = 0
    for x = 1:100
        birthdays = []
        for i = 1:people
        month = rand(days)[1]
        day = rand(days[month])
            push!(birthdays,(month,day))
        end
        if length(Set(birthdays)) != length(birthdays)
            repeat += 1
        end
    end
    return repeat/100
end
birthday(70)

function codons(ORF)
    codons_ = []
    i = 1
    while i < length(ORF) - 1
        push!(codons_, ORF[i:i+2])
        i += 3
        println(i)
    end
    return codons_
end
codons("ATAGCGAATATCTCTCATGAGAGGGAA")

function entropy(array_)
    array_ = [(x*log2(x)) for x in array_]
    sum_ = -1*sum(array_)
    return sum_
end
entropy([.1 .2 .3 .4])

function frame(DNA)
    for i = 1:length(DNA)
        x = string(i)
        y = string(i%3)
        println(x*" "*y*" "*DNA[i])
    end
end
frame("ATGGCCTTT")

function gc(DNA, x)
    count = 0
    for i = 1:length(DNA)
        if DNA[i]=='G'||DNA[i]=='C'
            count += 1
        end
    end
    percent = count/length(DNA)
    return round(percent, digits = x)
end
gc("ACAGAGCCAGCAGATATACAGCAGATACTAT", 2)

function gc_win(DNA, L)
    for i = 1: length(DNA)-(L-1)
        println(string(i)*" "*string(gc(DNA[i:i+L-1],4))*" "*DNA[i:i+L-1])
    end
end
gc_win("ACGACGCAGGAGGAGAGTTTCAGAGATCACGAATACATCCATATTACCCAGAGAGAG", 11)

function gc_win2(seq, L)
    out_ = Matrix{Union{Float64, String, Nothing}}(nothing, length(seq)-(L-1), 3)
    for i = 1:length(seq)-(L-1)
        out_[i, 1] = Float64(i)
        out_[i, 2] = seq[i:i+L-1]
        out_[i, 3] = gc(seq[i:i+L-1], 4)
    end
    return out_
end
gc_win2("ACGACGCAGGAGGAGAGTTTCAGAGATCACGAATACATCCATATTACCCAGAGAGAG", 11)

# no n50
using Statistics
function seqstats(file)
    file = open(file, "r")
    sequences = readlines(file)
    num_seq = 0
    length_ = 0
    seq_len = []
    for line in sequences
        if line[1] != ">"
            num_seq+=1
            length_ += length(line)
            push!(seq_len, length(line))
        end
    end
    min = sort(seq_len)[1]
    max = sort(seq_len)[length(seq_len)]
    avg = mean(seq_len)
    med = median(seq_len)
    return num_seq, length_, min, max, avg, med
end

print("hi")

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
__precompile__()

module GLFW

using Compat
import Compat.String

const lib = Libdl.find_library(["glfw3", "libglfw3", "glfw", "libglfw"], [Pkg.dir("GLFW/deps/usr$(Sys.WORD_SIZE)/lib")])
if isempty(lib)
	error("could not find GLFW library")
end
include_dependency(string(lib, ".", Libdl.dlext)) # Trigger recompilation if the library changes

function GetVersion()
	major, minor, rev = Ref{Cint}(), Ref{Cint}(), Ref{Cint}()
	ccall( (:glfwGetVersion, lib), Void, (Ref{Cint}, Ref{Cint}, Ref{Cint}), major, minor, rev)
	VersionNumber(major[], minor[], rev[])
end
const VERSION = GetVersion()

if VERSION.major == 3
	include("callback.jl")
	include("glfw3.jl")
else
	error("GLFW $VERSION is not supported")
end

function __init__()
	SetErrorCallback((code, description) -> error(description))
	GLFW.Init()
	atexit(GLFW.Terminate)
end

end


text:
function relu(Z::Array{T}) where {T<:AbstractFloat}
    return @. Z * (Z > 0)
end

function logistic(Z::Array{T}) where {T<:AbstractFloat}
    return @. 1 / (1 + exp(-Z))
end

function softmax(Z::Array{T}) where {T<:AbstractFloat}
    e = exp.(Z)
    return e ./ sum(e, dims = 1)
end

function propagate_forward!(NN::NeuralNetwork, mini_batch::MiniBatch)
    gs = NN.params.activation_functions

    NN.cache.Z[1] = NN.params.W[1] * mini_batch.X .+ NN.params.b[1]
    NN.cache.A[1] = gs[1](NN.cache.Z[1])

    for l = 2:NN.L-1
        NN.cache.Z[l] = NN.params.W[l] * NN.cache.A[l-1] .+ NN.params.b[l]
        NN.cache.A[l] = gs[l](NN.cache.Z[l])
    end

    NN.cache.Z[NN.L] = NN.params.W[NN.L] * NN.cache.A[NN.L-1] .+ NN.params.b[NN.L]
    NN.cache.A[NN.L] = gs[NN.L](NN.cache.Z[NN.L])
end

function cross_entropy_binary(Y::Array{T,2}, Y_hat::Array{T,2}) where {T<:AbstractFloat}
    return -sum(Y .* log.(Y_hat) .+ (1 .- Y) .* log.(1 .- Y_hat))
end

function cross_entropy_multi(Y::Array{T,2}, Y_hat::Array{T,2}) where {T<:AbstractFloat}
    return -sum(Y .* log.(Y_hat))
end

function mean_squared_error(Y::Array{T,2}, Y_hat::Array{T,2}) where {T<:AbstractFloat}
    return sum((Y .- Y_hat) .^ 2) / 2
end

function compute_cost!(NN::NeuralNetwork, mini_batch::MiniBatch)
    NN.J[NN.epoch] += NN.params.loss_function(mini_batch.Y, NN.cache.A[NN.L]) / NN.m

    # Add in regularization if using gradient descent
    if NN.hparams.optimization == :gd
        regularization = typeof(NN.m)(0)
        for l = 1:NN.L
            regularization += sum(NN.params.W[l] .^ 2)
        end
        NN.J[NN.epoch] += regularization * NN.hparams.λ / 2 / NN.m
    end
end


text:
###################################################################################################
#                                        Production Conditions
###################################################################################################
function can_wait(actr)
    c1(actr) = actr.visual_location.state.empty
    c2(actr) = actr.visual.state.empty
    c3(actr) = !actr.visual.state.busy
    c4(actr) = !actr.motor.state.busy
    return c1, c2, c3, c4
end

function can_attend(actr)
    c1(actr) = !actr.visual_location.state.empty
    c2(actr) = !actr.visual.state.busy
    return c1, c2
end

function can_respond(actr)
    c1(actr) = !actr.visual.state.empty
    c2(actr) = !actr.motor.state.busy
    return c1, c2
end
###################################################################################################
#                                        Production Actions
###################################################################################################
function wait_action(actr, args...)
    description = "Wait"
    register!(actr, ()->(), now; description)
    return nothing
end

function attend_action(actr, task::PVT, args...)
    buffer = actr.visual_location.buffer
    chunk = deepcopy(buffer[1])
    clear_buffer!(actr.visual_location)
    attending!(actr, chunk, task)
    return nothing
end

function respond_action(actr, task)
    clear_buffer!(actr.visual)
    key = "sb"
    responding!(actr, task, key)
    return nothing
end
###################################################################################################
#                                        Utilities
###################################################################################################
function init_model(scheduler, task, id=1, parms=())
    name = string("model", id)
    procedural = Procedural()
    T = vo_to_chunk() |> typeof
    visual_location = VisualLocation(buffer=T[])
    visual = Visual(buffer=T[])
    visicon = VisualObject[]
    motor = Motor()
    actr = ACTR(;name, scheduler, procedural, visual_location, visual, motor, visicon, parms...)
    rule1 = Rule(;conditions=can_attend, action=attend_action, actr, task, name="Attend")
    push!(procedural.rules, rule1)
    rule2 = Rule(;conditions=can_wait, action=wait_action, actr, task, name="Wait")
    push!(procedural.rules, rule2)
    rule3 = Rule(;conditions=can_respond, action=respond_action, actr, task, name="Respond")
    push!(procedural.rules, rule3)
    return actr
end


text:
#include("../benchmark.jl")

using JuMP, OnePhase
#using Ipopt

function create_tax_problem(na::Int64)
    A = 1:na;

    mu = ones(na);
    lambda = ones(na)

    wmin = 1;
    wmax = wmin + na - 1;
    w = [wmin + ((wmax-wmin)/(na-1))*(i-1) for i = 1:na]

    m = Model()
    @variable(m, c[i in A] >= 0.0)# for (i,j) in T)
    @variable(m, y[i in A] >= 0.0)# for (i,j) in T )

    @NLobjective(m, Min, -sum( lambda[i] * (log(c[i]) - (y[i]/w[i])^(mu[i]+1) / (mu[i]+1)) for i in A) )

    # Incentive
    for i in A
      for p in A
        if i != p
          @NLconstraint(m,(log(c[i]) - (y[i]/w[i])^(mu[i]+1) / (mu[i]+1))
        - (log(c[p]) - (y[p]/w[i])^(mu[i]+1) / (mu[i]+1)) >= 0.0)
        end
      end
    end

    @NLconstraint(m, sum(lambda[i] * (y[i] - c[i]) for i in A)  >= 0.0);

    return m
end

############
m = create_tax_problem(50);
@show m

nlp_raw = MathProgNLPModel(m);
my_pars = Class_parameters();
iter = one_phase_solve(nlp_raw, my_pars);

using Ipopt
setsolver(m,IpoptSolver())
solve(m)

############
m = create_tax_problem(100);
@show m

nlp_raw = MathProgNLPModel(m);
my_pars = Class_parameters();
iter = one_phase_solve(nlp_raw, my_pars);

m = create_tax_problem(250);
@show m

nlp_raw = MathProgNLPModel(m);
my_pars = Class_parameters();
iter = one_phase_solve(nlp_raw, my_pars);

using Ipopt
setsolver(m,IpoptSolver())
solve(m)

############
m = create_tax_problem(500);
@show m

nlp_raw = MathProgNLPModel(m);
my_pars = Class_parameters();
my_par.max_time = 60.0 * 60.0 * 10.0
iter = one_phase_solve(nlp_raw, my_pars);

using Ipopt
setsolver(m,IpoptSolver())
solve(m)

#=
m = create_tax_problem(800);
@show m

nlp_raw = MathProgNLPModel(m);
my_pars = Class_parameters();
my_par.max_time = 60.0 * 60.0 * 10.0
iter = one_phase_solve(nlp_raw, my_pars);
=#
#=
using Ipopt
setsolver(m,IpoptSolver())
solve(m)
=#

#=srand(1)
m = create_tax_problem(5)
one_phase_solve(m)

srand(1)
m = create_tax_problem(10)
one_phase_solve(m)

srand(1)
m = create_tax_problem(20)
one_phase_solve(m)=#
#=
ORG_STDOUT = STDOUT
file = open("test.txt", "w")
redirect_stdout(file)

srand(1)
m = create_tax_problem(30)
one_phase_solve(m)

STDOUT = ORG_STDOUT

close(file)
=#

#srand(1)
#m = create_tax_problem(40)
#one_phase_solve(m)

#=srand(1)
m = create_tax_problem(5)
status = solve(m)

srand(1)
m = create_tax_problem(10)
status = solve(m)

srand(1)
m = create_tax_problem(20)
status = solve(m)=#

#srand(1)
#m = create_tax_problem(40)
#status = solve(m)

#finalize(nlp_raw)


text:
immutable Options
    maxits::Integer
end

immutable PriorPars
    sigma2_beta::Float64
    a_sigma2::Float64
    b_sigma2::Float64
    a_tau::Array{Float64}
    b_tau::Array{Float64}
    a::Float64
    b::Array{Float64}
end

varBayes = function (dat::Dat, priors::PriorPars, options::Options)
	# load the dat
	y = dat.y
	X = dat.X
	K = dat.K
	Z = dat.Z
	z_P = dat.z_P
	# load the priors
	maxIter = options.maxits
	# load the priors
	sigma2_beta = priors.sigma2_beta
	a_sigma2 = priors.a_sigma2
	b_sigma2 = priors.b_sigma2
	a_tau = priors.a_tau
	b_tau = priors.b_tau
	a = priors.a
	b = priors.b

	# Estimation for Beta
	beta_E = rand(1:10, 2)
	# Estimation for sigma2
	sigma2_E = b_sigma2/(a_sigma2 - 1)
	# Estimation for tau
	tau_E = fill(2.1, K) #b_tau/(a_tau - 1)
	tau_ln_E = fill(0, K)
	# Estimation for Kernels and keep them fixed
	inv_Sigma = []
	for k in 1:K
	    jgl_results = inv(cov(Z[:, find(z_P .== k)]))
	    inv_Sigma  = push!(inv_Sigma, jgl_results)

	end
	Ks = zeros(K, n, n)
	for k in 1:K
	    Ks[k, :, :] = calcSigma(Z[:, find(z_P .== k)], inv_Sigma[k])
	end
	# Estimation for H
	H_E = zeros(n, K);
	H_V = zeros(K, n, n);
	for k in 1:K
	    sigma_k = tau_E[k] * Ks[k, :, :]
	    # add a diagonal matrix to the original matrix to make sigma p.d.
	    sigma_k = max(0, -minimum(eigvals(sigma_k))) * I
	    #sigma_k = cholfact(Hermitian(sigma_k))
	    h = rand(MvNormal(tau_E[k] * Ks[k, :, :]), 1)
	    H_E[:, k] = h
	    H_V[k, :, :] = eye(n)
	end
	# Estimation for gamma
	gamma_E = fill(0.5, K)

	"variational Bayes starts now"
	for iter in 1:maxIter
	    # Report progress of mcmc iter
	    if (mod(iter, 500) == 0)
	        println("iter = $iter\n")
	    end

	    # Update Beta
	    beta_V = inv(1/sigma2_beta * eye(p_beta) + 1/sigma2_E * X' * X)
	    beta_E = 1/sigma2_E * beta_V * X' * (y - sum(H_E, 2))
	    #println("beta_E = $beta_E\n")

	    # Update sigma2
	    sum_h = 0
	    for k in 1:K
	        sum_h = sum_h + H_E[:, k]' * H_E[:, k] +  trace(H_V[k, :, :]) - 2 * H_E[:, k]' * y;
	    end
	    idx =  collect(combinations(1:K, 2))
	    for j in 1:length(idx)
	        sum_h = sum_h + 2 * H_E[:, idx[j][1]]' * H_E[:, idx[j][2]]
	    end
	    sum_h = sum_h +  y' * y
	    A = (X * beta_E)' * (X * beta_E) + trace(X * beta_V * X') - 2 * beta_E' * X' * (y - sum(H_E, 2)) + sum_h
	    shape_sigma2 = n/2 + a_sigma2
	    rate_sigma2 = 1/2 * A + b_sigma2
	    sigma2_E = (rate_sigma2/(shape_sigma2 - 1))[1, 1];
	    #println("sigma2_E = $sigma2_E\n")

	    for k in 1:K
	        # Update gamma
	        num  = - 1/2 * 1/sigma2_E * A
	        col_index = setdiff(1:K, k)
	        sum_h = 0
	        for j in col_index
	            sum_h = sum_h + H_E[:, j]' * H_E[:, j] + trace(H_V[j, :, :]) - 2 * H_E[:, j]' * y
	        end
	        idx =  collect(combinations(1:K, 2))
	        for j in 1:length(idx)
	            sum_h = sum_h + 2 * H_E[:, idx[j][1]]' * H_E[:, idx[j][2]]
	        end
	        sum_h = sum_h + y' * y
	        B =  (X * beta_E)' * (X * beta_E) + trace(X * beta_V * X') - 2 * beta_E' * X' * (y - sum(H_E, 2)) + sum_h
	        den = - 1/2 * 1/sigma2_E * B
	        lik_odd = exp(num - den)
	        sigm_1 = exp(a + sum(b[k, col_index] .* gamma_E[col_index]))/ (1 + sum(b[k, col_index] .* gamma_E[col_index]))
	        sigm_0 = 1/ (1 + exp(a + sum(b[k, col_index] .* gamma_E[col_index])))
	        prior_odd = sigm_1/sigm_0
	        odd = (lik_odd * prior_odd)[1,1]
	        gamma_E[k] = odd/(odd + 1)

	        # Update tau
	        shape_tau = n/2 + a_tau[k]
	        B = H_E[:, k]' * inv(Ks[k, :, :]) * H_E[:, k]
	        + trace(inv(Ks[k, :, :]) * H_V[k, :, :])
	        rate_tau = 1/2 * B[1, 1] + b_tau[k]
	        tau_E[k] = rate_tau/(shape_tau - 1)
	        tau_E[k] = gamma_E[k] * tau_E[k] + 0

	        # Update H
	        col_index = setdiff(1:K, k)
	        V_h = inv((1/sigma2_E) * eye(n) + inv(tau_E[k] * Ks[k, :, :]))
	        m_h = 1/sigma2_E * V_h * (y - X * beta_E - sum(H_E[:, col_index], 2))
	        H_E[:, k] = gamma_E[k] * m_h + 0
	        H_V[k, :, :] = gamma_E[k] * gamma_E[k] * V_h + 0
	    end
	    #println("gamma_E = $gamma_E\n")
	    #println("tau_E = $tau_E\n")
	    #println("\n")
	end

	return beta_E, sigma2_E, gamma_E, H_E
end


text:
export Neurite, WhiteSteadyDist, ColourSteadyDist

@with_kw struct Neurite
  τ_v::Float64
  λ::Float64
  M::Int64 # number of spatial steps
  dx::Float64 # spatial step size
end

function x_axis(args_N::Neurite)
  @unpack M,dx = args_N
  x = 0.5*dx:dx:M*dx
end
export x_axis

abstract type DriveParams end

@with_kw struct WhiteSteadyDist <: DriveParams
  μ::Float64
  σ_w::Float64
end

@with_kw struct ColourSteadyDist <: DriveParams
  μ::Float64
  σ_s::Float64
  τ_s::Float64
end

function cfunc(x,ζ,L,λ)
  a = sqrt(1+ζ)
  return (2/a)*cosh((L-x)*a/λ)*cosh(x*a/λ)/sinh(L*a/λ)
end

function sealed_steady_var(N_args::Neurite, I_args::WhiteSteadyDist, x)
  @unpack λ,M,dx = N_args
  @unpack σ_w = I_args
  L = M*dx
  return σ_w^2 .* cfunc.(x,0,L,λ)
end
export sealed_steady_var

function sealed_steady_var(N_args::Neurite, I_args::ColourSteadyDist, x)
  @unpack λ,M,dx,τ_v = N_args
  @unpack σ_s,τ_s = I_args
  α_s = τ_s/τ_v
  L = M*dx
  return σ_s^2*α_s.*(cfunc.(x,0,L,λ).-cfunc.(x,1/α_s,L,λ))
end
export sealed_steady_var

function sealed_steady_dvar(N_args::Neurite, I_args::ColourSteadyDist, x)
  @unpack λ,M,dx,τ_v = N_args
  @unpack σ_s,τ_s = I_args
  α_s = τ_s/τ_v
  L = M*dx
  return σ_s^2 .* cfunc.(x,1/α_s,L,λ)./α_s
end
export sealed_colour_steady_dvar

function sealed_steady_var(N_args::Neurite, I_args::WhiteSteadyDist, T::TimeAxis; seed=0)
  if seed != 0
    Random.seed!(seed)
  end
  @unpack τ_v,λ,M,dx = N_args
  @unpack μ,σ_w = I_args
  @unpack dt, N = T

  V = zeros(M)
  dV = zeros(M+1)
  ζ = zeros(M)
  V2_sum = zeros(M)
  for i=1:N
    randn!(ζ)
    V[1] = (1-dt)*V[1]+λ^2*(dV[2]-dV[1])dt/dx+2*σ_w*sqrt(λ*dt/dx)*ζ[1]
    V2_sum[1] += V[1]^2
    @simd for j=2:M
      @inbounds V[j] = (1-dt)*V[j]+λ^2*(dV[j+1]-dV[j])*dt/dx+2*σ_w*sqrt(λ*dt/dx)*ζ[j]
      @inbounds dV[j] = (V[j]-V[j-1])/dx
      @inbounds V2_sum[j] += V[j]^2
    end
  end
  return V2_sum/N
end
export sealed_steady_var

function sealed_steady_var(N_args::Neurite, I_args::ColourSteadyDist, T::TimeAxis; seed=0)
  if seed != 0
    Random.seed!(seed)
  end
  @unpack τ_v,λ,M,dx = N_args
  @unpack μ,σ_s,τ_s = I_args
  α_s = τ_s/τ_v
  @unpack dt, N = T
  V2 = zeros(M) # new value of V
  V1 = zeros(M) # old value of V
  S = zeros(M) # synaptic drive term
  dV = zeros(M+1)
  ζ = zeros(M)
  V2_sum = zeros(M)
  dV2_sum = zeros(M)
  for i=1:N
    randn!(ζ)
    V2[1] = V1[1]*(1-dt)+λ^2*(dt/dx)*(dV[2]-dV[1])+S[1]*dt
    S[1] = (1-dt/α_s)*S[1]+2*σ_s*ζ[1]*sqrt(λ*dt/(dx*α_s))
    V2_sum[1] += V2[1]^2
    dV2_sum[1] += ((V2[1]-V1[1])/dt)^2
    V1[1] = V2[1]
    @simd for j=2:M
      @inbounds V2[j] = V1[j]*(1-dt)+λ^2*(dt/dx)*(dV[j+1]-dV[j])+S[j]*dt
      @inbounds S[j] = (1-dt/α_s)*S[j]+2*σ_s*ζ[j]*sqrt(λ*dt/(dx*α_s))
      @inbounds dV[j] = (V2[j]-V2[j-1])/dx
      @inbounds V2_sum[j] += V2[j]^2
      @inbounds dV2_sum[j] += ((V2[j]-V1[j])/dt)^2
      @inbounds V1[j] = V2[j]
    end;
  end
  return V2_sum/N, dV2_sum/N
end
export sealed_steady_var

function sealed_white_steady_var_test()
  T = TimeAxis(dt = 1/500, N = 5000000)
  example = Neurite(τ_v = 10.0,λ=200,M=50,dx=20)
  drive = WhiteSteadyDist(μ=5.0,σ_w = 1.5)
  @time Var_sim = sealed_steady_var(example, drive, T; seed=1000)
  x = x_axis(example)
  Var_th = sealed_steady_var(example,drive,x)
  return Var_sim,Var_th,x
end
export sealed_white_steady_var_test

function sealed_colour_steady_var_test()
  T = TimeAxis(dt = 1/500, N = 5000000)
  example = Neurite(τ_v = 10.0,λ=200,M=50,dx=20)
  drive = ColourSteadyDist(μ=5.0,σ_s = 1.5,τ_s=5.0)
  @time Var_sim, dVar_sim = sealed_steady_var(example, drive, T; seed=1000)
  x = x_axis(example)
  Var_th = sealed_steady_var(example,drive,x)
  dVar_th = sealed_steady_dvar(example,drive,x)
  return Var_sim, dVar_sim, Var_th, dVar_th, x
end
export sealed_colour_steady_var_test

max_stars_repo_name: rpgowers/NoisyNeuron

max_issues_repo_name: rpgowers/NoisyNeuron

max_forks_repo_name: rpgowers/NoisyNeuron


text:
"""
Minimum data required to describe dof ownership.
At this moment, the cell-wise ownership is compressed on cell types.
This can be relaxed in the future, to have an arbitrary cell-wise dof ownership.
"""
struct CellConformity{T} <: GridapType
  cell_ctype::T
  ctype_lface_own_ldofs::Vector{Vector{Vector{Int}}}
  ctype_lface_pindex_pdofs::Vector{Vector{Vector{Vector{Int}}}}
  d_ctype_num_dfaces::Vector{Vector{Int}}
end

function Base.getproperty(a::CellConformity, sym::Symbol)
  if sym == :d_ctype_offset
    _d_ctype_offset(a)
  elseif sym == :d_ctype_ldface_own_ldofs
    _d_ctype_ldface_own_ldofs(a)
  else
    getfield(a, sym)
  end
end

function Base.propertynames(x::CellConformity, private::Bool=false)
  (fieldnames(typeof(x))...,:d_ctype_offset,:d_ctype_ldface_own_ldofs)
end

function _d_ctype_offset(a::CellConformity)
  num_ctypes = length(a.ctype_lface_own_ldofs)
  num_ds = length(a.d_ctype_num_dfaces)
  d_ctype_offset = [ zeros(Int,num_ctypes) for d in 1:num_ds]
  for d in 2:num_ds
    for ctype in 1:num_ctypes
      d_ctype_offset[d][ctype] = d_ctype_offset[d-1][ctype] + a.d_ctype_num_dfaces[d-1][ctype]
    end
  end
  d_ctype_offset
end

function _d_ctype_ldface_own_ldofs(a::CellConformity)
  num_ctypes = length(a.ctype_lface_own_ldofs)
  num_ds = length(a.d_ctype_num_dfaces)
  [[[ a.ctype_lface_own_ldofs[ctype][ldface+a.d_ctype_offset[d][ctype]]
    for ldface in 1:a.d_ctype_num_dfaces[d][ctype] ]
    for ctype in 1:num_ctypes ]
    for d in 1:num_ds ]
end

"""
Minimum data required to build a conforming FE space.
At this moment, the some cell-wise info is compressed on cell types.
This can be relaxed in the future, and have an arbitrary cell-wise data.
"""
struct CellFE{T} <: GridapType
  cell_ctype::T
  ctype_num_dofs::Vector{Int}
  ctype_ldof_comp::Vector{Vector{Int}}
  cell_conformity::CellConformity{T}
  cell_shapefuns::AbstractArray{<:AbstractVector{<:Field}}
  cell_dof_basis::AbstractArray{<:AbstractVector{<:Dof}}
  cell_shapefuns_domain::DomainStyle
  cell_dof_basis_domain::DomainStyle
  max_order::Int
end
# If the shapefuns are not polynomials, max_order has to be understood as the order of a
# reasonable quadrature rule to integrate the shape functions. Only used by FESpace
# constructors that need to integrate the shape functions (e.g., ZeroMeanFESpace).

Geometry.num_cells(cell_fe::CellFE) = length(cell_fe.cell_ctype)

function CellConformity(cell_fe::CellFE)
  cell_fe.cell_conformity
end

"""
Generate A CellConformity from a vector of reference fes
"""
function CellConformity(cell_reffe::AbstractArray{<:ReferenceFE},conf::Conformity)
  ctype_reffe, cell_ctype = compress_cell_data(cell_reffe)
  ctype_lface_own_ldofs = map(reffe->get_face_own_dofs(reffe,conf),ctype_reffe)
  ctype_lface_pindex_pdofs = map(reffe->get_face_own_dofs_permutations(reffe,conf),ctype_reffe)
  D = num_dims(first(ctype_reffe))
  d_ctype_num_dfaces = [ map(reffe->num_faces(get_polytope(reffe),d),ctype_reffe) for d in 0:D]
  CellConformity(
    cell_ctype,
    ctype_lface_own_ldofs,
    ctype_lface_pindex_pdofs,
    d_ctype_num_dfaces)
end

"""
Generate a CellFE from a vector of reference fes
"""
function CellFE(
  model::DiscreteModel,
  cell_reffe::AbstractArray{<:ReferenceFE},
  conformity::Conformity
 )
  cell_conformity = CellConformity(cell_reffe,conformity)
  ctype_reffe, cell_ctype = compress_cell_data(cell_reffe)
  ctype_num_dofs = map(num_dofs,ctype_reffe)
  ctype_ldof_comp = map(reffe->get_dof_to_comp(reffe),ctype_reffe)
  cell_shapefuns = get_cell_shapefuns(model,cell_reffe,conformity)
  cell_dof_basis = get_cell_dof_basis(model,cell_reffe,conformity)
  cell_shapefuns_domain = ReferenceDomain()
  cell_dof_basis_domain = cell_shapefuns_domain
  max_order = maximum(map(get_order,ctype_reffe))

  CellFE(
    cell_ctype,
    ctype_num_dofs,
    ctype_ldof_comp,
    cell_conformity,
    cell_shapefuns,
    cell_dof_basis,
    cell_shapefuns_domain,
    cell_dof_basis_domain,
    max_order)
end

function get_cell_dof_basis(model::DiscreteModel,
                            cell_reffe::AbstractArray{<:ReferenceFE},
                            ::Conformity)
  lazy_map(get_dof_basis,cell_reffe)
end

function get_cell_shapefuns(model::DiscreteModel,
                            cell_reffe::AbstractArray{<:ReferenceFE},
                            ::Conformity)
  lazy_map(get_shapefuns,cell_reffe)
end

# Low level conforming FE Space constructor
# The user is not expected to call this function. Use the factory function instead
function _ConformingFESpace(
  vector_type::Type,
  model::DiscreteModel,
  face_labeling::FaceLabeling,
  cell_fe::CellFE,
  dirichlet_tags,
  dirichlet_components)

  grid_topology = get_grid_topology(model)
  ntags = length(dirichlet_tags)

  cell_dofs_ids, nfree, ndirichlet, dirichlet_dof_tag, dirichlet_cells = compute_conforming_cell_dofs(
    cell_fe,CellConformity(cell_fe),grid_topology,face_labeling,dirichlet_tags,dirichlet_components)

  trian = Triangulation(model)
  cell_shapefuns, cell_dof_basis = compute_cell_space(cell_fe,trian)

  cell_is_dirichlet = fill(false,num_cells(trian))
  cell_is_dirichlet[dirichlet_cells] .= true

  UnconstrainedFESpace(
    vector_type,
    nfree,
    ndirichlet,
    cell_dofs_ids,
    cell_shapefuns,
    cell_dof_basis,
    cell_is_dirichlet,
    dirichlet_dof_tag,
    dirichlet_cells,
    ntags)
end

function compute_cell_space(cell_shapefuns::AbstractArray{<:AbstractArray{<:Field}},
                            cell_dof_basis::AbstractArray{<:AbstractArray{<:Dof}},
                            cell_shapefuns_domain::DomainStyle,
                            cell_dof_basis_domain::DomainStyle,
                            trian::Triangulation)
    SingleFieldFEBasis(cell_shapefuns,trian,TestBasis(),cell_shapefuns_domain),
    CellDof(cell_dof_basis,trian,cell_dof_basis_domain)
end

function compute_cell_space(cellfe::CellFE,trian::Triangulation)
    compute_cell_space(cellfe.cell_shapefuns,
                       cellfe.cell_dof_basis,
                       cellfe.cell_shapefuns_domain,
                       cellfe.cell_dof_basis_domain,
                       trian)
end

function compute_cell_space(cell_fe,trian::Triangulation)
  cell_shapefuns, cell_dof_basis, d1, d2 = _compute_cell_space(cell_fe)
  SingleFieldFEBasis(cell_shapefuns,trian,TestBasis(),d1), CellDof(cell_dof_basis,trian,d2)
end

function _compute_cell_space(cell_fe::CellFE)
  ( cell_fe.cell_shapefuns,
    cell_fe.cell_dof_basis,
    cell_fe.cell_shapefuns_domain,
    cell_fe.cell_dof_basis_domain)
end

function _compute_cell_space(
  cell_reffe::AbstractArray{<:ReferenceFE})
  cell_shapefuns = lazy_map(get_shapefuns,cell_reffe)
  cell_dof_basis = lazy_map(get_dof_basis,cell_reffe)
  cell_shapefuns, cell_dof_basis, ReferenceDomain(), ReferenceDomain()
end

"""
The result is the tuple

    (cell_dofs, nfree, ndiri, dirichlet_dof_tag, dirichlet_cells)

If `dirichlet_components`  is given, then `get_dof_to_comp` has to be defined
for the reference elements in `reffes`.
"""
function compute_conforming_cell_dofs(
  cell_fe,cell_conformity,grid_topology,face_labeling,dirichlet_tags,dirichlet_components=nothing)

  cell_to_ctype = cell_fe.cell_ctype
  ctype_to_ldof_to_comp = cell_fe.ctype_ldof_comp
  ctype_to_num_dofs = cell_fe.ctype_num_dofs

  d_to_ctype_to_ldface_to_own_ldofs = cell_conformity.d_ctype_ldface_own_ldofs
  ctype_to_lface_to_own_ldofs = cell_conformity.ctype_lface_own_ldofs
  ctype_to_lface_to_pindex_to_pdofs = cell_conformity.ctype_lface_pindex_pdofs

  D = num_cell_dims(grid_topology)
  n_faces = num_faces(grid_topology)
  d_to_cell_to_dfaces = [ Table(get_faces(grid_topology,D,d)) for d in 0:D]
  d_to_dface_to_cells = [ Table(get_faces(grid_topology,d,D)) for d in 0:D]
  d_to_offset = get_offsets(grid_topology)

  face_to_own_dofs, ntotal, d_to_dface_to_cell, d_to_dface_to_ldface =  _generate_face_to_own_dofs(
     n_faces,
     cell_to_ctype,
     d_to_cell_to_dfaces,
     d_to_dface_to_cells,
     d_to_offset,
     d_to_ctype_to_ldface_to_own_ldofs)

  d_to_dface_to_tag = [ get_face_tag_index(face_labeling,dirichlet_tags,d)  for d in 0:D]
  cell_to_faces = Table(get_cell_faces(grid_topology))

  _dirichlet_components = _convert_dirichlet_components(dirichlet_tags,dirichlet_components)

  nfree, ndiri, diri_dof_tag = _split_face_own_dofs_into_free_and_dirichlet!(
    face_to_own_dofs,
    d_to_offset,
    d_to_dface_to_tag,
    d_to_dface_to_cell,
    d_to_dface_to_ldface,
    cell_to_ctype,
    d_to_ctype_to_ldface_to_own_ldofs,
    ctype_to_ldof_to_comp,
    _dirichlet_components)

  cell_to_lface_to_pindex = Table(get_cell_permutations(grid_topology))

  cell_dofs = CellDofsNonOriented(
    cell_to_faces,
    cell_to_lface_to_pindex,
    cell_to_ctype,
    ctype_to_lface_to_own_ldofs,
    ctype_to_num_dofs,
    face_to_own_dofs,
    ctype_to_lface_to_pindex_to_pdofs)

  diri_cells = _generate_diri_cells(
    d_to_dface_to_tag,
    d_to_cell_to_dfaces)

  (Table(cell_dofs), nfree, ndiri, diri_dof_tag, diri_cells)
end

function _convert_dirichlet_components(dirichlet_tags::AbstractArray,dirichlet_components)
  dirichlet_components
end

function _convert_dirichlet_components(dirichlet_tags::Integer,dirichlet_components)
  dirichlet_components==nothing ? nothing : [dirichlet_components,]
end

function _convert_dirichlet_components(dirichlet_tags::AbstractString,dirichlet_components)
  dirichlet_components==nothing ? nothing : [dirichlet_components,]
end

function _generate_face_to_own_dofs(
  n_faces,
  cell_to_ctype,
  d_to_cell_to_dfaces::Vector{Table{T,Vd,Vp}},
  d_to_dface_to_cells::Vector{Table{T,Vd,Vp}},
  d_to_offset,
  d_to_ctype_to_ldface_to_own_ldofs) where {T,Vd,Vp}

  P=eltype(Vp)
  face_to_own_dofs_ptrs = zeros(P,n_faces+1)

  D = length(d_to_offset)-1

  icell = 1
  d_to_dface_to_cell = [ get_local_item(d_to_dface_to_cells[d+1],icell)  for d in 0:D ]

  d_to_dface_to_ldface = [
    find_local_index(d_to_dface_to_cell[d+1],d_to_cell_to_dfaces[d+1]) for d in 0:D ]

  for d in 0:D
    cell_to_dfaces = d_to_cell_to_dfaces[d+1]
    dface_to_cells = d_to_dface_to_cells[d+1]
    offset = d_to_offset[d+1]
    ctype_to_ldface_to_own_ldofs = d_to_ctype_to_ldface_to_own_ldofs[d+1]
    ctype_to_ldface_to_num_own_ldofs = map( (x) -> length.(x) ,ctype_to_ldface_to_own_ldofs)
    dface_to_cell_owner = d_to_dface_to_cell[d+1]
    dface_to_ldface = d_to_dface_to_ldface[d+1]

    if any( ctype_to_ldface_to_num_own_ldofs .!= 0)
      _generate_face_to_own_dofs_count_d!(
        face_to_own_dofs_ptrs,
        offset,
        cell_to_ctype,
        dface_to_cell_owner,
        dface_to_ldface,
        ctype_to_ldface_to_num_own_ldofs)
    end
  end

  length_to_ptrs!(face_to_own_dofs_ptrs)

  n_dofs = face_to_own_dofs_ptrs[end]-1
  face_to_own_dofs_data = collect(T(1):T(n_dofs))

  face_to_own_dofs = Table(face_to_own_dofs_data,face_to_own_dofs_ptrs)
  (face_to_own_dofs, n_dofs, d_to_dface_to_cell, d_to_dface_to_ldface)
end

function  _generate_face_to_own_dofs_count_d!(
  face_to_own_dofs_ptrs,
  offset,
  cell_to_ctype,
  dface_to_cell_owner,
  dface_to_ldface,
  ctype_to_ldface_to_num_own_ldofs)

  n_dfaces = length(dface_to_ldface)
  for dface in 1:n_dfaces
    cell = dface_to_cell_owner[dface]
    ldface = dface_to_ldface[dface]
    ctype = cell_to_ctype[cell]
    n_own_ldofs = ctype_to_ldface_to_num_own_ldofs[ctype][ldface]
    face = dface + offset
    face_to_own_dofs_ptrs[face+1] = n_own_ldofs
  end
end

function _split_face_own_dofs_into_free_and_dirichlet!(
  face_to_own_dofs,
  d_to_offset,
  d_to_dface_to_tag,
  d_to_dface_to_cell,
  d_to_dface_to_ldface,
  cell_to_ctype,
  d_to_ctype_to_ldface_to_own_ldofs,
  ctype_to_ldof_to_comp,
  dirichlet_components::Nothing)

  _split_face_own_dofs_into_free_and_dirichlet_generic!(
    face_to_own_dofs,
    d_to_offset,
    d_to_dface_to_tag)

end

function _split_face_own_dofs_into_free_and_dirichlet!(
  face_to_own_dofs,
  d_to_offset,
  d_to_dface_to_tag,
  d_to_dface_to_cell,
  d_to_dface_to_ldface,
  cell_to_ctype,
  d_to_ctype_to_ldface_to_own_ldofs,
  ctype_to_ldof_to_comp,
  dirichlet_components)

  D = length(d_to_dface_to_cell)-1

  _split_face_own_dofs_into_free_and_dirichlet_with_components!(
    face_to_own_dofs,
    d_to_offset,
    d_to_dface_to_tag,
    d_to_dface_to_cell,
    d_to_dface_to_ldface,
    cell_to_ctype,
    ctype_to_ldof_to_comp,
    d_to_ctype_to_ldface_to_own_ldofs,
    dirichlet_components)
end

function _split_face_own_dofs_into_free_and_dirichlet_generic!(
  face_to_own_dofs,
  d_to_offset,
  d_to_dface_to_tag)

  D = length(d_to_offset)-1
  nfree = 0
  ndiri = 0
  dirichlet_dof_tag = Int8[]
  for d in 0:D
    dface_to_tag = d_to_dface_to_tag[d+1]
    offset = d_to_offset[d+1]
    for (dface, tag) in enumerate(dface_to_tag)
      face = dface + offset
      pini = face_to_own_dofs.ptrs[face]
      pend = face_to_own_dofs.ptrs[face+1]-1
      if tag == UNSET
        for p in pini:pend
          nfree += 1
          face_to_own_dofs.data[p] = nfree
        end
      else
        for p in pini:pend
          ndiri += 1
          face_to_own_dofs.data[p] = -ndiri
          push!(dirichlet_dof_tag,tag)
        end
      end
    end
  end

  (nfree, ndiri, dirichlet_dof_tag)
end

function _split_face_own_dofs_into_free_and_dirichlet_with_components!(
  face_to_own_dofs,
  d_to_offset,
  d_to_dface_to_tag,
  d_to_dface_to_cell,
  d_to_dface_to_ldface,
  cell_to_ctype,
  ctype_to_ldof_to_comp,
  d_to_ctype_to_ldface_to_own_ldofs,
  dirichlet_components)

  D = length(d_to_offset)-1
  nfree = 0
  ndiri = 0
  dirichlet_dof_tag = Int8[]
  for d in 0:D
    dface_to_tag = d_to_dface_to_tag[d+1]
    offset = d_to_offset[d+1]
    dface_to_cell = d_to_dface_to_cell[d+1]
    dface_to_ldface = d_to_dface_to_ldface[d+1]
    ctype_to_ldface_to_own_ldofs = d_to_ctype_to_ldface_to_own_ldofs[d+1]
    ctype_to_ldface_to_num_own_ldofs = map( (x) -> length.(x) ,ctype_to_ldface_to_own_ldofs)
    if all(ctype_to_ldface_to_num_own_ldofs .== 0)
      continue
    end
    for (dface, diritag) in enumerate(dface_to_tag)
      face = dface + offset
      pini = face_to_own_dofs.ptrs[face]
      pend = face_to_own_dofs.ptrs[face+1]-1
      if diritag == UNSET
        for p in pini:pend
          nfree += 1
          face_to_own_dofs.data[p] = nfree
        end
      else
        cell = dface_to_cell[dface]
        ldface = dface_to_ldface[dface]
        ctype = cell_to_ctype[cell]
        ldof_to_comp = ctype_to_ldof_to_comp[ctype]
        own_ldofs = ctype_to_ldface_to_own_ldofs[ctype][ldface]
        comp_to_isdiri = dirichlet_components[diritag]
        for (fldof, p) in enumerate(pini:pend)
          ldof = own_ldofs[fldof]
          comp = ldof_to_comp[ldof]
          isdiri = comp_to_isdiri[comp]
          if isdiri
            ndiri += 1
            face_to_own_dofs.data[p] = -ndiri
            push!(dirichlet_dof_tag,diritag)
          else
            nfree += 1
            face_to_own_dofs.data[p] = nfree
          end
        end
      end
    end
  end

  (nfree, ndiri, dirichlet_dof_tag)
end

function _generate_diri_cells(
  d_to_dface_to_tag,
  d_to_cell_to_dfaces)

  ncells = length(first(d_to_cell_to_dfaces))
  cell_visited = fill(false,ncells)

  diri_cells = Int32[]

  D = length(d_to_dface_to_tag)-1
  for d in 0:D
    dface_to_tag = d_to_dface_to_tag[d+1]
    cell_to_dfaces = Table(d_to_cell_to_dfaces[d+1])
    for cell in 1:length(cell_to_dfaces)
      if cell_visited[cell]
        continue
      end
      pini = cell_to_dfaces.ptrs[cell]
      pend = cell_to_dfaces.ptrs[cell+1]-1
      for p in pini:pend
        dface = cell_to_dfaces.data[p]
        tag = dface_to_tag[dface]
        if tag != UNSET
          push!(diri_cells,cell)
          cell_visited[cell] = true
          break
        end
      end
    end
  end

  diri_cells

end

struct CellDofsNonOriented <:AbstractVector{Vector{Int32}}
  cell_to_faces::Table{Int32,Vector{Int32},Vector{Int32}}
  cell_to_lface_to_pindex::Table{Int8,Vector{Int8},Vector{Int32}}
  cell_to_ctype::Vector{Int8}
  ctype_to_lface_to_own_ldofs::Vector{Vector{Vector{Int}}}
  ctype_to_num_dofs::Vector{Int}
  face_to_own_dofs::Table{Int32,Vector{Int32},Vector{Int32}}
  ctype_to_lface_to_pindex_to_pdofs::Vector{Vector{Vector{Vector{Int}}}}
end

Base.size(a::CellDofsNonOriented) = (length(a.cell_to_faces),)

Base.IndexStyle(::Type{<:CellDofsNonOriented}) = IndexLinear()

function array_cache(a::CellDofsNonOriented)
  n_dofs = testitem(a.ctype_to_num_dofs)
  T = eltype(eltype(a))
  v = zeros(T,n_dofs)
  CachedArray(v)
end

function getindex!(cache,a::CellDofsNonOriented,cell::Integer)
  ctype = a.cell_to_ctype[cell]
  n_dofs = a.ctype_to_num_dofs[ctype]
  setsize!(cache,(n_dofs,))
  dofs = cache.array
  lface_to_own_ldofs = a.ctype_to_lface_to_own_ldofs[ctype]
  p = a.cell_to_faces.ptrs[cell]-1
  for (lface, own_ldofs) in enumerate(lface_to_own_ldofs)
    face = a.cell_to_faces.data[p+lface]
    pindex = a.cell_to_lface_to_pindex.data[p+lface]
    pdofs = a.ctype_to_lface_to_pindex_to_pdofs[ctype][lface][pindex]

    q = a.face_to_own_dofs.ptrs[face]-1
    for (i,ldof) in enumerate(own_ldofs)
      j = pdofs[i]
      dof = a.face_to_own_dofs.data[q+j]
      dofs[ldof] = dof
    end
  end
  dofs
end

function Base.getindex(a::CellDofsNonOriented,cell::Integer)
  cache = array_cache(a)
  getindex!(cache,a,cell)
end


text:
using Base.Test
using Rotations
using StaticArrays

# TODO test mean()

include("rotation_tests.jl")
include("derivative_tests.jl")


text:
const p2 = 2
function pow_array2(x::Vector{Float64})
s = 0.0
for y in x
s = s + y^p2
end
return s
end

max_stars_repo_path: Module 2/Chapter04/pow_Array2.jl

max_issues_repo_path: Module 2/Chapter04/pow_Array2.jl

max_forks_repo_path: Module 2/Chapter04/pow_Array2.jl


text:
function undef_copy(v::Array{T, d}) where {T, d}
    Array{T, d}(undef, size(v))
end

function make_volterra_kernel(K::Array{T, 2}) where {T <: Number}
    n = size(K, 1)

    V = LowerTriangular(Array{T, 2}(undef, (n, n)))

    V[1,1] = 0

    for i  = 2:n
        V[i,1] = K[i,1]
        V[i,i] = K[i,i]

        for j = 2:(i-1)
            V[i,j] = 2 * K[i,j]
        end
    end

    return(V)
end


text:
module Reader
    using GZip
    using BinDeps

    export
        readimages,
        readlabels

    # Constants

    const IMAGEOFFSET = 16
    const LABELOFFSET = 8

    # Includes

    include("readheader.jl")
    include("readimages.jl")
    include("readlabels.jl")
end


text:
module DeconvOptim

export deconvolution

 # for optimization
using Optim
 #mean
using Statistics
using FFTW
 # for fast array regularizers
using Tullio
 # possible up_sampling
using Interpolations
 # to check wether CUDA is enabled
using Requires
# for defining custom derivatives
using ChainRulesCore

# via require we can check whether CUDA is loaded
# to enable CUDA support simply load CUDA before load DeconvOptim
 #to_gpu_or_not(x) = x
 #function __init__()
 #    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
 #        print("CUDA support is enabled")
 #        @eval using CUDA
 #        @eval to_gpu_or_not(x) = CuArray(x)
 #    end
 #end

include("forward_models.jl")
include("lossfunctions.jl")
include("mappings.jl")
include("regularizer.jl")
include("utils.jl")

# refresh Zygote to load the custom rrules defined with ChainRulesCore
using Zygote: @adjoint, gradient

"""
    deconvolution(measured, psf; <keyword arguments>)
Computes the deconvolution of `measured` and `psf`. Return parameter is a tuple with
two elements. The first entry is the deconvolved image. The second return parameter
is the output of the optimization of Optim.jl

Multiple keyword arguments can be specified for different loss functions,
regularizers and mappings.

# Arguments
- `loss=Poisson()`: the loss function taking a vector the same shape as measured.
- `regularizer=GR()`: A regularizer function, same form as `loss`.
- `λ=0.05`: A float indicating the total weighting of the regularizer with
    respect to the global loss function
- `background=0`: A float indicating a background intensity level.
- `mapping=Non_negative()`: Applies a mapping of the optimizer weight. Default is a
              parabola which achieves a non-negativity constraint.
- `iterations=20`: Specifies a number of iterations after the optimization.
    definitely should stop.
- `plan_fft=true`: Boolean whether plan_fft is used. Gives a slight speed improvement.
- `padding=0`: an float indicating the amount (fraction of the size in that dimension)
        of padded regions around the reconstruction. Prevents wrap around effects of the FFT.
        A array with `size(arr)=(400, 400)` with `padding=0.05` would result in reconstruction size of
        `(440, 440)`. However, we only return the reconstruction cropped to the original size.
        `padding=0` disables any padding.
- `optim_options=nothing`: Can be a options file required by Optim.jl. Will overwrite iterations.
- `optim_optimizer=LBFGS()`: The chosen Optim.jl optimizer.

# Example
```julia-repl
julia> using DeconvOptim, TestImages, Images, FFTW, Noise, ImageView;

julia> img = channelview(testimage("resolution_test_512"));

julia> psf = generate_psf(size(img), 30);

julia> img_b = conv_psf(img, psf);

julia> img_n = poisson(img_b, 300);

julia> @time res, o = deconvolution(img_n, psf);

julia> colorview(Gray, [img img_n res]) |> imshow
[...]
```
"""
function deconvolution(measured::AbstractArray{T, N}, psf;
        loss=Poisson(),
        regularizer=GR(),
        λ=0.05,
        background=0,
        mapping=Non_negative(),
        iterations=20,
        plan_fft=true,
        padding=0.00,
        optim_options=nothing,
        optim_optimizer=LBFGS(),
        ) where {T, N}

    # do some type conversion to ensure same type everywhere
    # provides speed-up
    λ = convert(eltype(measured), λ)
    psf = convert(Array{eltype(measured)}, psf)
    background = convert(eltype(measured), background)

    # Get the mapping functions to achieve constraints
    # like non negativity
    if mapping != nothing
        mf, m_invf = mapping
    else
        m_invf = identity
    end

    # rec0 will be an array storing the final reconstruction
    # we choose it larger than the measured array to reduce
    # wrap around artifacts of the Fourier Transform
    # we create a array size_padded which stores a new array size
    # our reconstruction array will be larger than measured
    # to prevent wrap around artifacts
    size_padded = []
    for i = 1:ndims(measured)
        # if the size of the i-th dimension is 1
        # don't do any padding because there won't be no
        # convolution in that dimension
        if  size(measured)[i] == 1
            push!(size_padded, 1)
        else
            # only pad, if padding is true
            if ~iszero(padding)
                # 2 * ensures symmetric padding
                # minimum padding is 2 (4 in total) on each side
                x = max(4, 2 * round(Int, size(measured)[i] * padding))
            else
                x = 0
            end
            push!(size_padded, size(measured)[i] + x)
        end
    end

    # the dimensions we do the Fourier Transform over
    fft_dims = collect(1:ndims(psf))

    # we divide by the maximum to normalize
    rescaling = maximum(measured)
    measured = measured ./ rescaling
    # create rec0 which will be the initial guess for the reconstruction
    rec0 = zeros(T, (size_padded...))

    # alternative rec0_center, unused at the moment
    #rec0_center = m_invf(abs.(conv_psf(measured, psf, fft_dims)))
    #
    # take the mean as the initial guess
    # therefore has the same total energy at the initial guess as
    # measured
    rec0_center = mean(measured) .* ones(T, size(measured))
    center_set!(rec0, rec0_center)

    # psf_n is the psf with the same size as rec0 but only in that dimensions
    # that were supported by the initial psf. Broadcasting of psf with less
    # dimensions is still supported
    # we put the small psf into the new one
    # it is important to pad the PSF instead of the OTF

    psf_new_size = Array{Int}(undef, 0)
    for i = 1:ndims(psf)
        push!(psf_new_size, size(rec0)[i])
    end

    psf_new_size = tuple(psf_new_size...)
    psf_n = zeros(eltype(rec0), psf_new_size)
    psf_n = center_set!(psf_n, fftshift(psf))
    psf = ifftshift(psf_n)

    # the psf should be normalized to 1
    psf ./= sum(psf)

    # use plan_fft (function of the FFTW.jl has the same name)
    # for speed improvement
    if plan_fft
        # otf is obtained by rfft(psf)
        # therefore size(psf) != size(otf)
        otf, conv = plan_conv_r(psf, rec0, fft_dims)
    else
        otf = rfft(psf, fft_dims)
        conv(rec, otf) = conv_otf_r(rec, otf, fft_dims)
    end

    # forward model is a convolution
    # due to numerics, we need to clip at 0
    # analytically it's a convolution psf ≥ 0 and image ≥ 0
    # so it must be conv(psf, image) ≥ 0
    forward(x) = (conv_aux(conv, x, otf)) .+ background
    # create the loss function which depends simply on the current rec
    function total_loss(rec)
        # handle if there is a provided mapping function
        if mapping != nothing
            mf_rec = mf(rec)
        else
            mf_rec = rec
        end
        forward_v = forward(mf_rec)
        loss_v = loss(forward_v, measured)
        # handle if there is a regularizer
        if regularizer != nothing
            reg_v = regularizer(mf_rec)
            out = loss_v + λ * reg_v
        else
            out = loss_v
        end
        return out
    end

    # this is the function which will be provided to Optimize
    # check Optim's documentation for the purpose of F and Get
    # but simply speaking F is the loss value and G it's gradient
    # depending whether one of them is nothing, we skip some computations
    # we need to call Base.invokelatest because the regularizer is a function
    # generated at runtime with eval.
    # This leads to the common "world age problem" in Julia
    # for more details on that check:
    # https://discourse.julialang.org/t/dynamically-create-a-function-initial-idea-with-eval-failed-due-to-world-age-issue/49139/17
    function f!(F, G, rec)
        # Zygote calculates both derivative and loss, therefore do everything in one step
        if G != nothing
            y, back = Base.invokelatest(Zygote._pullback, total_loss, rec)
            # calculate gradient
            G .= Base.invokelatest(back, 1)[2]
            if F != nothing
                return y
            end
        end
        if F != nothing
            return Base.invokelatest(total_loss, rec)
        end
    end

    # if not special options are given, just restrict iterations
    if optim_options == nothing
        optim_options = Optim.Options(iterations=iterations)
    end

    # do the optimization with LBGFS
    res = Optim.optimize(Optim.only_fg!(f!), rec0, optim_optimizer, optim_options)

    # apply the mapping to get the real result
    if mapping != nothing
        res_out = mf(Optim.minimizer(res))
    else
        res_out = Optim.minimizer(res)
    end

    res_out .*= rescaling
    # since we do some padding we need to extract the center part
    res_out = center_extract(res_out, size(measured))
    return res_out, res
end

# end module
end