text:
"""
$(SIGNATURES)

Equilibrium part of gas kinetic flux which works with fluid-particle simulations
"""
function flux_equilibrium!(
    fw::T1,
    wL::T2,
    wR::T2,
    inK::Real,
    γ::Real,
    μᵣ::Real,
    ω::Real,
    Pr::Real,
    dt::Real,
    dxL::Real,
    dxR::Real,
    swL = zeros(eltype(fw), axes(wL))::T1,
    swR = zeros(eltype(fw), axes(wR))::T1,
) where {T1<:AA{<:FN,1},T2<:AA{<:Real,1}} # 1D

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    Mu1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

    w =
        primL[1] .* moments_conserve(MuL1, Mxi1, 0, 0) .+
        primR[1] .* moments_conserve(MuR2, Mxi2, 0, 0)
    prim = conserve_prim(w, γ)

    # positivity check
    if prim[end] < 0.0
        primL = conserve_prim(wL .- swL .* dxL, γ)
        primR = conserve_prim(wR .+ swR .* dxR, γ)

        Mu1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
        Mu2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

        w =
            primL[1] .* moments_conserve(MuL1, Mxi1, 0, 0) .+
            primR[1] .* moments_conserve(MuR2, Mxi2, 0, 0)
        prim = conserve_prim(w, γ)

        tau = vhs_collision_time(prim, μᵣ, ω) #+
        #2.0 * dt * abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        #(primL[1] / primL[end] + primR[1] / primR[end])

        # time-integration constants
        Mt = zeros(2)
        Mt[2] = tau * (1.0 - exp(-dt / tau))
        Mt[1] = dt - Mt[2]

        # flux related to central distribution
        Mu, Mxi, MuL, MuR = gauss_moments(prim, inK)
        Muv = moments_conserve(Mu, Mxi, 1, 0)
        fw .= Mt[1] .* prim[1] .* Muv

    else
        tau =
            vhs_collision_time(prim, μᵣ, ω) +
            2.0 * dt * abs(primL[1] / primL[end] - primR[1] / primR[end]) /
            (primL[1] / primL[end] + primR[1] / primR[end])

        Mu, Mxi, MuL, MuR = gauss_moments(prim, inK)
        sw0L = (w .- wL) ./ dxL
        sw0R = (wR .- w) ./ dxR
        gaL = pdf_slope(prim, sw0L, inK)
        gaR = pdf_slope(prim, sw0R, inK)
        sw =
            -prim[1] .* (
                moments_conserve_slope(gaL, MuL, Mxi, 1) .+
                moments_conserve_slope(gaR, MuR, Mxi, 1)
            )
        # ga = pdf_slope(prim, sw, inK)
        # sw = -prim[1] .* moments_conserve_slope(ga, Mu, Mxi, 1)
        gaT = pdf_slope(prim, sw, inK)

        # time-integration constants
        Mt = zeros(5)
        Mt[4] = tau * (1.0 - exp(-dt / tau))
        Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]
        Mt[1] = dt - Mt[4]
        Mt[2] = -tau * Mt[1] + Mt[5]
        Mt[3] = 0.5 * dt^2 - tau * Mt[1]

        # flux related to central distribution
        Muv = moments_conserve(Mu, Mxi, 1, 0)
        MauL = moments_conserve_slope(gaL, MuL, Mxi, 2)
        MauR = moments_conserve_slope(gaR, MuR, Mxi, 2)
        # Mau = moments_conserve_slope(ga, MuR, Mxi, 2)
        MauT = moments_conserve_slope(gaT, Mu, Mxi, 1)

        fw .=
            Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* (MauL .+ MauR) .+
            Mt[3] .* prim[1] .* MauT
        # fw .= Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* Mau .+ Mt[3] .* prim[1] .* MauT
    end

    return nothing

end

"""
$(SIGNATURES)

2D
"""
function flux_equilibrium!(
    fw::T1,
    wL::T2,
    wR::T2,
    inK::Real,
    γ::Real,
    μᵣ::Real,
    ω::Real,
    Pr::Real,
    dt::Real,
    dxL::Real,
    dxR::Real,
    dy::Real,
    swL = zeros(eltype(fw), axes(wL))::T1,
    swR = zeros(eltype(fw), axes(wR))::T1,
) where {T1<:AA{<:FN,1},T2<:AA{<:Real,1}}

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

    w =
        primL[1] .* moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0) .+
        primR[1] .* moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)
    prim = conserve_prim(w, γ)

    # positivity check
    if prim[end] < 0.0
        primL = conserve_prim(wL .- swL .* dxL, γ)
        primR = conserve_prim(wR .+ swR .* dxR, γ)

        Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
        Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

        w =
            primL[1] .* moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0) .+
            primR[1] .* moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)
        prim = conserve_prim(w, γ)

        swL .= 0.0
        swR .= 0.0
    end

    tau =
        vhs_collision_time(prim, μᵣ, ω) +
        2.0 * dt * abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end])

    faL = pdf_slope(primL, swL, inK)
    sw = -primL[1] .* moments_conserve_slope(faL, Mu1, Mv1, Mxi1, 1, 0)
    faTL = pdf_slope(primL, sw, inK)

    faR = pdf_slope(primR, swR, inK)
    sw = -primR[1] .* moments_conserve_slope(faR, Mu2, Mv1, Mxi2, 1, 0)
    faTR = pdf_slope(primR, sw, inK)

    Mu, Mv, Mxi, MuL, MuR = gauss_moments(prim, inK)
    sw0L = (w .- wL) ./ dxL
    sw0R = (wR .- w) ./ dxR
    gaL = pdf_slope(prim, sw0L, inK)
    gaR = pdf_slope(prim, sw0R, inK)
    sw =
        -prim[1] .* (
            moments_conserve_slope(gaL, MuL, Mv, Mxi, 1, 0) .+
            moments_conserve_slope(gaR, MuR, Mv, Mxi, 1, 0)
        )
    # ga = pdf_slope(prim, sw, inK)
    # sw = -prim[1] .* moments_conserve_slope(ga, Mu, Mv, Mxi, 1, 0)
    gaT = pdf_slope(prim, sw, inK)

    # time-integration constants
    Mt = zeros(5)
    Mt[4] = tau * (1.0 - exp(-dt / tau))
    Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]
    Mt[1] = dt - Mt[4]
    Mt[2] = -tau * Mt[1] + Mt[5]
    Mt[3] = 0.5 * dt^2 - tau * Mt[1]

    # flux related to central distribution
    Muv = moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    MauL = moments_conserve_slope(gaL, MuL, Mv, Mxi, 2, 0)
    MauR = moments_conserve_slope(gaR, MuR, Mv, Mxi, 2, 0)
    # Mau = moments_conserve_slope(ga, Mu, Mv, Mxi, 2, 0)
    MauT = moments_conserve_slope(gaT, Mu, Mv, Mxi, 1, 0)

    @. fw =
        (Mt[1] * prim[1] * Muv + Mt[2] * prim[1] * (MauL + MauR) + Mt[3] * prim[1] * MauT) *
        dy
    # fw .= (Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* Mau .+ Mt[3] .* prim[1] .* MauT) .* dy

    return nothing

end


text:
# ------------------------------------------------------------------
# Copyright (c) 2018, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    pointgradients(img, points; [method])

Compute the gradients along all dimensions of N-dimensional `img`
at `points` using `method`. Default method is `:ando3`.
"""
function pointgradients(img::AbstractArray, points::AbstractVector; method=:ando3)
  extent = size(img)
  ndirs = length(extent)
  npoints = length(points)

  # smoothing weights
  weights = (method == :sobel ? [1,2,1] :
             method == :ando3 ? [.112737,.274526,.112737] :
             error("Unknown gradient method: $method"))

  # pad input image
  padimg = padarray(img, Pad(:replicate, ones(Int, ndirs), ones(Int, ndirs)))

  # gradient matrix
  G = zeros(npoints, ndirs)

  # compute gradient for all directions at specified points
  for i in 1:ndirs
    # kernel = centered difference + perpendicular smoothing
    if extent[i] > 1
      # centered difference
      idx = ones(Int, ndirs); idx[i] = 3
      kern = reshape([-1,0,1], idx...)
      # perpendicular smoothing
      for j in setdiff(1:ndirs, i)
        if extent[j] > 1
          idx = ones(Int, ndirs); idx[j] = 3
          kern = broadcast(*, kern, reshape(weights, idx...))
        end
      end

      A = zeros(size(kern))
      shape = size(kern)
      for (k, icenter) in enumerate(points)
        i1 = CartesianIndex(ntuple(i->1, ndirs))
        for ii in CartesianIndices(shape)
          A[ii] = padimg[ii + icenter - i1]
        end

        G[k,i] = sum(kern .* A)
      end
    end
  end

  G
end

ext: jl


text:
module BarotropicQGQL

export
  Problem,
  set_zeta!,
  updatevars!,

  energy,
  enstrophy,
  dissipation,
  work,
  drag

using
  FFTW,
  CUDA,
  Reexport

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!
using FourierFlows: parsevalsum, parsevalsum2
import FFTW: rfft

abstract type BarotropicQGQLVars <: AbstractVars end

nothingfunction(args...) = nothing

"""
    Problem(dev=CPU(); parameters...)

Construct a BarotropicQGQL turbulence problem on device `dev`.
"""
function Problem(dev::Device=CPU();
  # Numerical parameters
            nx = 256,
            Lx = 2π,
            ny = nx,
            Ly = Lx,
            dt = 0.01,
  # Physical parameters
             β = 0.0,
           eta = nothing,
  # Drag and/or hyper-/hypo-viscosity
             ν = 0.0,
            nν = 1,
             μ = 0.0,
  # Timestepper and equation options
       stepper = "RK4",
         calcF = nothingfunction,
    stochastic = false,
             T = Float64)

  # the grid
  grid = TwoDGrid(dev, nx, Lx, ny, Ly; T=T)
  x, y = gridpoints(grid)

  # topographic PV
  eta === nothing && ( eta = zeros(dev, T, (nx, ny)) )

  params = !(typeof(eta)<:ArrayType(dev)) ?
           Params(grid, β, eta, μ, ν, nν, calcF) :
           Params(β, eta, rfft(eta), μ, ν, nν, calcF)

  vars = calcF == nothingfunction ?
         Vars(dev, grid) : stochastic ?
                           StochasticForcedVars(dev, grid) :
                           ForcedVars(dev, grid)

  equation = BarotropicQGQL.Equation(params, grid)

  FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev)
end

# ----------
# Parameters
# ----------

"""
    Params

A struct that contains all parameter values for a two-dimensional barotropic QG QL problem.
"""
struct Params{T, Aphys, Atrans} <: AbstractParams
       β :: T          # Planetary vorticity y-gradient
     eta :: Aphys      # Topographic PV
    etah :: Atrans     # FFT of Topographic PV
       μ :: T          # Linear drag
       ν :: T          # Viscosity coefficient
      nν :: Int        # Hyperviscous order (nν=1 is plain old viscosity)
  calcF! :: Function   # Function that calculates the forcing on QGPV q
end

"""
    Params(grid, β, eta::Function, μ, ν, nν, calcF)

Constructor for `params` that accepts a generating function for the topographic PV, `eta`.
"""
function Params(grid::AbstractGrid{T, A}, β, eta::Function, μ, ν, nν, calcF) where {T, A}
  x, y = gridpoints(grid)
  eta_on_grid = A([eta(grid.x[i], grid.y[j]) for i=1:grid.nx, j=1:grid.ny])
  etah_on_grid = rfft(eta_on_grid)
  Params(β, eta_on_grid, etah_on_grid, μ, ν, nν, calcF)
end

# ---------
# Equations
# ---------

"""
    Equation(params, grid)

Returns the equation for two-dimensional barotropic QG QL problem with parameters `params` and on `grid`.
"""
function Equation(params::Params, grid::AbstractGrid)
  L = @. -params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr * grid.invKrsq
  CUDA.@allowscalar L[1, 1] = 0
  FourierFlows.Equation(L, calcN!, grid)
end

# ----
# Vars
# ----

"""
    Vars
A struct that contains all variables for a two-dimensional barotropic QG QL problem.
"""
mutable struct Vars{Aphys, Atrans, F, P} <: BarotropicQGQLVars
        u :: Aphys
        v :: Aphys
        U :: Aphys
    uzeta :: Aphys
    vzeta :: Aphys
     zeta :: Aphys
     Zeta :: Aphys
      psi :: Aphys
      Psi :: Aphys
       Nz :: Atrans
       NZ :: Atrans
       uh :: Atrans
       vh :: Atrans
       Uh :: Atrans
    zetah :: Atrans
    Zetah :: Atrans
     psih :: Atrans
     Psih :: Atrans
       Fh :: F
  prevsol :: P
end

const ForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, Nothing}
const StochasticForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, <:AbstractArray}

"""
    Vars(dev, grid)

Returns the vars for unforced two-dimensional barotropic QG problem with `grid`.
"""
function Vars(dev::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)
  @devzeros Dev T (grid.nx, grid.ny) u v U uzeta vzeta zeta Zeta psi Psi
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) N NZ uh vh Uh zetah Zetah psih Psih
  Vars(u, v, U, uzeta, vzeta, zeta, Zeta, psi, Psi, N, NZ, uh, vh, Uh, zetah, Zetah, psih, Psih, nothing, nothing)
end

"""
    ForcedVars(dev, grid)

Returns the `vars` for forced two-dimensional barotropic QG QL problem on device `dev` and `grid`.
"""
function ForcedVars(dev::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)
  @devzeros Dev T (grid.nx, grid.ny) u v U uzeta vzeta zeta Zeta psi Psi
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) N NZ uh vh Uh zetah Zetah psih Psih Fh
  Vars(u, v, U, uzeta, vzeta, zeta, Zeta, psi, Psi, N, NZ, uh, vh, Uh, zetah, Zetah, psih, Psih, Fh, nothing)
end

"""
    StochasticForcedVars(dev, grid)

Returns the `vars` for stochastically forced two-dimensional barotropic QG QL problem on
device `dev` and `grid`.
"""
function StochasticForcedVars(dev::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)
  @devzeros Dev T (grid.nx, grid.ny) u v U uzeta vzeta zeta Zeta psi Psi
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) N NZ uh vh Uh zetah Zetah psih Psih Fh prevsol
  Vars(u, v, U, uzeta, vzeta, zeta, Zeta, psi, Psi, N, NZ, uh, vh, Uh, zetah, Zetah, psih, Psih, Fh, prevsol)
end

# -------
# Solvers
# -------

function calcN_advection!(N, sol, t, clock, vars, params, grid)
  @. vars.zetah = sol
  CUDA.@allowscalar @. vars.zetah[1, :] = 0
  @. vars.Zetah = sol
  CUDA.@allowscalar @. vars.Zetah[2:end, :] = 0

  @. vars.uh =  im * grid.l  * grid.invKrsq * vars.zetah
  @. vars.vh = -im * grid.kr * grid.invKrsq * vars.zetah
  @. vars.Uh =  im * grid.l  * grid.invKrsq * vars.Zetah

  ldiv!(vars.zeta, grid.rfftplan, vars.zetah)
  ldiv!(vars.u, grid.rfftplan, vars.uh)
  ldiv!(vars.v, grid.rfftplan, vars.vh)

  ldiv!(vars.Zeta, grid.rfftplan, vars.Zetah)
  ldiv!(vars.U, grid.rfftplan, vars.Uh)

  @. vars.uzeta = vars.u * vars.zeta           # u*ζ
  @. vars.vzeta = vars.v * vars.zeta           # v*ζ

  mul!(vars.uh, grid.rfftplan, vars.uzeta)     # \hat{u*q}
  @. vars.NZ = -im * grid.kr * vars.uh         # -∂[u*q]/∂x
  mul!(vars.vh, grid.rfftplan, vars.vzeta)     # \hat{v*q}
  @. vars.NZ += - im * grid.l * vars.vh        # -∂[v*q]/∂y
  CUDA.@allowscalar @. vars.NZ[2:end, :] = 0

  @. vars.U = vars.U * vars.zeta                # U*ζ
  @. vars.u = vars.u * vars.Zeta                # u*Ζ
  @. vars.v = vars.v * vars.Zeta                # v*Ζ

  mul!(vars.uh, grid.rfftplan, vars.U + vars.u) # \hat{U*ζ + u*Ζ}
  @. vars.Nz = -im * grid.kr*vars.uh            # -∂[U*ζ + u*Ζ]/∂x
  mul!(vars.vh, grid.rfftplan, vars.v)          # \hat{v*Z}
  @. vars.Nz += - im * grid.l*vars.vh           # -∂[v*Z]/∂y
  CUDA.@allowscalar @. vars.Nz[1, :] = 0

  @. N = vars.NZ + vars.Nz

  return nothing
end

function calcN!(N, sol, t, clock, vars, params, grid)
  calcN_advection!(N, sol, t, clock, vars, params, grid)
  addforcing!(N, sol, t, clock, vars, params, grid)
  return nothing
end

addforcing!(N, sol, t, cl, vars::Vars, params, grid) = nothing

function addforcing!(N, sol, t, clock, vars::ForcedVars, params, grid)
  params.calcF!(vars.Fh, sol, t, clock, vars, params, grid)
  @. N += vars.Fh

  return nothing
end

function addforcing!(N, sol, t, clock, vars::StochasticForcedVars, params, grid)
  if t == clock.t # not a substep
    @. vars.prevsol = sol # sol at previous time-step is needed to compute budgets for stochastic forcing
    params.calcF!(vars.Fh, sol, t, clock, vars, params, grid)
  end
  @. N += vars.Fh

  return nothing
end

# ----------------
# Helper functions
# ----------------

"""
    updatevars!(sol, vars, params, grid)
    updatevars!(prob)

Update the `vars` of a problem `prob` that has `grid` and `params` with the solution in `sol`.
"""
function updatevars!(sol, vars, params, grid)
  CUDA.@allowscalar sol[1, 1] = 0
  @. vars.zetah = sol
  CUDA.@allowscalar @. vars.zetah[1, :] = 0
  @. vars.Zetah = sol
  CUDA.@allowscalar @. vars.Zetah[2:end, :] = 0

  @. vars.Psih = -vars.Zetah * grid.invKrsq
  @. vars.psih = -vars.zetah * grid.invKrsq
  @. vars.uh = -im * grid.l  * vars.psih
  @. vars.vh =  im * grid.kr * vars.psih
  @. vars.Uh =  im * grid.l  * vars.Zetah * grid.invKrsq

  ldiv!(vars.zeta, grid.rfftplan, deepcopy(vars.zetah))
  ldiv!(vars.Zeta, grid.rfftplan, deepcopy(vars.Zetah))
  ldiv!(vars.psi, grid.rfftplan, vars.psih)
  ldiv!(vars.Psi, grid.rfftplan, vars.Psih)
  ldiv!(vars.u, grid.rfftplan, deepcopy(vars.uh))
  ldiv!(vars.v, grid.rfftplan, deepcopy(vars.vh))
  ldiv!(vars.U, grid.rfftplan, deepcopy(vars.Uh))

  return nothing
end

updatevars!(prob) = updatevars!(prob.sol, prob.vars, prob.params, prob.grid)

"""
    set_zeta!(prob, zeta)
    set_zeta!(sol, vars, grid, zeta)

Set the solution `sol` as the transform of `zeta` and update variables `vars` on the `grid`.
"""
function set_zeta!(sol, vars, params, grid, zeta)
  mul!(vars.zetah, grid.rfftplan, zeta)
  CUDA.@allowscalar vars.zetah[1, 1] = 0.0
  @. sol = vars.zetah

  updatevars!(sol, vars, params, grid)
  return nothing
end

set_zeta!(prob, zeta) = set_zeta!(prob.sol, prob.vars, prob.params, prob.grid, zeta)

"""
    energy(sol, grid)
    energy(prob)

Returns the domain-averaged kinetic energy of `sol`.
"""
@inline energy(sol, grid::AbstractGrid) =
  0.5 * (parsevalsum2(grid.kr .* grid.invKrsq .* sol, grid)
        + parsevalsum2(grid.l .* grid.invKrsq .* sol, grid)) / (grid.Lx * grid.Ly)

energy(prob) = energy(prob.sol, prob.grid)

"""
    enstrophy(sol, grid, vars)
    enstrophy(prob)

Returns the domain-averaged enstrophy of `sol`.
"""
function enstrophy(sol, grid::AbstractGrid, vars::AbstractVars)
  @. vars.uh = sol
  CUDA.@allowscalar vars.uh[1, 1] = 0
  return 0.5 * parsevalsum2(vars.uh, grid) / (grid.Lx * grid.Ly)
end
enstrophy(prob) = enstrophy(prob.sol, prob.grid, prob.vars)

"""
    dissipation(prob)
    dissipation(sol, v, p, g)

Returns the domain-averaged dissipation rate. `nν` must be >= 1.
"""
@inline function dissipation(sol, vars, params, grid)
  @. vars.uh = grid.Krsq^(params.nν - 1) * abs(sol)^2
  CUDA.@allowscalar vars.uh[1, 1] = 0
  return params.ν / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid)
end

@inline dissipation(prob) = dissipation(prob.sol, prob.vars, prob.params, prob.grid)

"""
    work(prob)
    work(sol, vars, params, grid)

Returns the domain-averaged rate of work of energy by the forcing, `params.Fh`.
"""
@inline function work(sol, vars::ForcedVars, grid)
  @. vars.uh = grid.invKrsq * sol * conj(vars.Fh)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid)
end

@inline function work(sol, vars::StochasticForcedVars, grid)
  @. vars.uh = grid.invKrsq * (vars.prevsol + sol) / 2 * conj(vars.Fh) # Stratonovich
  # @. vars.uh = grid.invKrsq * vars.prevsol * conj(vars.Fh)             # Ito
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid)
end

@inline work(prob) = work(prob.sol, prob.vars, prob.grid)

"""
    drag(prob)
    drag(sol, vars, params, grid)

Returns the extraction of domain-averaged energy by drag `μ`.
"""
@inline function drag(prob)
  sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid
  @. vars.uh = grid.invKrsq * abs(sol)^2
  CUDA.@allowscalar vars.uh[1, 1] = 0
  return params.μ / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid)
end

end # module


text:
@coeff function C_B(reactor::Reactor)
  _gamma(value) = Gamma(value) / Gamma(nu_p+value)

  nu_p = nu_n + nu_T

  cur_C_B = _gamma(5/4)

  cur_C_B += 2*(1-3*nu_J)*_gamma(9/4)
  cur_C_B += (1-4*nu_J+9*nu_J^2)*_gamma(13/4)
  cur_C_B += 2*nu_J*(1-3*nu_J)*_gamma(17/4)
  cur_C_B += nu_J^2*_gamma(21/4)

  cur_C_B *= Gamma(nu_p) / ((1-nu_J)^2)
end


text:
"This abstract type needs to be implemented by all equations."
abstract type Equation end

"This abstract type needs to be implementred by all scenarios"
abstract type Scenario end

"""
This function initializes the discrete solution at t = 0 for an equation `eq` and a scenario `scenario`.
It evaluates get\\_initial\\_values for `eq` and `scenario` at the global coordinates of our nodes.
It then interpolates this on the `basis`.
"""
function interpolate_initial_dofs(eq::Equation, scenario::Scenario, celldofs, cell)
    # Note: It would be better to compute spatial average here
    # but then special care would be needed in some cases
    # as arithmetic mean is not always the correct "average"

    # Instead evaluate at center
    pxg, pyg = globalposition(cell, (0.5, 0.5))
    celldofs .= get_initial_values(eq, scenario, (pxg, pyg))
end

"This function returns the number of variables for equation `eq::Equation`"
function get_nvars(eq::Equation)
    get_ndofs(eq) - get_nparams(eq)
end

"This function returns the number of material parameters for equation `eq::Equation`."
function get_nparams(eq::Equation)
    0
end

"""

    get_variable_name(eq::Equation, variable_index)

Return the name of the variable with index `variable_index` for equation `eq::Equation`.
# Implementation
Custom equations need to define `get_variable_names(eq)` for this to work.
"""
function get_variable_name(eq::Equation, variable_index)
    get_variable_names(eq)[variable_index]
end

"""
    is_periodic_boundary(equation::Equation, scenario::Scenario)

Returns true if all boundary conditions should be periodic for equation `eq::Equation` and
scenario `scenario::Scenario`.
"""
function is_periodic_boundary(equation::Equation, scenario::Scenario)
    false
end

"""
    evaluate_boundary(eq::Equation, scenario::Equation, face, normalidx, dofsface, dofsfaceneigh)

Evaluate the boundary condition, defined for equation `eq` and `scenario`, on face `face` with `normalidx`.
Store dofs of ghost cell in `dofsfaceneigh`.
You can use inner (face)-dofs `dofsface` to construct to boundary.
"""
function evaluate_boundary(eq::Equation, scenario::Equation, face, normalidx, dofsface, dofsfaceneigh)
    error("evaluate_boundary called but not define.")
end

"""
    get_initial_values(eq::Equation, scenario::Scenario, global_position; t=0.0)

Return a vector of initial values for equation `eq` and scenario `scenario` at
`global_positon` and at time `t`.
Values other than `t = 0` only have to be supported if
`is_analytical_solution` returns true.
"""
function get_initial_values(eq::Equation, scenario::Scenario, global_position; t=0.0)
    error("Function get_initial_values not defined for equation/scenario")
end

"""
    is_analytical_solution(equation::Equation, scenario::Scenario)

Returns true if the initial condition is also an analytical solution
for equation `equation` and scenario `scenario`.
"""
function is_analytical_solution(equation::Equation, scenario::Scenario)
    false
end

"""
    evaluate_flux(eq::Equation, celldofs, cellflux)

Evaluate the flux for equation `eq` and dofs `celldofs`,
Stores them in `cellflux`.
"""
function evaluate_flux(eq::Equation, celldofs, cellflux)
    error("Function evaluate_flux not defined for equation/scenario")
end

"""
    max_eigenval(eq::Equation, celldata, normalidx)

Return the maximum eigenvalue in direction of `normalidx` for
`celldata`.

"""
function max_eigenval(eq::Equation, celldata, normalidx)
    error("Function max_eigenval not defined for equation/scenario")
end

"""
    declare_dofs(eq, dof_names, num_params=0)

Generate boilerplate code for equation type `eq` which has
degrees of freedom named `dof_names` (should be iterable of symbols)
and `num_params` parameters.
Generates

    get_variable_names(eq)
    get_ndofs(eq)
    get_nvars(eq)
    get_nparams(eq)

Also generates a struct called [typename]Shortcuts.

# Example

    struct Advection <: Equation end
    @declare_dofs Advection [:ρ_1, :ρ_2, :ρ_3]

Generates:

    get_variable_names(eq::Advection) = [:ρ_1, :ρ_2, :ρ_3]
    get_ndofs(eq::Advection) = 3
    get_nvars(eq::Advection) = 3
    get_nparams(eq::Advection) = 0

You can use the struct in this way:

    s = AdvectionShortcuts()
    @assert s.ρ_1 == 1
    @assert s.ρ_2 == 2
    @assert s.ρ_3 == 3
"""
macro declare_dofs(eq, dof_names, num_params=0)
    local expressions::Array{Expr, 1} = []
    local shortcut_struct = Symbol(eq, "Shortcuts")
    push!(expressions, quote
        function $(esc(Symbol("get_variable_names")))(eq::$eq)
            $dof_names
        end

        function $(esc(Symbol("get_ndofs")))(eq::$eq)
            $(length(dof_names.args))
        end

        function $(esc(Symbol("get_nvars")))(eq::$eq)
            $(length(dof_names.args) - num_params)
        end

        function $(esc(Symbol("get_nparams")))(eq::$eq)
            $num_params
        end

        struct $(esc(shortcut_struct))

        end

        @inline function Base.getproperty(x::$(esc(shortcut_struct)), name::Symbol)
            get_variable_shortcut(x, Val(name))
        end
    end
    )

    for (i, name) in enumerate(dof_names.args)
        push!(expressions, quote
            @inline function $(esc(Symbol("get_variable_shortcut")))(eq::$(esc(shortcut_struct)), ::Val{$name})
                return $i
            end

            @inline function $(esc(Symbol("get_variable_shortcut")))(eq::$eq, ::Val{$name})
                return $i
            end
        end)
        end

    quote $(expressions...) end
end

include("equations/advection.jl")
include("equations/telegraphers.jl")

"""
    make_equation(config::Configuration)

Returns equation object for configuration `config`.
"""
function make_equation(config::Configuration)
    if config.equation_name == "advection"
        return Advection()
    elseif config.equation_name == "telegraphers"
        return Telegraphers()
    else
        error(string("Unknown equation name: ", config.equation_name))
    end
end

"""
    make_scenario(config::Configuration)

Returns scenario object for configuration `config`.
"""
function make_scenario(config::Configuration)
    if config.scenario_name == "planar_waves"
        return PlanarWaves()
    else
        error(string("Unknown scenario name: ", scenario_name))
    end
end


text:
function _issubtype(a, b, store)
    _isany(b) && return true
    _type_compare(a, b) && return true
    sup_a = _super(a, store)
    _type_compare(sup_a, b) && return true
    !_isany(sup_a) && return _issubtype(sup_a, b, store)
    return false
end

_isany(x::SymbolServer.FakeTypeName) = x.name == VarRef(VarRef(nothing, :Core), :Any)
_isany(x::SymbolServer.DataTypeStore) = x.name.name == VarRef(VarRef(nothing, :Core), :Any)
_isany(x) = false

_type_compare(a::SymbolServer.DataTypeStore, b::SymbolServer.DataTypeStore) = a.name == b.name
_type_compare(a::SymbolServer.FakeTypeName, b::SymbolServer.FakeTypeName) = a == b
_type_compare(a::SymbolServer.FakeTypeName, b::SymbolServer.DataTypeStore) = a == b.name
_type_compare(a::SymbolServer.DataTypeStore, b::SymbolServer.FakeTypeName) = a.name == b
_type_compare(a::SymbolServer.DataTypeStore, b::SymbolServer.FakeUnion) = _type_compare(a, b.a) ||
_type_compare(a, b.b)

function _type_compare(a::SymbolServer.DataTypeStore, b::SymbolServer.FakeTypeVar)
    if b.ub isa SymbolServer.FakeUnion
        return _type_compare(a, b.ub)
    end
    a == b
end

_type_compare(a, b) = a == b

_super(a::SymbolServer.DataTypeStore, store) = SymbolServer._lookup(a.super.name, store)
_super(a::SymbolServer.FakeTypeVar, store) = a.ub
_super(a::SymbolServer.FakeUnionAll, store) = a.body
_super(a::SymbolServer.FakeTypeName, store) = _super(SymbolServer._lookup(a.name, store), store)
@static if !(Vararg isa Type)
    _super(a::SymbolServer.FakeTypeofVararg, store) = CoreTypes.Any
end

function _super(b::Binding, store)
    StaticLint.CoreTypes.isdatatype(b.type) || error()
    b.val isa Binding && return _super(b.val, store)
    sup = _super(b.val, store)
    if sup isa EXPR && StaticLint.hasref(sup)
        StaticLint.refof(sup)
    else
        store[:Core][:Any]
    end
end

function _super(x::EXPR, store)::Union{EXPR,Nothing}
    if x.head === :struct
        _super(x.args[2], store)
    elseif x.head === :abstract || x.head === :primtive
        _super(x.args[1], store)
    elseif CSTParser.issubtypedecl(x)
        x.args[2]
    elseif CSTParser.isbracketed(x)
        _super(x.args[1], store)
    end
end

function subtypes(T::Binding)
    @assert CSTParser.defines_abstract(T.val)
    subTs = []
    for r in T.refs
        if r isa EXPR && r.parent isa EXPR && CSTParser.issubtypedecl(r.parent) && r.parent.parent isa EXPR && CSTParser.defines_datatype(r.parent.parent)
            push!(subTs, r.parent.parent)
        end
    end
    subTs
end


text:
# Implementation of
#  "Table-driven Implementation of the Logarithm Function in IEEE Floating-point Arithmetic"
#  Tang, Ping-Tak Peter
#  ACM Trans. Math. Softw. (1990), 16(4):378--400
#  http://dx.doi.org/10.1145/98267.98294
# Copyright (c): Simon Byrne

# Does not currently handle floating point flags (inexact, div-by-zero, etc).

# nested methods inside let blocks cause problems (#18201)
# let
# global log_tang
# global log1p_tang

# Float64 lookup table.
# to generate values:
  # N=39 # (can be up to N=42).
  # sN = 2.0^N
  # isN = 1.0/sN
  # s7 = 2.0^7
  # is7 = 1.0/s7
  # for j=0:128
  #   l_big = Base.log(big(1.0+j*is7))
  #   l_hi = isN*Float64(round(sN*l_big))
  #   l_lo = Float64(l_big-l_hi)
  #   j % 2 == 0 && print("\n    ")
  #   print("(",l_hi,",",l_lo,"),")
  # end
const _log_tang_table_Float64 = [(0.0,0.0),(0.007782140442941454,-8.865052917267247e-13),
    (0.015504186536418274,-4.530198941364935e-13),(0.0231670592820592,-5.248209479295644e-13),
    (0.03077165866670839,4.529814257790929e-14),(0.0383188643027097,-5.730994833076631e-13),
    (0.04580953603181115,-5.16945692881222e-13),(0.053244514518155484,6.567993368985218e-13),
    (0.06062462181580486,6.299848199383311e-13),(0.06795066190898069,-4.729424109166329e-13),
    (0.07522342123775161,-1.6408301585598662e-13),(0.08244366921098845,8.614512936087814e-14),
    (0.08961215869021544,-5.283050530808144e-13),(0.09672962645890948,-3.5836667430094137e-13),
    (0.10379679368088546,7.581073923016376e-13),(0.11081436634049169,-2.0157368416016215e-13),
    (0.11778303565552051,8.629474042969438e-13),(0.1247034785010328,-7.556920687451337e-14),
    (0.1315763577895268,-8.075373495358435e-13),(0.13840232285838283,7.363043577087051e-13),
    (0.14518200984457508,-7.718001336828099e-14),(0.15191604202664166,-7.996871607743758e-13),
    (0.15860503017574956,8.890223439724663e-13),(0.16524957289584563,-5.384682618788232e-13),
    (0.17185025692742784,-7.686134224018169e-13),(0.17840765747314435,-3.2605717931058157e-13),
    (0.18492233849428885,-2.7685884310448306e-13),(0.1913948530000198,-3.903387893794952e-13),
    (0.1978257433293038,6.160755775588723e-13),(0.20421554142922105,-5.30156516006026e-13),
    (0.21056476910780475,-4.55112422774782e-13),(0.21687393830143264,-8.182853292737783e-13),
    (0.22314355131493357,-7.238189921749681e-13),(0.22937410106533207,-4.86240001538379e-13),
    (0.23556607131286,-9.30945949519689e-14),(0.24171993688651128,6.338907368997553e-13),
    (0.24783616390413954,4.4171755371315547e-13),(0.25391520998164196,-6.785208495970588e-13),
    (0.25995752443668607,2.3999540484211735e-13),(0.2659635484978935,-7.555569400283742e-13),
    (0.27193371548310097,5.407904186145515e-13),(0.2778684510030871,3.692037508208009e-13),
    (0.28376817313073843,-9.3834172236637e-14),(0.28963329258294834,9.43339818951269e-14),
    (0.29546421289342106,4.148131870425857e-13),(0.3012613305781997,-3.7923164802093147e-14),
    (0.3070250352957373,-8.25463138725004e-13),(0.31275571000333,5.668653582900739e-13),
    (0.318453731119007,-4.723727821986367e-13),(0.32411946865431673,-1.0475750058776541e-13),
    (0.32975328637257917,-1.1118671389559323e-13),(0.33535554192167183,-5.339989292003297e-13),
    (0.3409265869704541,1.3912841212197566e-13),(0.3464667673470103,-8.017372713972018e-13),
    (0.35197642315688427,2.9391859187648e-13),(0.3574558889213222,4.815896111723205e-13),
    (0.3629054936900502,-6.817539406325327e-13),(0.36832556115950865,-8.009990055432491e-13),
    (0.3737164097929053,6.787566823158706e-13),(0.37907835293481185,1.5761203773969435e-13),
    (0.3844116989112081,-8.760375990774874e-13),(0.38971675114044046,-4.152515806343612e-13),
    (0.3949938082405424,3.2655698896907146e-13),(0.40024316412745975,-4.4704265010452445e-13),
    (0.4054651081078191,3.452764795203977e-13),(0.4106599249844294,8.390050778518307e-13),
    (0.4158278951435932,1.1776978751369214e-13),(0.4209692946442374,-1.0774341461609579e-13),
    (0.42608439531068143,2.186334329321591e-13),(0.43117346481813,2.413263949133313e-13),
    (0.4362367667745275,3.90574622098307e-13),(0.44127456080423144,6.437879097373207e-13),
    (0.44628710262804816,3.713514191959202e-13),(0.45127464413963025,-1.7166921336082432e-13),
    (0.4562374334818742,-2.8658285157914353e-13),(0.4611757151214988,6.713692791384601e-13),
    (0.46608972992544295,-8.437281040871276e-13),(0.4709797152190731,-2.821014384618127e-13),
    (0.4758459048698569,1.0701931762114255e-13),(0.4806885293455707,1.8119346366441111e-13),
    (0.4855078157816024,9.840465278232627e-14),(0.49030398804461583,5.780031989454028e-13),
    (0.49507726679803454,-1.8302857356041668e-13),(0.4998278695566114,-1.620740015674495e-13),
    (0.5045560107519123,4.83033149495532e-13),(0.5092619017905236,-7.156055317238212e-13),
    (0.5139457511013461,8.882123951857185e-13),(0.5186077642083546,-3.0900580513238243e-13),
    (0.5232481437651586,-6.10765519728515e-13),(0.5278670896204858,3.565996966334783e-13),
    (0.532464798869114,3.5782396591276384e-13),(0.5370414658973459,-4.622608700154458e-13),
    (0.5415972824321216,6.227976291722515e-13),(0.5461324375974073,7.283894727206574e-13),
    (0.5506471179523942,2.680964661521167e-13),(0.5551415075406112,-1.0960825046059278e-13),
    (0.5596157879353996,2.3119493838005378e-14),(0.5640701382853877,-5.846905800529924e-13),
    (0.5685047353526897,-2.1037482511444942e-14),(0.5729197535620187,-2.332318294558741e-13),
    (0.5773153650352469,-4.2333694288141915e-13),(0.5816917396350618,-4.3933937969737843e-13),
    (0.5860490450031648,4.1341647073835564e-13),(0.590387446602108,6.841763641591467e-14),
    (0.5947071077462169,4.758553400443064e-13),(0.5990081896452466,8.367967867475769e-13),
    (0.6032908514389419,-8.576373464665864e-13),(0.6075552502243227,2.1913281229340092e-13),
    (0.6118015411066153,-6.224284253643115e-13),(0.6160298772156239,-1.098359432543843e-13),
    (0.6202404097512044,6.531043137763365e-13),(0.6244332880123693,-4.758019902171077e-13),
    (0.6286086594227527,-3.785425126545704e-13),(0.6327666695706284,4.0939233218678666e-13),
    (0.636907462236195,8.742438391485829e-13),(0.6410311794206791,2.521818845684288e-13),
    (0.6451379613736208,-3.6081313604225574e-14),(0.649227946625615,-5.05185559242809e-13),
    (0.6533012720119586,7.869940332335532e-13),(0.6573580727090302,-6.702087696194906e-13),
    (0.6613984822452039,1.6108575753932459e-13),(0.6654226325445052,5.852718843625151e-13),
    (0.6694306539429817,-3.5246757297904794e-13),(0.6734226752123504,-1.8372084495629058e-13),
    (0.6773988235909201,8.860668981349492e-13),(0.6813592248072382,6.64862680714687e-13),
    (0.6853040030982811,6.383161517064652e-13),(0.6892332812385575,2.5144230728376075e-13),
    (0.6931471805601177,-1.7239444525614835e-13)]

# Float32 lookup table
# to generate values:
  # N=16
  # sN = 2f0^N
  # isN = 1f0/sN
  # s7 = 2.0^7
  # is7 = 1.0/s7
  # for j=0:128
  #   j % 4 == 0 && print("\n    ")
  #   print(float64(Base.log(big(1.0+j*is7))),",")
  # end

const _log_tang_table_Float32 = [0.0,0.007782140442054949,0.015504186535965254,0.02316705928153438,
    0.030771658666753687,0.0383188643021366,0.0458095360312942,0.053244514518812285,
    0.06062462181643484,0.06795066190850775,0.07522342123758753,0.08244366921107459,
    0.08961215868968714,0.09672962645855111,0.10379679368164356,0.11081436634029011,
    0.11778303565638346,0.12470347850095724,0.13157635778871926,0.13840232285911913,
    0.1451820098444979,0.15191604202584197,0.15860503017663857,0.16524957289530717,
    0.17185025692665923,0.1784076574728183,0.184922338494012,0.19139485299962947,
    0.19782574332991987,0.2042155414286909,0.21056476910734964,0.21687393830061436,
    0.22314355131420976,0.22937410106484582,0.2355660713127669,0.24171993688714516,
    0.24783616390458127,0.25391520998096345,0.25995752443692605,0.26596354849713794,
    0.27193371548364176,0.2778684510034563,0.2837681731306446,0.28963329258304266,
    0.2954642128938359,0.3012613305781618,0.3070250352949119,0.3127557100038969,
    0.3184537311185346,0.324119468654212,0.329753286372468,0.3353555419211378,
    0.3409265869705932,0.34646676734620857,0.3519764231571782,0.3574558889218038,
    0.3629054936893685,0.3683255611587076,0.37371640979358406,0.37907835293496944,
    0.38441169891033206,0.3897167511400252,0.394993808240869,0.4002431641270127,
    0.4054651081081644,0.4106599249852684,0.415827895143711,0.42096929464412963,
    0.4260843953109001,0.4311734648183713,0.43623676677491807,0.4412745608048752,
    0.44628710262841953,0.45127464413945856,0.4562374334815876,0.46117571512217015,
    0.46608972992459924,0.470979715218791,0.4758459048699639,0.4806885293457519,
    0.4855078157817008,0.4903039880451938,0.4950772667978515,0.4998278695564493,
    0.5045560107523953,0.5092619017898079,0.5139457511022343,0.5186077642080457,
    0.5232481437645479,0.5278670896208424,0.5324647988694718,0.5370414658968836,
    0.5415972824327444,0.5461324375981357,0.5506471179526623,0.5551415075405016,
    0.5596157879354227,0.564070138284803,0.5685047353526688,0.5729197535617855,
    0.5773153650348236,0.5816917396346225,0.5860490450035782,0.5903874466021763,
    0.5947071077466928,0.5990081896460834,0.6032908514380843,0.6075552502245418,
    0.6118015411059929,0.616029877215514,0.6202404097518576,0.6244332880118935,
    0.6286086594223741,0.6327666695710378,0.6369074622370692,0.6410311794209312,
    0.6451379613735847,0.6492279466251099,0.6533012720127457,0.65735807270836,
    0.661398482245365,0.6654226325450905,0.6694306539426292,0.6734226752121667,
    0.6773988235918061,0.6813592248079031,0.6853040030989194,0.689233281238809,
    0.6931471805599453]

# Procedure 1
@inline function _log_tang_proc1(y::Float64,mf::Float64,F::Float64,f::Float64,jp::Int)
    ## Steps 1 and 2
    @inbounds hi,lo = _log_tang_table_Float64[jp]
    l_hi = mf* 0.6931471805601177 + hi
    l_lo = mf*-1.7239444525614835e-13 + lo

    ## Step 3
    # @inbounds u = f*c_invF[jp]
    # u = f/F
    # q = u*u*@horner(u,
    #                 -0x1.0_0000_0000_0001p-1,
    #                 +0x1.5_5555_5550_9ba5p-2,
    #                 -0x1.f_ffff_ffeb_6526p-3,
    #                 +0x1.9_99b4_dfed_6fe4p-3,
    #                 -0x1.5_5576_6647_2e04p-3)

    ## Step 3' (alternative)
    u = (2*f)/(y+F)
    v = u*u
    q = u*v*@horner(v,
                    0.08333333333303913,
                    0.012500053168098584)

    ## Step 4
    l_hi + (u + (q + l_lo))
end

# Procedure 2
@inline function _log_tang_proc2(f::Float64)
    ## Step 1
    g = 1/(2+f)
    u = 2*f*g
    v = u*u

    ## Step 2
    q = u*v*@horner(v,
                    0.08333333333333179,
                    0.012500000003771751,
                    0.0022321399879194482,
                    0.0004348877777076146)

    ## Step 3
    # based on:
    #   2(f-u) = 2(f(2+f)-2f)/(2+f) = 2f^2/(2+f) = fu
    #   2(f-u1-u2) - f*(u1+u2) = 0
    #   2(f-u1) - f*u1 = (2+f)u2
    #   u2 = (2(f-u1) - f*u1)/(2+f)
    if is_fma_fast(Float64)
        return u + fma(fma(-u,f,2(f-u)), g, q)
    else
        u1 = trunclo(u)
        f1 = trunclo(f)
        f2 = f-f1
        u2 = ((2*(f-u1)-u1*f1)-u1*f2)*g
        ## Step 4
        return u1 + (u2 + q)
    end
end

@inline function _log_tang_proc1(y::Float32,mf::Float32,F::Float32,f::Float32,jp::Int)
    ## Steps 1 and 2
    @inbounds hi = _log_tang_table_Float32[jp]
    l = mf*0.6931471805599453 + hi

    ## Step 3
    # @inbounds u = f*c_invF[jp]
    # q = u*u*@horner(u,
    #                 Float32(-0x1.00006p-1),
    #                 Float32(0x1.55546cp-2))

    ## Step 3' (alternative)
    u = (2*f)/(y+F)
    v = u*u
    q = u*v*0.08333351f0

    ## Step 4
    Float32(l + (u + q))
end

@inline function _log_tang_proc2(f::Float32)
    ## Step 1
    # compute in higher precision
    u64 = Float64(2*f)/(2+Float64(f))
    u = Float32(u64)
    v = u*u

    ## Step 2
    q = u*v*@horner(v,
                    0.08333332f0,
                    0.012512346f0)

    ## Step 3: not required

    ## Step 4
    Float32(u64 + q)
end

function log_tang(x::T) where T<:FloatTypes

    if x > 0
        x == T(Inf) && return x

        # Step 2
        # exp(-1/16) < x < exp(1/16)
        if T(0.9394130628134757) < x < T(1.0644944589178595)
            f = x-1
            return _log_tang_proc2(f)
        end

        # Step 3
        xu = reinterpret(Unsigned, x)
        m = (xu >> significand_bits(T))%Int
        if m == 0 # x is subnormal
            x *= T(2)^significand_bits(T) # normalise significand
            xu = reinterpret(Unsigned,x)
            m = (xu >> significand_bits(T))%Int - significand_bits(T)
        end
        m -= exponent_bias(T)
        y = reinterpret(T,(xu & significand_mask(T)) | exponent_one(T))

        mf = T(m)

        k = T(2)^(significand_bits(T) - 7)
        F = (y + k) - k # round(y,7,2)
        f = y-F
        jp = unsafe_trunc(Int, T(2)^7 * F) - 127

        return _log_tang_proc1(y,mf,F,f,jp)
    elseif x == 0
        -T(Inf)
    elseif isnan(x)
        T(NaN)
    else
        throw(DomainError())
    end
end

function log1p_tang(x::T) where T<:FloatTypes
    if x > -1
        x == Inf && return x
        if -eps(T)/2 < x < eps(T)/2
            return x # Inexact
        end

        # Step 2
        # expm1(-1/16) < x < expm1(1/16)
        if T(-0.06058693718652422) < x < T(0.06449445891785943)
            return _log_tang_proc2(x)
        end

        # Step 3
        z = 1 + x # z > 0, and not subnormal

        zu = reinterpret(Unsigned, z)
        m = (zu >> significand_bits(T))%Int - exponent_bias(T)           # exponent(z)
        s = reinterpret(T, (exponent_one(T)<<1) - (zu & exponent_mask(T))) # 2^-m (except for largest normal binade, where it gives 0, but that doesn't matter)
        c = m > 0 ? 1-(z-x) : x-(z-1)                                    # 1+x = z+c exactly
        y = reinterpret(T,(zu & significand_mask(T)) | exponent_one(T))  # y = 2 ∈ [1,2)
        mf = T(m)

        k = T(2)^(significand_bits(T) - 7)
        F = (y + k) - k # round(y,7,2)
        f = (y - F) + c*s # F+f = (1+x)*2^-m = (z+c)**2^-m
        jp = unsafe_trunc(Int, T(2)^7 * F) - 127

        _log_tang_proc1(y,mf,F,f,jp)
    elseif x == -1
        -T(Inf)
    elseif isnan(x)
        T(NaN)
    else
        throw(DomainError())
    end
end
#end


text:
@enum Direction FFT_FORWARD=-1 FFT_BACKWARD=1
@enum Pow24 POW2=2 POW4=1

abstract type AbstractFFTType end

# Represents a Composite Cooley-Tukey FFT
struct CompositeFFT <: AbstractFFTType end

# Represents a Radix-2 Cooley-Tukey FFT
struct Pow2FFT <: AbstractFFTType end

# Represents a Radix-3 Cooley-Tukey FFT
struct Pow3FFT <: AbstractFFTType end

# Represents a Radix-4 Cooley-Tukey FFT
struct Pow4FFT <: AbstractFFTType end

# Represents an O(N²) DFT
struct DFT <: AbstractFFTType end

"""
$(TYPEDSIGNATURES)
Node of a call graph

# Arguments
`left`: Offset to the left child node
`right`: Offset to the right child node
`type`: Object representing the type of FFT
`sz`: Size of this FFT

"""
struct CallGraphNode
    left::Int
    right::Int
    type::AbstractFFTType
    sz::Int
    s_in::Int
    s_out::Int
end

"""
$(TYPEDSIGNATURES)
Object representing a graph of FFT Calls

# Arguments
`nodes`: Nodes keeping track of the graph
`workspace`: Preallocated Workspace

"""
struct CallGraph{T<:Complex}
    nodes::Vector{CallGraphNode}
    workspace::Vector{Vector{T}}
end

# Get the node in the graph at index i
Base.getindex(g::CallGraph{T}, i::Int) where {T} = g.nodes[i]

"""
$(TYPEDSIGNATURES)
Check if `N` is a power of `base`

"""
function _ispow(N, base)
    while N % base == 0
        N = N/base
    end
    return N == 1
end

"""
$(TYPEDSIGNATURES)
Check if `N` is a power of 2 or 4

"""
function _ispow24(N::Int)
    N < 1 && return nothing
    while N & 0b11 == 0
        N >>= 2
    end
    return N < 3 ? Pow24(N) : nothing
end

"""
$(TYPEDSIGNATURES)
Recursively instantiate a set of `CallGraphNode`s

# Arguments
`nodes`: A vector (which gets expanded) of `CallGraphNode`s
`N`: The size of the FFT
`workspace`: A vector (which gets expanded) of preallocated workspaces
`s_in`: The stride of the input
`s_out`: The stride of the output

"""
function CallGraphNode!(nodes::Vector{CallGraphNode}, N::Int, workspace::Vector{Vector{T}}, s_in::Int, s_out::Int)::Int where {T}
    if iseven(N)
        pow = _ispow24(N)
        if !isnothing(pow)
            push!(workspace, T[])
            push!(nodes, CallGraphNode(0, 0, pow == POW2 ? Pow2FFT() : Pow4FFT(), N, s_in, s_out))
            return 1
        end
    end
    if N % 3 == 0
        if _ispow(N, 3)
            push!(workspace, T[])
            push!(nodes, CallGraphNode(0, 0, Pow3FFT(), N, s_in, s_out))
            return 1
        end
    end
    if isprime(N)
        push!(workspace, T[])
        push!(nodes, CallGraphNode(0,0, DFT(),N, s_in, s_out))
        return 1
    end
    Ns = [first(x) for x in collect(factor(N)) for _ in 1:last(x)]
    if Ns[1] == 2
        N1 = prod(Ns[Ns .== 2])
    elseif Ns[1] == 3
        N1 = prod(Ns[Ns .== 3])
    else
        # Greedy search for closest factor of N to sqrt(N)
        Nsqrt = sqrt(N)
        N_cp = cumprod(Ns[end:-1:1])[end:-1:1]
        N_prox = abs.(N_cp .- Nsqrt)
        _,N1_idx = findmin(N_prox)
        N1 = N_cp[N1_idx]
    end
    N2 = N ÷ N1
    push!(nodes, CallGraphNode(0,0,DFT(),N,s_in,s_out))
    sz = length(nodes)
    push!(workspace, Vector{T}(undef, N))
    left_len = CallGraphNode!(nodes, N1, workspace, N2, N2*s_out)
    right_len = CallGraphNode!(nodes, N2, workspace, N1*s_in, 1)
    nodes[sz] = CallGraphNode(1, 1 + left_len, CompositeFFT(), N, s_in, s_out)
    return 1 + left_len + right_len
end

"""
$(TYPEDSIGNATURES)
Instantiate a CallGraph from a number `N`

"""
function CallGraph{T}(N::Int) where {T}
    nodes = CallGraphNode[]
    workspace = Vector{Vector{T}}()
    CallGraphNode!(nodes, N, workspace, 1, 1)
    CallGraph(nodes, workspace)
end


text:
using Metatheory

taylor = @theory x a b begin
  exp(x) --> Σ(x^:n / factorial(big(:n)))
  cos(x) --> Σ((-1)^:n * x^2(:n) / factorial(big(2 * :n)))
  Σ(a) + Σ(b) --> Σ(a + b)
end

function expand(iters)
  RewriteRule(
    PatTerm(:call, :Σ, [PatVar(:a)], @__MODULE__),
    PatTerm(:call, :sum, [PatTerm(:(->), :(->), [:n, PatVar(:a)], @__MODULE__), 0:iters], @__MODULE__),
  )
end

a = rewrite(:(exp(x) + cos(x)), taylor)

r = expand(5000)
bexpr = rewrite(a, [r])

# you may want to do algebraic simplification
# with egraphs here

x = big(42)

b = eval(bexpr)
# 1.739274941520501044994695988622883932193276720547806372656638132701531037200611e+18

exp(x) + cos(x)
# 1.739274941520501046994695988622883932193276720547806372656638132701531037200651e+18

@testset "Infinite Series Approximation" begin
  @test b ≈ (exp(x) + cos(x))
end


text:
#Посчитать число всех горизонтальных прямолинейных перегородок (вертикальных - нет)

function ugol(r,side1,side2)
    while !isborder(r,side1) || !isborder(r,side2)
        if !isborder(r,side1)
            move!(r,side1)
        else
            move!(r,side2)
        end
    end
end

function kray(r,side) #функция перемещает робота к краю
    while isborder(r,side) == false
        move!(r,side)
    end
end

function num_vert(r)
    ugol(r, West, Sud)
    n = 0
    flag = 1
    while !isborder(r, Nord)
        while !isborder(r, Ost)
            move!(r, Ost)
            if(flag == 0 && isborder(r, Nord))
                flag = 1
                n += 1
            elseif(flag == 1 && !isborder(r, Nord))
                flag = 0
            end
        end
        move!(r, Nord)
        kray(r,West)
    end
print(n)
end

max_stars_repo_name: gardend/mirea-progs

max_issues_repo_name: gardend/mirea-progs

max_forks_repo_name: gardend/mirea-progs


text:
#
# Dimension: translates sequences of Symbols, Strings, or Ints to
# ordinal values, so the Dimension keys can be used for lookup in
# Arrays in Mimi. A Dimension can be declared one of several ways:
# dim = Dimension(:foo, :bar, :baz)      # varargs
# dim = Dimension([:foo, :bar, :baz])    # Vector
# dim = Dimension(2010:2100)             # Range
# dim = Dimension(4)                     # Same as 1:4
#
# Similarly, values can be referenced several ways:
# value  = Dimension[:foo]               # single lookup
# values = Dimension[(:foo, :bar)]       # Tuple of values
# values = Dimension[[2010, 2020, 2030]] # Vector of values
# values = Dimension[2010, 2020, 2030]   # Varargs
# values = Dimension[2010:10:2030]       # AbstractRange
#
using DataStructures

key_type(dim::Dimension{T}) where {T <: DimensionKeyTypes} = T

#
# Iteration and basic dictionary methods are delegated to the internal dict
#
Base.length(dim::Dimension) = length(dim.dict)
Base.iterate(dim::Dimension, state...) = iterate(dim.dict, state...)

Base.keys(dim::Dimension)   = keys(dim.dict)
Base.values(dim::Dimension) = values(dim.dict)

Base.get(dim::Dimension, key::Union{Number, Symbol, String}, default::Any) = get(dim.dict, key, default)

Base.getindex(dim::Dimension, key::Colon) = collect(values(dim.dict))

Base.getindex(dim::AbstractDimension, key::Union{Number, Symbol, String}) = getindex(dim.dict, key)

# Support dim[[:foo, :bar, :baz]], dim[(:foo, :bar, :baz)], and dim[2010:2020]
Base.getindex(dim::AbstractDimension, keys::Union{Vector{T} where T, Tuple, AbstractRange}) = [getindex(dim.dict, key) for key in keys]

Base.getindex(dim::AbstractDimension, keys...) = getindex(dim, keys)

Base.length(dim::RangeDimension) = length(dim.range)
Base.iterate(dim::RangeDimension, state...)  = iterate(dim.range, state...)

Base.keys(dim::RangeDimension)   = collect(dim.range)
Base.values(dim::RangeDimension) = collect(1:length(dim.range))

# Get first/last value from OrderedDict of keys
Base.firstindex(dim::AbstractDimension) = dim.dict.keys[1]
Base.lastindex(dim::AbstractDimension)  = dim.dict.keys[length(dim)]

#
# Compute the index of a "key" (e.g., a year) in the range.
#
function Base.get(dim::RangeDimension, key::Int, default::Any=0)
    r = dim.range
    i = key - r.start
    return i == 0 ? 1 : (i % r.step != 0 ? default : 1 + div(i, r.step))
end

# Support dim[[2010, 2020, 2030]], dim[(:foo, :bar, :baz)], and dim[2010:2050]
Base.getindex(dim::RangeDimension, keys::Union{Vector{Int}, Tuple, AbstractRange}) = [get(dim, key, 0) for key in keys]

# Symbols are added to the dim_dict in @defcomp (with value of nothing), but are set later using set_dimension!
has_dim(obj::AbstractCompositeComponentDef, name::Symbol) = (haskey(obj.dim_dict, name) && obj.dim_dict[name] !== nothing)

isuniform(obj::AbstractCompositeComponentDef) = obj.is_uniform

set_uniform!(obj::AbstractCompositeComponentDef, value::Bool) = (obj.is_uniform = value)

dimension(obj::AbstractCompositeComponentDef, name::Symbol) = obj.dim_dict[name]

dim_names(obj::AbstractCompositeComponentDef, dims::Vector{Symbol}) = [dimension(obj, dim) for dim in dims]

dim_count_dict(obj::AbstractCompositeComponentDef) = Dict([name => length(value) for (name, value) in dim_dict(obj)])

dim_counts(obj::AbstractCompositeComponentDef, dims::Vector{Symbol}) = [length(dim) for dim in dim_names(obj, dims)]
dim_count(obj::AbstractCompositeComponentDef, name::Symbol) = length(dimension(obj, name))

dim_keys(obj::AbstractCompositeComponentDef, name::Symbol)   = collect(keys(dimension(obj, name)))
dim_values(obj::AbstractCompositeComponentDef, name::Symbol) = collect(values(dimension(obj, name)))

"""
    set_dimension!(ccd::CompositeComponentDef, name::Symbol, keys::Union{Int, Vector, Tuple, AbstractRange})

Set the values of `ccd` dimension `name` to integers 1 through `count`, if `keys` is
an integer; or to the values in the vector or range if `keys` is either of those types.
"""
function set_dimension!(ccd::AbstractCompositeComponentDef, name::Symbol, keys::Union{Int, Vector, Tuple, AbstractRange})

    redefined = has_dim(ccd, name)
    dim = Dimension(keys)

    if name == :time

        # if we are redefining the time dimension of a model, the timestep length
        # must match the timestep length of the old time dimension
        redefined && _check_time_redefinition(ccd, keys)

        # propagate the time dimension through all sub-components
        _propagate_time_dim!(ccd, dim)
        set_uniform!(ccd, isuniform(keys))

        # if we are redefining the time dimension for a Model Definition
        # pad the time arrays with missings and update their time labels
        redefined && (ccd isa ModelDef) && _pad_parameters!(ccd)

    end

    return set_dimension!(ccd, name, dim)
end

"""
    set_dimension!(obj::AbstractComponentDef, name::Symbol, dim::Dimension)

Set the dimension `name` in `obj` to `dim`.
"""
function set_dimension!(obj::AbstractComponentDef, name::Symbol, dim::Dimension)
    dirty!(obj)
    obj.dim_dict[name] = dim

    if name == :time
        for subcomp in compdefs(obj)
            set_dimension!(subcomp, :time, dim)
        end
    end
    return dim
end

"""
    function add_dimension!(comp::AbstractComponentDef, name)

Add a dimension of name `name` to `comp`, where the dimension will be `nothing`
unless `name` is an Int in which case we create an "anonymouse" dimension on
the fly with keys `1` through `count` where `count` = `name`.
"""
function add_dimension!(comp::AbstractComponentDef, name)
    # generally, we add dimension name with nothing instead of a Dimension instance,
    # but in the case of an Int name, we create the "anonymous" dimension on the fly.
    dim = (name isa Int) ? Dimension(name) : nothing
    comp.dim_dict[Symbol(name)] = dim                         # TBD: test this
end

# Note that this operates on the registered comp, not one added to a composite
add_dimension!(comp_id::ComponentId, name) = add_dimension!(compdef(comp_id), name)

"""
    dim_names(ccd::AbstractCompositeComponentDef)

Return a list of the dimension names of `ccd`.
"""
function dim_names(ccd::AbstractCompositeComponentDef)
    dims = OrderedSet{Symbol}()             # use a set to eliminate duplicates
    for cd in compdefs(ccd)
        union!(dims, keys(dim_dict(cd)))    # TBD: test this
    end

    return collect(dims)
end

"""
    dim_names(comp_def::AbstractComponentDef, datum_name::Symbol)

Return a list of the dimension names of datum `datum_name` in `comp_def`.
"""
dim_names(comp_def::AbstractComponentDef, datum_name::Symbol) = dim_names(datumdef(comp_def, datum_name))

"""
    dim_count(def::AbstractDatumDef)

Return number of dimensions in `def`.
"""
dim_count(def::AbstractDatumDef) = length(dim_names(def))

"""
    _check_time_redefinition(obj::AbstractCompositeComponentDef, keys::Union{Int, Vector, Tuple, AbstractRange})

Run through all necesssary safety checks for redefining `obj`'s time dimenson to
a new dimension with keys `keys`.
"""
function _check_time_redefinition(obj::AbstractCompositeComponentDef, keys::Union{Int, Vector, Tuple, AbstractRange}) where T

    # get useful variables
    curr_keys = time_labels(obj)
    curr_first = obj.first
    curr_last = obj.last

    new_keys = [keys...]
    new_first = first(new_keys)
    new_last = last(new_keys)

    # (1) check that the shift is legal
    isa(obj, ModelDef) ? obj_name = "model" : obj_name = "component $(nameof(obj))"
    new_first > curr_first && error("Cannot redefine the time dimension to start at $new_first because it is after the $obj_name's current start $curr_first.")
    curr_first > new_last && error("Cannot redefine the time dimension to end at $new_last because it is before the $obj_name's current start $curr_first")

    # (2) check first and last
    !(curr_first in new_keys) && error("The current first index ($curr_first) must exist within the model's new time dimension $new_keys.") # can be assumed since we cannot move the time forward
    curr_last >= new_last && !(new_last in curr_keys) && error("The new last index ($new_last) must exist within the model's current time dimension $curr_keys, since the time redefinition contracts to an earlier year.")
    curr_last < new_last && !(curr_last in new_keys) && error("The current last index ($curr_last) must exist within the model's redefined time dimension $new_keys, since the time redefinition expands to a later year.")

    # (3) check that the overlap region between the current keys and new keys holds same keys
    if length(curr_keys) > 1 && length(new_keys) > 1
        if isuniform(curr_keys) # fixed timesteps
            step_size(curr_keys) != step_size(new_keys) && error("Cannot redefine the time dimension to have a timestep size of $(step_size(new_keys)), must match the timestep size of current time dimension, $(step_size(curr_keys))")

        else # variable timesteps
            start_idx = 1 # can be assumed since we cannot move the time forward
            new_last < curr_last ? end_idx = findfirst(isequal(new_last), curr_keys) : end_idx = length(curr_keys)
            expected_overlap = curr_keys[start_idx:end_idx]

            start_idx = findfirst(isequal(curr_first), new_keys)
            end_idx = start_idx + length(expected_overlap) - 1
            observed_overlap = new_keys[start_idx:end_idx]

            expected_overlap != observed_overlap && error("Cannot redefine the time dimension, the overlapping portion of the current and new times must be identical.")
        end
    end

end


text:
using Genie.Renderer.Html, Genie.Requests

greeting = "Welcome"
name = "Genie"

function htmlviewfile()
  "

$greeting

This is a $name test

  "
end

function htmltemplatefile()
  "
  <!DOCTYPE HTML>

    <title>$name test</title>

    <% @yield %>

    Just a footer

  "
end

@testset "HTML Rendering" begin
  @testset "WebRenderable constructors" begin
    wr = Genie.Renderer.WebRenderable("hello")
    @test wr.body == "hello"
    @test wr.content_type == Genie.Renderer.DEFAULT_CONTENT_TYPE
    @test wr.status == 200
    @test wr.headers == Genie.Renderer.HTTPHeaders()

    wr = Genie.Renderer.WebRenderable("hello", :json)
    @test wr.body == "hello"
    @test wr.content_type == :json
    @test wr.status == 200
    @test wr.headers == Genie.Renderer.HTTPHeaders()

    wr = Genie.Renderer.WebRenderable()
    @test wr.body == ""
    @test wr.content_type == Genie.Renderer.DEFAULT_CONTENT_TYPE
    @test wr.status == 200
    @test wr.headers == Genie.Renderer.HTTPHeaders()

    wr = Genie.Renderer.WebRenderable(body = "bye", content_type = :javascript, status = 301, headers = Genie.Renderer.HTTPHeaders("Location" => "/bye"))
    @test wr.body == "bye"
    @test wr.content_type == :javascript
    @test wr.status == 301
    @test wr.headers["Location"] == "/bye"

    wr = Genie.Renderer.WebRenderable(Genie.Renderer.WebRenderable(body = "good morning", content_type = :javascript), 302, Genie.Renderer.HTTPHeaders("Location" => "/morning"))
    @test wr.body == "good morning"
    @test wr.content_type == :javascript
    @test wr.status == 302
    @test wr.headers["Location"] == "/morning"
  end;

  @testset "String HTML rendering" begin
    r = Requests.HTTP.Response()

    @testset "String no layout" begin
      rm("build", force = true, recursive = true)
      r = html(htmlviewfile(), forceparse = true)

      @test String(r.body) == "
$greeting

This is a $name test

"
    end;

    @testset "String with layout" begin
      rm("build", force = true, recursive = true)
      r = html(htmlviewfile(), layout = htmltemplatefile())

      @test String(r.body) == "<title>$name test</title>

$greeting

This is a $name test

\n
Just a footer"
    end;

    @test r.status == 200
    @test r.headers[1]["Content-Type"] == "text/html; charset=utf-8"

    rm("build", force = true, recursive = true)
  end;
end;


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -29.23432407510253 <= q <= 25.005662503891426)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.9584800931585375 + 0.8291526833328144*1/(1 + exp(-(0.45833851557375827 + -0.07893784992449904*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.05269389376863742*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.9157122134504139*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.14234137821989812*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))))))))) + -0.47308336785041494*1/(1 + exp(-(0.8182392255084516 + 0.5838043616725095*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.4648893036822681*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.819336300807536*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.1444146117738878*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))))))))) + 0.6108443940953072*1/(1 + exp(-(-0.7015763069538381 + -0.6574215519391942*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.13767395548269734*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.21605537563661414*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.3991258949611094*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))))))))) + -0.4993930565410927*1/(1 + exp(-(-0.21931284624450997 + 0.9567749529797225*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.8630566689068395*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.8170122079846647*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.11471605851894884*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))))))))) + 1/(1 + exp(-(0.475745519197857 + -0.20728788527253217*1/(1 + exp(-(0.45833851557375827 + -0.07893784992449904*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.05269389376863742*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.9157122134504139*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.14234137821989812*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))))))))) + 0.5536591664216424*1/(1 + exp(-(0.8182392255084516 + 0.5838043616725095*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.4648893036822681*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.819336300807536*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.1444146117738878*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))))))))) + 0.5772469178059145*1/(1 + exp(-(-0.7015763069538381 + -0.6574215519391942*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.13767395548269734*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.21605537563661414*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.3991258949611094*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))))))))) + 0.14258534816037427*1/(1 + exp(-(-0.21931284624450997 + 0.9567749529797225*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.8630566689068395*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.8170122079846647*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.11471605851894884*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))))))))) + 1/(1 + exp(-(-0.19885854025662475 + 0.9363781836864602*1/(1 + exp(-(0.45833851557375827 + -0.07893784992449904*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.05269389376863742*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.9157122134504139*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.14234137821989812*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))))))))) + -0.2868856241387201*1/(1 + exp(-(0.8182392255084516 + 0.5838043616725095*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.4648893036822681*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.819336300807536*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.1444146117738878*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))))))))) + 0.6989952318088708*1/(1 + exp(-(-0.7015763069538381 + -0.6574215519391942*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.13767395548269734*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.21605537563661414*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.3991258949611094*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))))))))) + 0.2430276118759056*1/(1 + exp(-(-0.21931284624450997 + 0.9567749529797225*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.8630566689068395*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.8170122079846647*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.11471605851894884*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))))))))) + 1/(1 + exp(-(-0.4748268724301301 + -0.943475605090133*1/(1 + exp(-(0.45833851557375827 + -0.07893784992449904*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.05269389376863742*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.9157122134504139*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.14234137821989812*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))))))))) + -0.504486324487766*1/(1 + exp(-(0.8182392255084516 + 0.5838043616725095*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.4648893036822681*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.819336300807536*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.1444146117738878*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))))))))) + -0.8621266685661562*1/(1 + exp(-(-0.7015763069538381 + -0.6574215519391942*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.13767395548269734*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.21605537563661414*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.3991258949611094*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))))))))) + -0.8772206856506037*1/(1 + exp(-(-0.21931284624450997 + 0.9567749529797225*1/(1 + exp(-(-0.22994750449881307 + -0.3590609593879259*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.10337862983418544*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.321850031280575*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.612120130564894*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.8630566689068395*1/(1 + exp(-(0.19384674006188307 + -0.9239566552941776*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + -0.4571782009479768*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + -0.7401514915717242*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.3444562251116623*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + 0.8170122079846647*1/(1 + exp(-(0.4065836733363821 + 0.7200826994240335*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.11568748742914314*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.5958961905666036*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + -0.0402986153245215*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))) + -0.11471605851894884*1/(1 + exp(-(-0.18406387526761447 + 0.28290146144347617*1/(1 + exp(-(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])))) + 0.5991916639069861*1/(1 + exp(-(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])))) + 0.7799932719867426*1/(1 + exp(-(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])))) + 0.858692422726886*1/(1 + exp(-(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
export outascii_scalar,outascii_mesh

function outascii_mesh(fname::String,g::Geometry,t::AbstractTopology)

  topo_ids      = Dict(VertexTopology=>1,IntervalTopology=>2,TriangleTopology=>3,TetrahedronTopology=>4)
  topo_id       = topo_ids[typeof(t)]
  element_len   = [1;2;3;4][topo_id]

  fid = open(fname*".node","w");
  xx=zeros(Float64,3)
  for k=1:nvertices(g)
    for d=1:dim(g)
      xx[d] = xyz(g,k,d)
    end
    println(fid,"$(xx[1]) $(xx[2]) $(xx[3])")
  end
  close(fid)

  fid = open(fname*".ele","w");
  for e in t
    println(fid,"$(join(e[2:(1+element_len)]," "))")
  end
  close(fid)

end

function outascii_scalar(fname::String,g::Geometry,t::AbstractTopology,sname::String,u::Vector{Float64})

  topo_ids      = Dict(VertexTopology=>1,IntervalTopology=>2,TriangleTopology=>3,TetrahedronTopology=>4)
  topo_id       = topo_ids[typeof(t)]
  element_len   = [1;2;3;4][topo_id]

  # fid = open(fname*".node","w");
  # xx=zeros(Float64,3)
  # for k=1:nvertices(g)
  #   for d=1:dim(g)
  #     xx[d] = xyz(g,k,d)
  #   end
  #   println(fid,"$(xx[1]) $(xx[2]) $(xx[3])")
  # end
  # close(fid)
  #
  # fid = open(fname*".ele","w");
  # for e in t
  #   println(fid,"$(join(e[2:(1+element_len)]," "))")
  # end
  # close(fid)

  fid = open(fname*".sol","w");
  for i=1:length(u)
    println(fid,"$(u[i])")
    #  @printf(fid,"%4.3f\n",u[i])
  end
  close(fid)

end

max_stars_repo_name: dpeschka/jPDE

max_issues_repo_name: dpeschka/jPDE

max_forks_repo_name: dpeschka/jPDE