text:
"""
$(SIGNATURES)

Returns the number of arguments of `f` for the method which has the most arguments.
"""
function numargs(f)
  typ = Tuple{Any, Val{:analytic}, Vararg}
  typ2 = Tuple{Any, Type{Val{:analytic}}, Vararg} # This one is required for overloaded types
  typ3 = Tuple{Any, Val{:jac}, Vararg}
  typ4 = Tuple{Any, Type{Val{:jac}}, Vararg} # This one is required for overloaded types
  typ5 = Tuple{Any, Val{:tgrad}, Vararg}
  typ6 = Tuple{Any, Type{Val{:tgrad}}, Vararg} # This one is required for overloaded types
  numparam = maximum([(m.sig<:typ || m.sig<:typ2 || m.sig<:typ3 || m.sig<:typ4 || m.sig<:typ5 || m.sig<:typ6) ? 0 : num_types_in_tuple(m.sig) for m in methods(f)])
  return (numparam-1) #-1 in v0.5 since it adds f as the first parameter
end

"""
$(SIGNATURES)

Get the number of parameters of a Tuple type, i.e. the number of fields.
"""
function num_types_in_tuple(sig)
  length(sig.parameters)
end

function num_types_in_tuple(sig::UnionAll)
  length(Base.unwrap_unionall(sig).parameters)
end

"""
    isinplace(f, inplace_param_number)
    isinplace(f::AbstractDiffEqFunction[, inplace_param_number])

Check whether a function operates in place by comparing its number of arguments
to the expected number. The second parameter is optional if `f` is an
[`AbstractDiffEqFunction`](@ref DiffEqBase.AbstractDiffEqFunction).

# See also
* [`numargs`](@ref DiffEqBase.numargs)
"""
function isinplace(f,inplace_param_number)
  numargs(f)>=inplace_param_number
end

isinplace(f::AbstractDiffEqFunction{iip}) where {iip} = iip
isinplace(f::AbstractDiffEqFunction{iip}, inplace_param_number) where {iip} = iip

"""
    @CSI_str cmd

Create an ANSI escape sequence string for the CSI command `cmd`.
"""
macro CSI_str(str)
    return :(string("\x1b[", $(esc(str)), "m"))
end

const TYPE_COLOR = CSI"36"  # Cyan
const NO_COLOR = CSI"0"

"""
    @def name definition
"""
macro def(name, definition)
    return quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

using Compat.TypeUtils: typename

Base.@pure parameterless_type(T::Type) = typename(T).wrapper
Base.@pure parameterless_type(x) = parameterless_type(typeof(x))

# support functions
export check_keywords, warn_compat
function check_keywords(alg, kwargs, warnlist)
    flg = false
    for (kw, val) in kwargs
        if kw in warnlist
            if val !== nothing
                flg = true
                @warn(string("The ", kw, " argument is ignored by ", alg, "."))
            end
        end
    end
    flg
end

"""
$(SIGNATURES)

Emit a warning with a link to the solver compatibility chart in the documenation.
"""
warn_compat() =
    @warn("Please see http://docs.juliadiffeq.org/latest/basics/compatibility_chart.html")

"""
    @add_kwonly function_definition

Define keyword-only version of the `function_definition`.

    @add_kwonly function f(x; y=1)
        ...
    end

expands to:

    function f(x; y=1)
        ...
    end
    function f(; x = error("No argument x"), y=1)
        ...
    end
"""
macro add_kwonly(ex)
  esc(add_kwonly(ex))
end

add_kwonly(ex::Expr) = add_kwonly(Val{ex.head}, ex)

function add_kwonly(::Type{<: Val}, ex)
  error("add_only does not work with expression $(ex.head)")
end

function add_kwonly(::Union{Type{Val{:function}},
                            Type{Val{:(=)}}}, ex::Expr)
  body = ex.args[2:end]  # function body
  default_call = ex.args[1]  # e.g., :(f(a, b=2; c=3))
  kwonly_call = add_kwonly(default_call)
  if kwonly_call === nothing
    return ex
  end

  return quote
    begin
      $ex
      $(Expr(ex.head, kwonly_call, body...))
    end
  end
end

function add_kwonly(::Type{Val{:where}}, ex::Expr)
  default_call = ex.args[1]
  rest = ex.args[2:end]
  kwonly_call = add_kwonly(default_call)
  if kwonly_call === nothing
    return nothing
  end
  return Expr(:where, kwonly_call, rest...)
end

function add_kwonly(::Type{Val{:call}}, default_call::Expr)
  # default_call is, e.g., :(f(a, b=2; c=3))
  funcname = default_call.args[1]  # e.g., :f
  required = []  # required positional arguments; e.g., [:a]
  optional = []  # optional positional arguments; e.g., [:(b=2)]
  default_kwargs = []
  for arg in default_call.args[2:end]
    if isa(arg, Symbol)
      push!(required, arg)
    elseif arg.head == :(::)
      push!(required, arg)
    elseif arg.head == :kw
      push!(optional, arg)
    elseif arg.head == :parameters
      @assert default_kwargs == []  # can I have :parameters twice?
      default_kwargs = arg.args
    else
      error("Not expecting to see: $arg")
    end
  end
  if isempty(required) && isempty(optional)
    # If the function is already keyword-only, do nothing:
    return nothing
  end
  if isempty(required)
    # It's not clear what should be done.  Let's not support it at
    # the moment:
    error("At least one positional mandatory argument is required.")
  end

  kwonly_kwargs = Expr(:parameters, [
    Expr(:kw, pa, :(error($("No argument $pa"))))
    for pa in required
  ]..., optional..., default_kwargs...)
  kwonly_call = Expr(:call, funcname, kwonly_kwargs)
  # e.g., :(f(; a=error(...), b=error(...), c=1, d=2))

  return kwonly_call
end

"""
$(SIGNATURES)

List symbols `export`'ed but not actually defined.
"""
function undefined_exports(mod)
  undefined = []
  for name in names(mod)
    if ! isdefined(mod, name)
      push!(undefined, name)
    end
  end
  return undefined
end

function wrapfun_oop(ff, inputs::Tuple)
  IT = map(typeof, inputs)
  FunctionWrapper{IT[1], Tuple{IT...}}((args...)->(ff(args...)))
end

function wrapfun_iip(ff, inputs::Tuple)
  IT = map(typeof, inputs)
  FunctionWrapper{Nothing, Tuple{IT...}}((args...)->(ff(args...); nothing))
end

function unwrap_fw(fw::FunctionWrapper)
  fw.obj[]
end
_vec(v) = vec(v)
_vec(v::Number) = v
_vec(v::AbstractVector) = v

_reshape(v, siz) = reshape(v, siz)
_reshape(v::Number, siz) = v

macro tight_loop_macros(ex)
   :($(esc(ex)))
end

# Overloaded in other repositories
function unwrap_cache end


text:
using ForwardDiff, Distributions, FiniteDifferences
using Tracker, Zygote, ReverseDiff, Random, LinearAlgebra, PDMats
using DistributionsAD, Test, LinearAlgebra, Combinatorics
using ForwardDiff: Dual
using StatsFuns: binomlogpdf, logsumexp
const FDM = FiniteDifferences
using DistributionsAD: TuringMvNormal, TuringMvLogNormal, TuringUniform
using Distributions: meanlogdet

include("test_utils.jl")
include("distributions.jl")
include("others.jl")

ext: jl


text:
using CuArrays
using Flux
using Random

a = CuArray([ 1 2 3
    4 5 6f0])

na = FewShotAnomalyDetection.normalizecolumns(a)

b = param(a)
FewShotAnomalyDetection.normalizecolumns(b)

c = Array(a)
FewShotAnomalyDetection.normalizecolumns(c)

a = CuArray(randn(Float32, 1000, 1000))
b = param(a)
c = Array(a)

@elapsed FewShotAnomalyDetection.normalizecolumns(a)
@elapsed FewShotAnomalyDetection.normalizecolumns(b)
@elapsed FewShotAnomalyDetection.normalizecolumns(c)

max_stars_repo_name: honzabim/FewShotAnomalyDetection

max_issues_repo_name: honzabim/FewShotAnomalyDetection

max_forks_repo_name: honzabim/FewShotAnomalyDetection


text:
"""
    ds_loadPkg(path_RLibrary)

Function to load DataSHIELD-Client-Packages needed for Boosting.

The functions uses the package RCall.
The packages need to be installed seperately into the user R-library
and the corresponding server-side packages on the OPAL-Server
where DataSHIELD is running.

Needed DataSHIELD-Client-Packages:
opal, dsBaseClient, dsStatsClient, dsModellingClient

# Arguments
- `path_RLibrary`::String: Path to R library where R packages opal, dsBaseClient, dsStatsClient and dsModellingClient are saved.

# Returned objects
- `RCall.RObject{RCall.StrSxp}`: Overview of loaded R-packages.

# Examples
```julia-repl
julia> ds_loadPkg("C:/Users/Username/Documents/R/win-library/R-Version")
RCall.RObject{RCall.StrSxp}
 [1] "dsModellingClient" "dsStatsClient"     "dsBaseClient"
 [4] "fields"            "maps"              "spam"
 [7] "grid"              "dotCall64"         "opal"
[10] "mime"              "rjson"             "RCurl"
[13] "bitops"            "stats"             "graphics"
[16] "grDevices"         "utils"             "datasets"
[19] "methods"           "base"
```
"""
function ds_loadPkg(path_RLibrary::String)
    # using RCall
    R"""
    .libPaths($path_RLibrary)

    library(opal)
    library(dsBaseClient)
    library(dsStatsClient)
    library(dsModellingClient)
    """
end

"""
    ds_login(url, user, password, table, servernames)

Function to login into DataSHIELD-Server.

The function uses the package RCall.
The DataSHIELD function datashield.login() is used.
The function ds_loadPkg() needs to be called first.

# Arguments
- `url::Array{String,1}`: Array containing URLs to all DataSHIELD-Server as Strings.
- `user::String`: String containing the User-Name for login
- `password::String`: String containin either the password or the name of the private key file.
- `table::Array{String,1}`: Array containing the table names (either one dimensional if always the same or of the same length as the URL array)
- `servernames::Array{String,1}`: Array containing the server names (same dimension as url).

# Returned objects
- `RCall.RObject{RCall.StrSxp}`: Overview of loaded DataSHIELD-server.
- Trace of DataSHIELD (e.g. loaded variables).

# Examples
```julia-repl
julia> ds_login(["https://server_url1:port", "https://server_url2:port"],
"user","password",["Projectname.Tablename"],["Server1","Server2"])
RCall.RObject{RCall.VecSxp}
Server1
url: https://server_url1:port
name: Server1
version: x.x.x.-xxxx
username: user
rid: xxxx-xxxx-xxx

Server2
url: https://server_url2:port
name: Server2
version: x.x.x.-xxxx
username: user
rid: xxxx-xxxx-xxx
```
"""
function ds_login(url::Array{String,1},user::String,password::String,table::Array{String,1},servernames::Array{String,1})
    # Compliance with DataSHIELD rules
    if(length(table)>1 && length(table)!=length(url))
        error("table must be either one-dimensional or of the same dimension as url")
    end
    if(length(servernames)!=length(url))
        error("servernames must be of the same dimension as url")
    end

    # Login into DataSHIELD using datashield.login(). All variables are assigned.
    R"""
    my_logindata <- data.frame(server=$servernames, url=$url, user=$user, password=$password, table=$table)

    opal_login <- datashield.login(logins=my_logindata, assign=TRUE)
    """
end

"""
    ds_check(serveranz, ignore=false, tolerance=0.0000000000000005, labels=Array{String,1}())

Function to check if the variables are standardized.

The function uses the package RCall.
The DataSHIELD functions ds.colnames(), ds.dim(), ds.mean(), ds.var().
Firstly, the mean will be checked (should be 0)
and only if the check is valid, the sd will be checked (should be 1).
The check is interrupted when the first unstandardized variable is found.
The function ds_login() needs to be called first.

# Arguments
- `serveranz::Int`: Number of used servers.
- `tolerance::Float64=0.0000000005`: Tolerance for comparison of mean and sd.
- `ignore::Bool=false`: If true: the algorithm will continue even if the variables are not standardizes, if false the algorithm will stop. The default is false.
- `labels::Array{String,1}=Array{String,1}()`: Labels of variables that should be checked. By default all variables are checked.

# Returned objects
- No return objects.
- Trace includes overview of already checked variables.

# Examples
```julia-repl
julia> ds_check(2)
```
"""
function ds_check(serveranz::Int, ignore::Bool=false, tolerance::Float64=0.0000000005, labels::Array{String,1}=Array{String,1}())
    if(isempty(labels))
		# Get all names of the assigned table D
        names_cov = rcopy(R"ds.colnames('D')[[1]]")
    else
        names_cov = deepcopy(labels)
    end

    # Check if mean approx 0, loop is stopped if the first is not
    mean_values = R"""
    tolerance <-  $tolerance
    dimension <- length($names_cov)
    print($names_cov)
    res <- FALSE
    cat("Check of mean: ")
    for (i in 1:dimension) {
        cat($names_cov[i],", ")
        if(any(abs(unlist(ds.mean(paste("D",$names_cov[i],sep="$"), type = "split"))) > tolerance)) {
            res <- TRUE
            break
        }
    }
    res
    """

    if(rcopy(mean_values))
        # mean unqual to 0, thus no standarization was performed
        if(ignore)
            @warn("Values need to be standardized, but mean value is unqual to 0 (variance is not checked and check for mean is stopped). The results of the algorithm might be unvalid.")
        else
            error("Values need to be standardized, but mean value is unqual to 0 (variance is not checked and check for mean is stopped). The results of the algorithm might be unvalid and the evaluation has been stopped.")
        end
    else
        # mean approx. 0, so sd is checked for approx 1, loop is stopped if the first is not
        sd_values = R"""
        tolerance <-  $tolerance
        dimension <- length($names_cov)
        res <- FALSE
        cat("\n Check of sd: ")
        for (i in 1:dimension) {
            cat($names_cov[i],", ")
            if(any(abs(unlist(ds.var(paste("D",$names_cov[i],sep="$"), type = "split"))-1) > tolerance)) {
                res <- TRUE
                break
            }
        }
        res
        """

        if(rcopy(sd_values))
            # sd unqual to 1, thus no standardization was performed
            if(ignore)
                @warn("Values need to be standardized, but sd value is unqual to 1 (mean is equal to 0 but check of variance is stopped). The results of the algorithm might be unvalid.")
            else
                error("Values need to be standardized, but sd value is unqual to 1 (mean is equal to 0 but check of variance is stopped). The results of the algorithm might be unvalid and the evaluation has been stopped.")
            end
        end
    end
end

"""
    ds_start(path_RLibrary, url, user, password, table, servernames, check=true, ignore=false, labels::Array{String,1}=Array{String,1}())

Function to wrap the complete login process.

Calls the functions ds_loadPKG(), ds_login, and optionally ds_check().

# Arguments
- `path_RLibrary::String`: Path to R library where R packages opal, dsBaseClient, dsStatsClient and dsModellingClient are saved
- `url::Array{String,1}`: Array containing URLs to all DataSHIELD-Server as Strings.
- `user::String`: String containing the User-Name for login
- `password::String`: String containin either the password or the name of the private key file.
- `table::Array{String,1}`: Array containing the table names (either one dimensional if always the same or of the same length as the URL array)
- `servernames::Array{String,1}`: Array containing the server names (same dimension as url).
- `check::Bool=true`: If true: It is checked if the variables are standardized, if false not. The default is true.
- `ignore::Bool=false`: Only relevant if check=true. If true: the algorithm will continue even if the variables are not standardizes, if false the algorithm will stop. The default is false.
- `labels::Array{String,1}=Array{String,1}()`: Labels of variables that should be checked.

# Returned objects
- No objects are returned.
- Trace of DataSHIELD (e.g. loaded variables).

# Examples
```julia-repl
julia> ds_start("C:/Users/Username/Documents/R/win-library/R-Version",["https://server_url1:port", "https://server_url2:port"],
"user","password",["Projectname.Tablename"],["Server1","Server2"],false)
```
"""
function ds_start(path_RLibrary::String, url::Array{String,1},user::String,password::String,table::Array{String,1},servernames::Array{String,1}, check::Bool=true, ignore::Bool=false, tolerance::Float64=0.0000000005, labels::Array{String,1}=Array{String,1}())
    # Load needed R-package
    ds_loadPkg(path_RLibrary)

    # Login to DataSHIELD
    ds_login(url, user, password, table, servernames)

    if(check)
        # Checking of standardization
        @warn("The standardization of the variables will be checked (mean=0 (firstly), sd=1). All variables are checked subsequently and the algorithm stops when at least one variable is not standardized. The algrotihm might take some time.")
        ds_check(length(servernames),ignore,tolerance,labels)
    end
end

"""
    ds_logout()

Function to log out of DataSHIELD, where the connection is made with ds_login().

# Arguments
No arguments.

# Returned objects
- `RCall.RObject{RCall.StrSxp}`: Showing the logout server.

# Examples
```julia-repl
julia> ds_logout()
RCall.RObject{RCall.VecSxp}
Server1
NULL

Server2
NULL
```
"""
function ds_logout()
    R"""
    datashield.logout(opal_login)
    """
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 944bd174-5e89-42c0-aef0-a3c4450e8194
begin

	# Setup the Julia environment -
	using PlutoUI
	using MAT

	# setup paths -
	const _PATH_TO_NOTEBOOK = pwd()
	const _PATH_TO_DATA = joinpath(_PATH_TO_NOTEBOOK,"data")
	const _PATH_TO_FIGS = joinpath(_PATH_TO_NOTEBOOK,"figs")

	# return -
	nothing
end

# ╔═╡ 3a51e7a0-48e6-4b02-8649-9708ccb6ffa5
md"""
## Introduction to Metabolic Engineering and Flux Balance Analysis

Metabolic engineering is the practice of optimizing genetic and regulatory processes within cells to increase the cell's production of a desired molecule or protein of interest.

Metabolic engineers manipulate the biochemical networks used by cells to convert raw materials into molecules necessary for the cell's survival. Metabolic engineering specifically seeks to:

1. Mathematically model biochemical networks, calculate the yield (product divided substrate) of useful products, and identify parts of the network that constrain the production of the products of interest.
1. Use genetic engineering techniques to modify the biochemical network in order to relieve constraints limiting production. The modified network can then be modeled to calculate the new product yield, and to identify new constraints (back to 1).

Resources for biochemical network information:
* [Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000 Jan 1;28(1):27-30. doi: 10.1093/nar/28.1.27. PMID: 10592173; PMCID: PMC102409.](https://www.genome.jp/kegg/)

* [Karp, Peter D et al. “The BioCyc collection of microbial genomes and metabolic pathways.” Briefings in bioinformatics vol. 20,4 (2019): 1085-1093. doi:10.1093/bib/bbx085](https://pubmed.ncbi.nlm.nih.gov/29447345/)

* [Gama-Castro, Socorro et al. “RegulonDB version 9.0: high-level integration of gene regulation, coexpression, motif clustering and beyond.” Nucleic acids research vol. 44,D1 (2016): D133-43. doi:10.1093/nar/gkv1156](https://pubmed.ncbi.nlm.nih.gov/26527724/)
"""

# ╔═╡ 7299f66a-2ca4-4a53-b5f2-bf249fe436f3
md"""
__Fig 1.__ The overall metabolic map from the KEGG database. Each dot (_node_) is a metabolite, each line (_edge_) is a metabolic reaction.
"""

# ╔═╡ 328bbf10-c6ec-4d3e-8ffb-0a09c327920e
PlutoUI.LocalResource("./figs/KEGG-map01100.png")

# ╔═╡ 3642b338-d819-4c0b-a93f-59006ce04174
md"""
## What is Flux Balance Analysis (FBA)?
Flux balance analysis (FBA) is a mathematical modeling and analysis approach that computes the flux (units: $\star$mol/gDW-time or $\star$mol/L-time for cell-free networks) of carbon and energy throughout a metabolic network. FBA, a member of the [constraint-based family of mathematical modeling tools](https://pubmed.ncbi.nlm.nih.gov/26000478/), is a widely used approach to compute metabolic flux. However, there are alternatives to FBA, such as metabolic flux analysis (MFA), but these alternatives vary more in the solution approach than the structure of the estimation problem.

Let's look at the following reference to better understand the different components of a flux balance analysis problem:

* [Orth, J., Thiele, I. & Palsson, B. What is flux balance analysis?. Nat Biotechnol 28, 245–248 (2010). https://doi.org/10.1038/nbt.1614](https://www.ncbi.nlm.nih.gov/labs/pmc/articles/PMC3108565/)

Computational tools for working with flux balance analysis models:

* [Heirendt, Laurent et al. “Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0.” Nature protocols vol. 14,3 (2019): 639-702. doi:10.1038/s41596-018-0098-2](https://pubmed.ncbi.nlm.nih.gov/30787451/)

#### Flux balance analysis problem structure
The FBA problem is typically encoded as a linear programming (LP) problem of the form:

$$\max_{v}\sum_{i=1}^{\mathcal{R}}c_{i}v_{i}$$

subject to the constraints:

$$\begin{eqnarray}
\mathbf{S}\mathbf{v} &=& \mathbf{0}\\
\mathcal{L}\leq&\mathbf{v}&\leq\mathcal{U}\\
\dots &=& \dots
\end{eqnarray}$$

where $\mathbf{S}$ denotes the stoichiometric matrix, $c_{i}$ denote the objective coefficients,
$\mathbf{v}$ denotes the metabolic flux (the unknown that we are trying to estimate), and
$\mathcal{L}$ ($\mathcal{U}$) denote the permissible lower (upper) bounds on the _unkown_ metabolic flux. The first set of constraints enforces conservation of mass, while the second imparts thermodynamic and kinetic information into the calculation. Finally, there are potentially other types of constraints (both linear and nonlinear) that can be used in this type of problem (we will not cover these here, but these additional constraints may be important in specific applications).

### What is a stoichiometric matrix?
The basis for flux balance calculations is the stoichiometric matrix.
The stochiometric matrix is a digitial representation of the biochemistry that can occur insides cells (such as that shown in Fig. 1). The stochiometric matrix, denoted by $\mathbf{S}$, is a $\mathcal{M}\times\mathcal{R}$ array of stoichiometric coefficients, denoted by the symbol $\sigma_{ij}$ ($i$ is the row index, $j$ is the col index). Each of the $\mathcal{M}$ rows of $\mathbf{S}$ describes a particular metabolite (_node_ in the metabolic network) while each of the $\mathcal{R}$ columns corresponds to a metabolic reaction (_edge_ in the metabolic network). Thus:

* A stoichiometric coefficient $\sigma_{ij}$ > 0 implies that metabolite $i$ is __produced__ by reaction $j$
* A stoichiometric coefficient $\sigma_{ij}$  = 0 implies that metabolite $i$ is __not connected__ to reaction $j$
* A stoichiometric coefficient $\sigma_{ij}$ < 0 implies that metabolite $i$ is __consumed__ by reaction $j$

"""

# ╔═╡ d92b38f1-1604-4b4d-8288-4c94f58d9869
md"""

__Example__: Core stoichiometric matrix from _Escherichia coli_: [Orth, Jeffrey D et al. “Reconstruction and Use of Microbial Metabolic Networks: the Core Escherichia coli Metabolic Model as an Educational Guide.” EcoSal Plus vol. 4,1 (2010): 10.1128/ecosalplus.10.2.1. doi:10.1128/ecosalplus.10.2.1](https://pubmed.ncbi.nlm.nih.gov/26443778/)
"""

# ╔═╡ 4b426b24-f416-4079-9951-275164d792f3
begin

	# what is the model name -
	model_file_name = "modelReg.mat"
	model_name = "modelReg"

	# where is model file?
	_PATH_TO_MODEL_FILE = joinpath(_PATH_TO_DATA,model_file_name)

	# load the mat file -> get the cobra_dictionary
	file = matopen(_PATH_TO_MODEL_FILE)
    cobra_dictionary = read(file, model_name)
    close(file)
end

# ╔═╡ 14ccf90a-7e77-4112-914f-7b81917389ed
keys(cobra_dictionary)

# ╔═╡ 1edb1823-001d-4e16-bc66-01f8136fbf1c
# get the stoichiometric matrix
stm_sparse = cobra_dictionary["S"]

# ╔═╡ ec8ed359-94e2-4a3a-a08c-386bd4166d00
stm_full = Matrix(stm_sparse)

# ╔═╡ ca63dac5-7fa4-498e-af3e-52c917be3dc2
rxns = cobra_dictionary["rxns"]

# ╔═╡ a5b8c0b0-fd80-47b9-ad52-6267e50b0539
mets = cobra_dictionary["metFormulas"]

# ╔═╡ 5d81d380-df70-4480-aaa6-6bde2ce60826
md"""
Additional resources for the stoichiometric matrix

* [Lecture 3: The genome reconstruction process. Systems Biology: Constraint-based Reconstruction and Analysis, Cambridge University Press, 2015](https://www.youtube.com/watch?v=jcRVSJj3XP8&feature=youtu.be)

* [Lecture 9: Properties of the Stoichiometric Matrix. Systems Biology: Constraint-based Reconstruction and Analysis, Cambridge University Press, 2015](https://www.youtube.com/watch?v=0nsjKx6a7X4&feature=youtu.be)

"""

# ╔═╡ 181c2a31-8196-4f00-87aa-b828144d1963
md"""
### What are flux bounds constraints?
Flux bounds constraints control the permissible ranges that a metabolic reaction rate (flux) can have. Flux bounds constraints can incorporate thermodynamic information (whether a reaction is reversible) as well as kinetic information.

Resources for thermodynamics in the constraint-based world:

* [Peres, Sabine, and Vincent Fromion. “Thermodynamic Approaches in Flux Analysis.” Methods in molecular biology (Clifton, N.J.) vol. 2088 (2020): 359-367. doi:10.1007/978-1-0716-0159-4_17](https://pubmed.ncbi.nlm.nih.gov/31893383/)

* [Flamholz, Avi et al. “eQuilibrator--the biochemical thermodynamics calculator.” Nucleic acids research vol. 40,Database issue (2012): D770-5. doi:10.1093/nar/gkr874](https://pubmed.ncbi.nlm.nih.gov/22064852/)

Let's do some sample calculations using [eQuilibrator](https://equilibrator.weizmann.ac.il)

Thermodynamic information tells you the directionality of the reaction (is the reaction reversible), but to get the permissible magnitudes of the reaction we need to compute the kinetics. We (and others) formulate flux bounds as the product of succesive corrections t the maximum possible rate e.g., for an irrereversible reaction:

$$0\leq{v_{i}}\leq{V_{max,j}^{\circ}}\left(\frac{\epsilon}{\epsilon^{\circ}}\right)\theta_{j}\left(\dots\right){f_{j}\left(\dots\right)}$$

where $V_{max,j}^{\circ}$ denotes the maximum reaction velocity computed for enzyme $j$ at some characteristic enzyme concetration (and full activity), the ratio $\epsilon/\epsilon^{\circ}$ is a correction for enzyme concentration, $\theta_{j}\left(\dots\right)\in\left[0,1\right]$ is an enzyme activity function (or measurement) and $f_{j}\left(\dots\right)$ is a function describing the substrate dependence of the reaction rate $j$. Both $\theta_{j}\left(\dots\right)$ and $f_{j}\left(\dots\right)$ could have associated parameters, e.g., saturation or binding constants, etc.

"""

# ╔═╡ 4b7dcaec-f7a0-44ad-898e-fbca9f5271e5
md"""
### What is the objective of an _E. coli_ cell?

__Short answer__: maximize the growth rate

__Longer answer__: This question was raging in the metabolic engineering community twenty years ago, but is mostly forgotten today. First, there are long-standing arguments from ecology (or even from the perspective of microeconomics) as to why _E.coli_ would be goal-oriented. However, Palsson and coworkers largely put this question to rest (at least for _E.coli_) with the publication:

* [Ibarra, Rafael U et al. “Escherichia coli K-12 undergoes adaptive evolution to achieve in silico predicted optimal growth.” Nature vol. 420,6912 (2002): 186-9. doi:10.1038/nature01149](https://pubmed.ncbi.nlm.nih.gov/12432395/)

__What about cell-free networks__?
No cells, but most of the biochemistry. Is there an objective? This is an open question. However, initial work in this area (for networks where gene expression is occurring) suggests the maximization of translation, at least for single gene systems:

* [Vilkhovoy, Michael et al. “Sequence-Specific Modeling of E. coli Cell-Free Protein Synthesis.” ACS synthetic biology vol. 7,8 (2018): 1844-1857. doi:10.1021/acssynbio.7b00465](https://pubmed.ncbi.nlm.nih.gov/29944340/)

"""

# ╔═╡ d2c038fe-486d-11ec-1841-8972a51688e6
html"""
<style>
main {
    max-width: 1200px;
    width: 85%;
    margin: auto;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
MAT = "~0.10.1"
PlutoUI = "~0.7.19"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "0bc60e3006ad95b4bb7497698dd7c6d649b9bc06"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Blosc]]
deps = ["Blosc_jll"]
git-tree-sha1 = "217da19d6f3a94753e580a8bc241c7cbefd9281f"
uuid = "a74b3585-a348-5f62-a45c-50e91977d574"
version = "0.7.1"

[[Blosc_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Lz4_jll", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "e747dac84f39c62aff6956651ec359686490134e"
uuid = "0b7ba130-8d10-5ba8-a3d6-c5182647fed9"
version = "1.21.0+0"

[[BufferedStreams]]
deps = ["Compat", "Test"]
git-tree-sha1 = "5d55b9486590fdda5905c275bb21ce1f0754020f"
uuid = "e1450e63-4bb3-523b-b2a4-4ffa8c0fd77d"
version = "1.0.0"

[[CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dce3e3fea680869eaa0b774b2e8343e9ff442313"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.40.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[HDF5]]
deps = ["Blosc", "Compat", "HDF5_jll", "Libdl", "Mmap", "Random", "Requires"]
git-tree-sha1 = "698c099c6613d7b7f151832868728f426abe698b"
uuid = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
version = "0.15.7"

[[HDF5_jll]]
deps = ["Artifacts", "JLLWrappers", "LibCURL_jll", "Libdl", "OpenSSL_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "fd83fa0bde42e01952757f01149dd968c06c4dba"
uuid = "0234f1f7-429e-5d53-9886-15a909be8d59"
version = "1.12.0+1"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Lz4_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "5d494bc6e85c4c9b626ee0cab05daa4085486ab1"
uuid = "5ced341a-0733-55b8-9ab6-a4889d929147"
version = "1.9.3+0"

[[MAT]]
deps = ["BufferedStreams", "CodecZlib", "HDF5", "SparseArrays"]
git-tree-sha1 = "5c62992f3d46b8dce69bdd234279bb5a369db7d5"
uuid = "23992714-dd62-5051-b70f-ba57cb901cac"
version = "0.10.1"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.10+0"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "ae4bbcadb2906ccc085cf52ac286dc1377dceccc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "e071adf21e165ea0d904b595544a8e514c8bb42c"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.19"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─3a51e7a0-48e6-4b02-8649-9708ccb6ffa5
# ╟─7299f66a-2ca4-4a53-b5f2-bf249fe436f3
# ╟─328bbf10-c6ec-4d3e-8ffb-0a09c327920e
# ╟─3642b338-d819-4c0b-a93f-59006ce04174
# ╟─d92b38f1-1604-4b4d-8288-4c94f58d9869
# ╠═4b426b24-f416-4079-9951-275164d792f3
# ╠═14ccf90a-7e77-4112-914f-7b81917389ed
# ╠═1edb1823-001d-4e16-bc66-01f8136fbf1c
# ╠═ec8ed359-94e2-4a3a-a08c-386bd4166d00
# ╠═ca63dac5-7fa4-498e-af3e-52c917be3dc2
# ╠═a5b8c0b0-fd80-47b9-ad52-6267e50b0539
# ╟─5d81d380-df70-4480-aaa6-6bde2ce60826
# ╟─181c2a31-8196-4f00-87aa-b828144d1963
# ╟─4b7dcaec-f7a0-44ad-898e-fbca9f5271e5
# ╟─944bd174-5e89-42c0-aef0-a3c4450e8194
# ╟─d2c038fe-486d-11ec-1841-8972a51688e6
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
module PandorasBoxGate

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/gate" Gate(x::Integer, y::Integer, inverted::Bool=false, flag::String="", texture::String="objects/pandorasBox/gate/")

const placements = Ahorn.PlacementDict(
    "Gate (Pandora's Box)" => Ahorn.EntityPlacement(
        Gate
    )
)

function getSprite(gate::Gate)
    texture = get(gate, "texture", "objects/pandorasBox/gate/")

    return texture * "gate0"
end

function Ahorn.selection(gate::Gate)
    x, y = Ahorn.position(gate)
    sprite = getSprite(gate)

    customSprite = Ahorn.getSprite(sprite, "Gameplay")

    if customSprite.width == 0 || customSprite.height == 0
        return Ahorn.Rectangle(x - 4, y - 4, 8, 8)

    else
        return Ahorn.getSpriteRectangle(sprite, x - 4, y, jx=0.0, jy=0.0)
    end
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, gate::Gate, room::Maple.Room)
    sprite = getSprite(gate)

    Ahorn.drawSprite(ctx, sprite, -4, 0, jx=0.0, jy=0.0)
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
import Base: size

const NDIM = 2

import CartesianGrids: cellsize, origin
import RigidBodyTools: assign_velocity!
import ImmersedLayers: normals, areas

"""
$(TYPEDEF)

A system type that utilizes a grid of `NX` x `NY` dual cells and `N` Lagrange forcing
points to solve the discrete Navier-Stokes equations in vorticity form. The
parameter `static_points` specifies whether the forcing points remain static in the
grid. It should be set to `false` if a supplied motion requires that the points move.

# Constructors:

`NavierStokes(Re,Δx,xlimits,ylimits,Δt
              [,freestream = (0.0, 0.0)]
              [,pulses=nothing])` specifies the Reynolds number `Re`, the grid
              spacing `Δx`, the dimensions of the domain in the tuples `xlimits`
              and `ylimits` (excluding the ghost cells), and the time step size `Δt`.
              The other arguments are optional. The `freestream` argument can be
              passed as either a tuple (a static freestream) or a `RigidBodyMotion`
              for a time-varying freestream. The `pulses` argument can be
              used to pass in one or more spatiotemporal pulses.

`NavierStokes(Re,Δx,xlimits,ylimits,Δt,bodies::Body/BodyList
              [,flow_side=ExternalInternalFlow]
              [,ddftype=CartesianGrids.Yang3])` passes the body
              information. This constructor
              sets the motions of the body/ies to be stationary.
              The same optional arguments used for the basic constructor
              also apply for this one. In addition, the `flow_side` can be set to
              `ExternalFlow` (default), `InternalFlow`, or `ExternalInternalFlow`.
              However, it is forced to `ExternalInternalFlow` for open Bodies
              (like `Plate` type).

`NavierStokes(Re,Δx,xlimits,ylimits,Δt,bodies::Body/BodyList,
              motions::RigidBodyMotion/RigidMotionList
              [,static_points=false])`
              passes the body and associated motion information.
              The list of motions must be the same length as the list of bodies.
              The same optional arguments used for the other constructors
              also apply for this one. In addition, `static_points` can
              be set to `true` if the supplied motion should not cause the
              points to move.

"""
mutable struct NavierStokes{NX, NY, N, MT<:PointMotionType, FS<:FreestreamType, SD<:FlowSide, DDF<:CartesianGrids.DDFType, FT, SP, FST} #, RKT}
    # Physical Parameters
    "Reynolds number"
    Re::Float64
    "Free stream velocities"
    U∞::FST
    "Bodies"
    bodies::Union{BodyList,Nothing}
    "Body motions"
    motions::Union{RigidMotionList,Nothing}
    "Pulses"
    pulses::Union{Vector{ModulatedField},Nothing}

    # Discretization
    "Grid metadata"
    grid::CartesianGrids.PhysicalGrid{2}
    "Time step"
    Δt::Float64

    # Operators
    "Laplacian operator"
    L::CartesianGrids.Laplacian

    # Layers
    dlf::Union{DoubleLayer,Nothing} # used for viscous surface terms
    slc::Union{SingleLayer,Nothing} # used for scalar potential field in velocity
    sln::Union{SingleLayer,Nothing} # might not be used

    # Body coordinate data, if present
    # if a static problem, these coordinates are in inertial coordinates
    # if a non-static problem, in their own coordinate systems
    points::VectorData{N,Float64}

    # Pre-stored regularization and interpolation matrices (if present)
    Rf::Union{RegularizationMatrix,Nothing} # faces (edges)
    Ef::Union{InterpolationMatrix,Nothing}
    Cf::Union{AbstractMatrix,Nothing}
    Rc::Union{RegularizationMatrix,Nothing} # cell centers
    Ec::Union{InterpolationMatrix,Nothing}
    Rn::Union{RegularizationMatrix,Nothing} # cell nodes
    En::Union{InterpolationMatrix,Nothing}

    # Operators
    f :: FT

    # state vector
    state_prototype :: SP

    # Scratch space

    ## Pre-allocated space for intermediate values
    Vb::VectorData{N,Float64}
    Sb::ScalarData{N,Float64}
    Δus::VectorData{N,Float64}
    τ::VectorData{N,Float64}
    Vf::Edges{Primal, NX, NY, Float64}
    Vv::Edges{Primal, NX, NY, Float64}
    Vn::Edges{Primal, NX, NY, Float64}
    Sc::Nodes{Primal, NX, NY,Float64}
    Sn::Nodes{Dual, NX, NY,Float64}
    Wn::Nodes{Dual, NX, NY,Float64}
    Vtf::EdgeGradient{Primal,Dual,NX,NY,Float64}
    DVf::EdgeGradient{Primal,Dual,NX,NY,Float64}
    VDVf::EdgeGradient{Primal,Dual,NX,NY,Float64}

end

function NavierStokes(Re::Real, Δx::Real, xlimits::Tuple{Real,Real},ylimits::Tuple{Real,Real}, Δt::Real;
                       freestream::FST = (0.0, 0.0),
                       bodies::Union{BodyList,Nothing} = nothing,
                       motions::Union{RigidMotionList,Nothing} = nothing,
                       pulses::PT = nothing,
                       static_points = true,
                       flow_side::Type{SD} = ExternalFlow,
                       ddftype=CartesianGrids.Yang3) where {FST,PT,SD<:FlowSide}

    g = PhysicalGrid(xlimits,ylimits,Δx)
    NX, NY = size(g)

    α = Δt/(Re*Δx^2)

    # Set up buffers
    Vf = Edges{Primal,NX,NY,Float64}()
    Vv = Edges{Primal,NX,NY,Float64}()
    Vn = Edges{Primal,NX,NY,Float64}()
    Sc = Nodes{Primal,NX,NY,Float64}()
    Sn = Nodes{Dual,NX,NY,Float64}()
    Wn = Nodes{Dual,NX,NY,Float64}()
    Vtf = EdgeGradient{Primal,Dual,NX,NY,Float64}()
    DVf = EdgeGradient{Primal,Dual,NX,NY,Float64}()
    VDVf = EdgeGradient{Primal,Dual,NX,NY,Float64}()

    L = plan_laplacian(Sn,with_inverse=true)

    pulsefields = _process_pulses(pulses,Sn,g)

    # for now, if there are any bodies that are Open,
    # then force flow_side to ExternalInternalFlow.
    # but should be more flexible here
    flow_side_internal = _any_open_bodies(bodies) ? ExternalInternalFlow : flow_side

    N = numpts(bodies)

    Vb = VectorData(N)
    Sb = ScalarData(N)
    Δus = VectorData(N)
    τ = VectorData(N)

    points, dlf, slc, sln, Rf, Ef, Cf, Rc, Ec, Rn, En =
              _immersion_operators(bodies,g,flow_side_internal,ddftype,Vf,Sc,Vb)

    viscous_L = plan_laplacian(Sn,factor=1/(Re*Δx^2))

    if isnothing(bodies)
      state_prototype = solvector(state=Sn)
      f = ConstrainedODEFunction(ns_rhs!,viscous_L,_func_cache=state_prototype)
    else
      if static_points
        state_prototype = solvector(state=Sn,constraint=τ)
        f = ConstrainedODEFunction(ns_rhs!,bc_constraint_rhs!,
                                      ns_op_constraint_force!,bc_constraint_op!,
                                      viscous_L,_func_cache=state_prototype)
      else
        state_prototype = solvector(state=Sn,constraint=τ,aux_state=zero_body_state(bodies))
        rhs! = ConstrainedSystems.r1vector(state_r1 = ns_rhs!,aux_r1 = rigid_body_rhs!)
        f = ConstrainedODEFunction(rhs!,bc_constraint_rhs!,
                                      ns_op_constraint_force!,bc_constraint_op!,
                                      viscous_L,_func_cache=state_prototype,
                                      param_update_func=update_immersion_operators!)
      end

    end

    NavierStokes{NX, NY, N, _motiontype(static_points), _fstype(FST), flow_side_internal, ddftype, typeof(f), typeof(state_prototype),FST}( #,typeof(rk)}(
                          Re, freestream, bodies, motions, pulsefields,
                          g, Δt, # rk,
                          L,
                          dlf,slc,sln,
                          points, Rf, Ef, Cf, Rc, Ec, Rn, En,
                          f,state_prototype,
                          Vb, Sb, Δus, τ,
                          Vf, Vv, Vn, Sc, Sn, Wn, Vtf, DVf, VDVf)
end

NavierStokes(Re,Δx,xlim,ylim,Δt,bodies::BodyList;
        motions=RigidMotionList(map(x -> RigidBodyMotion(0.0,0.0),bodies)),kwargs...) =
        NavierStokes(Re,Δx,xlim,ylim,Δt;bodies=bodies,motions=motions,kwargs...)

NavierStokes(Re,Δx,xlim,ylim,Δt,body::Body;kwargs...) =
        NavierStokes(Re,Δx,xlim,ylim,Δt,BodyList([body]);kwargs...)

function NavierStokes(Re,Δx,xlim,ylim,Δt,bodies::BodyList,motions::RigidMotionList;static_points=false,kwargs...)
    length(bodies) == length(motions) || error("Inconsistent lengths of bodies and motions lists")
    NavierStokes(Re,Δx,xlim,ylim,Δt,bodies;motions=motions,static_points=static_points,kwargs...)
end

NavierStokes(Re,Δx,xlim,ylim,Δt,body::Body,motion::RigidBodyMotion;static_points=false,kwargs...) =
        NavierStokes(Re,Δx,xlim,ylim,Δt,BodyList([body]),RigidMotionList([motion]);static_points=static_points,kwargs...)

function Base.show(io::IO, sys::NavierStokes{NX,NY,N,MT,FS,SD}) where {NX,NY,N,MT,FS,SD}
    mtype = (MT == StaticPoints) ? "static" : "moving"
    fsmsg = (FS == StaticFreestream) ? "Static freestream = $(sys.U∞)" : "Variable freestream"
    sdmsg = (N == 0) ? "Unbounded" : ((SD == ExternalFlow) ? "External flow" : ((SD == InternalFlow) ? "Internal flow" : "External/internal"))
    println(io, "$sdmsg Navier-Stokes system on a grid of size $NX x $NY and $N $mtype immersed points")
    println(io, "   $fsmsg")
    if N > 0
      bdmsg = (length(sys.bodies) == 1) ? "1 body" : "$(length(sys.bodies)) bodies"
      println(io, "   $bdmsg")
    end
end

# Routines to set up the immersion operators

function _immersion_operators(bodies::BodyList,g::PhysicalGrid,flow_side::Type{SD},ddftype,Vf::GridData{NX,NY},Sc::GridData{NX,NY},Vb::PointData{N}) where {NX,NY,N,SD<:ViscousFlow.FlowSide}

  points = VectorData(collect(bodies))
  numpts(bodies) == N || error("Inconsistent size of bodies")

  body_areas = areas(bodies)
  body_normals = normals(bodies)

  if !(flow_side==ExternalInternalFlow)
    dlf = DoubleLayer(bodies,g,Vf)
    slc = SingleLayer(bodies,g,Sc)
    #sln = SingleLayer(bodies,g,Sn)
    sln = nothing
  else
    dlf = nothing
    slc = nothing
    sln = nothing
  end

  #regop = Regularize(points,cellsize(g);I0=CartesianGrids.origin(g),weights=body_areas.data,ddftype=ddftype)
  regop = _regularization(points,g,bodies,ddftype)

  Rf = RegularizationMatrix(regop,Vb,Vf) # Used by B₁ᵀ
  Ef = InterpolationMatrix(regop,Vf,Vb) # Used by constraint_rhs! and B₂

  #Rc = RegularizationMatrix(regop,Sb,Sc)
  #Ec = InterpolationMatrix(regop,Sc,Sb)
  #Rn = RegularizationMatrix(regop,Sb,Sn)
  #En = InterpolationMatrix(regop,Sn,Sb)
  Rc = nothing
  Ec = nothing
  Rn = nothing
  En = nothing

  regopfilt = Regularize(points,cellsize(g);I0=CartesianGrids.origin(g),filter=true,weights=cellsize(g)^2,ddftype=ddftype)
  Ẽf = InterpolationMatrix(regopfilt,Vf,Vb)
  Cf = sparse(Ẽf*Rf)

  return points, dlf, slc, sln, Rf, Ef, Cf, Rc, Ec, Rn, En
end

_immersion_operators(::Nothing,a...) =
    VectorData(0), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing

# For updating the system with body data

function update_immersion_operators!(sys::NavierStokes{NX,NY,N,MT,FS,SD,DDF},bodies::BodyList) where {NX,NY,N,MT,FS,SD,DDF<:CartesianGrids.DDFType}
    sys.bodies = deepcopy(bodies)
    sys.points, sys.dlf, sys.slc, sys.sln, sys.Rf, sys.Ef, sys.Cf, sys.Rc, sys.Ec, sys.Rn, sys.En =
      _immersion_operators(sys.bodies,sys.grid,SD,DDF,sys.Vf,sys.Sc,sys.Vb)
    return sys
end

update_immersion_operators!(sys::NavierStokes,body::Body) = update_immersion_operators!(sys,BodyList([body]))

function update_immersion_operators!(sys::NavierStokes,x::AbstractVector)
    tl! = RigidTransformList(x)
    tl!(sys.bodies)
    update_immersion_operators!(sys,sys.bodies)
end

# The form passed to ConstrainedODEFunction
update_immersion_operators!(sys::NavierStokes,u,sys_old::NavierStokes,t) =
    update_immersion_operators!(sys,aux_state(u))

"""
    setstepsizes(Re[,gridRe=2][,cfl=0.5][,fourier=0.5]) -> Float64, Float64

Set the grid cell spacing and time step size based on the Reynolds number `Re`,
the grid Reynolds number `gridRe`, cfl number `cfl`, and grid Fourier number `fourier`.
The last three parameters all have default values.

# Example

Here is an example of setting parameters based on Reynolds number 100 (with
  default choices for grid Reynolds number, CFL number, and Fourier number):
```jldoctest
julia> Δx, Δt = setstepsizes(100)
(0.02, 0.01)
```
"""
function setstepsizes(Re::Real; gridRe = 2.0, cfl = 0.5, fourier = 0.5)
    Δx = gridRe/Re
    Δt = min(fourier*Δx,cfl*Δx^2*Re)
    return Δx, Δt
end

# some convenience functions
"""
    size(sys::NavierStokes,d::Int) -> Int

Return the number of indices of the grid used by `sys` along dimension `d`.
"""
size(sys::NavierStokes{NX,NY},d::Int) where {NX,NY} = d == 1 ? NX : NY

"""
    size(sys::NavierStokes) -> Tuple{Int,Int}

Return a tuple of the number of indices of the grid used by `sys`
"""
size(sys::NavierStokes{NX,NY}) where {NX,NY} = (size(sys,1),size(sys,2))

"""
    cellsize(sys::NavierStokes) -> Float64

Return the grid cell size of system `sys`
"""
cellsize(sys::NavierStokes) = cellsize(sys.grid)

"""
    timestep(sys::NavierStokes) -> Float64

Return the time step size of system `sys`
"""
timestep(sys::NavierStokes) = sys.Δt

"""
    origin(sys::NavierStokes) -> Tuple{Int,Int}

Return a tuple of the indices of the primal node that corresponds to the
physical origin of the coordinate system used by `sys`. Note that these
indices need not lie inside the range of indices occupied by the grid.
For example, if the range of physical coordinates occupied by the grid
is (1.0,3.0) x (2.0,4.0), then the origin is not inside the grid.
"""
origin(sys::NavierStokes) = origin(sys.grid)

"""
    timerange(tf,sys::NavierStokes)

Create a range of times, starting at the t = Δt (the time step of `sys`),
and ending at t = `tf`.
"""
timerange(tf,sys) = timestep(sys):timestep(sys):tf

# Wrap the output of the motion evaluation in VectorData
@inline assign_velocity!(u::VectorData,a...) = (assign_velocity!(u.u,u.v,a...); u)

@inline normals(sys::NavierStokes) = (!isnothing(sys.bodies)) ? normals(sys.bodies) : nothing

"""
    freestream(t,sys::NavierStokes) -> Tuple

Return the value of the freestream in `sys` at time `t` as a tuple.
"""
freestream(t::Real,sys::NavierStokes{NX,NY,N,MT,StaticFreestream}) where {NX,NY,N,MT} = sys.U∞

function freestream(t::Real,sys::NavierStokes{NX,NY,N,MT,VariableFreestream}) where {NX,NY,N,MT}
    _,ċ,_,_,_,_ = sys.U∞(t)
    return reim(ċ)
end

"""
    newstate(sys::NavierStokes)

Return a new (zero) instance of the state vector for `sys`.
"""
newstate(sys::NavierStokes) = zero(sys.state_prototype)

"""
    newstate(s::AbstractSpatialField,sys::NavierStokes)

Return an instance of the state vector for `sys`, assigned the
data in the spatial field `s`.
"""
function newstate(s::AbstractSpatialField,sys::NavierStokes)
  u = newstate(sys)
  gf = GeneratedField(state(u),s,sys.grid)
  state(u) .= cellsize(sys)*gf()
  return u
end

# Other functions
_hasfilter(sys::NavierStokes) = !isnothing(sys.Cf)
_motiontype(isstatic::Bool) = isstatic ? StaticPoints : MovingPoints
  _fstype(F) = F <: Union{RigidBodyMotion,Kinematics} ? VariableFreestream : StaticFreestream

_body_closure_type(b::T) where {T<:Body{N,C}} where {N,C} = C

_any_open_bodies(nothing) = false
_any_open_bodies(bodies::BodyList) =  any(b -> _body_closure_type(b) == RigidBodyTools.OpenBody,bodies)

_regularization(sys::NavierStokes{NX, NY, N, MT, FS, SD, DDF}) where {NX,NY,N,MT,FS,SD,DDF} =
        _regularization(sys.points,sys.grid,sys.bodies,DDF)

_regularization(points,g,bodies,ddftype) = Regularize(points,cellsize(g),
                                I0=CartesianGrids.origin(g),weights=areas(bodies).data,ddftype=ddftype)

include("navierstokes/surfacevelocities.jl")
include("navierstokes/fields.jl")
include("navierstokes/pointforce.jl")
include("navierstokes/basicoperators.jl")
include("navierstokes/rigidbodyoperators.jl")
include("navierstokes/movingbodyoperators.jl")
include("navierstokes/timemarching.jl")


text:
using Pkg, BinDeps, PyCall
using PyCall: python, conda, Conda
using BinDeps: generate_steps, getallproviders, lower, PackageManager

!Sys.islinux() && exit()

if isnothing(Sys.which("sudo")) # in docker
    try run(`apt update`) catch end
    try run(`yum update`) catch end
end

@BinDeps.setup

gcc = library_dependency("gcc")
dot = library_dependency("dot")
python3 = library_dependency("python3")

provides(AptGet, Dict("g++" => gcc, "python3-dev" => python3, "graphviz" => dot))
provides(Yum, Dict("gcc-c++" => gcc, "python3-devel" => python3, "graphviz" => dot))

for dep in bindeps_context.deps
    dp, opts = getallproviders(dep, PackageManager)[1]
    cmd = lower(generate_steps(dep, dp, opts)).steps[1]
    i = findfirst(x -> x == "install", cmd.exec)
    insert!(cmd.exec, i + 1, "-y")
    println(cmd)
    try run(cmd) catch end
end

for pkg in ["numpy<=1.20", "pandas", "scikit-learn<=0.22.2", "matplotlib",
            "lightgbm", "ipython", "shap", "keras", "tzlocal", "PyPDF2",
            "unidecode", "pdpbox", "git+https://github.com/oracle/Skater.git"]
    try
        run(`$python -m pip install "$pkg"`)
    catch e
        println(e)
    end
end

buildsh = joinpath(@__DIR__, "build.sh")
ENV["JULIA_DEPOT_PATH"] = DEPOT_PATH[1]
run(`bash $buildsh`)


text:
import JuLIP.Chemistry: rnn

export rattle!, r_sum, r_dot,
      swapxy!, swapxz!, swapyz!,
      dist, displacement, rmin

############################################################
###   Some useful utility functions
############################################################
###   Robust Summation and Dot Products
############################################################

"Robust summation. Uses `sum_kbn`."
r_sum(a) = sum_kbn(a)

## NOTE: if I see this correctly, then r_dot allocates a temporary
##       vector, which is likely quite a performance overhead.
##       probably, we want to re-implement this.
## TODO: new version without the intermediate allocation
"Robust inner product. Defined as `r_dot(a, b) = r_sum(a .* b)`"
r_dot(a, b) = r_sum(a[:] .* b[:])

"""
`rattle!(at::AbstractAtoms, r::Float64; rnn = 1.0, respect_constraint = true)
  -> at`

randomly perturbs the atom positions

* `r`: magnitude of perturbation
* `rnn` : nearest-neighbour distance
* `respect_constraint`: set false to also perturb the constrained atom positions
"""
function rattle!(at::AbstractAtoms, r::Float64; rnn = 1.0, respect_constraint = true)
   # if there is no constraint, then revert to respect_constraint = false
   if isa(constraint(at), NullConstraint)
      respect_constraint = false
   end
   if respect_constraint
      x = dofs(at)
      x += r * rnn * 2.0/sqrt(3) * (rand(length(x)) - 0.5)
      set_dofs!(at, x)
   else
      X = positions(at) |> mat
      X += r * rnn * 2.0/sqrt(3) * (rand(size(X)) - 0.5)
      set_positions!(at, X)
   end
   return at
end

"""
use this instead of `warn`, then warnings can be turned off by setting
`Main.JULIPWARN=false`
"""
function julipwarn(s)
   if isdefined(Main, :JULIPWARN)
      if Main.JULIPWARN == false
         return false
      end
   end
   warn(s)
end

"""
swap x and y position coordinates
"""
function swapxy!(at::AbstractAtoms)
   X = positions(at) |> mat
   X[[2,1],:] = X[[1,2],:]
   set_positions!(at, X)
   return at
end

"""
swap x and z position coordinates
"""
function swapxz!(at::AbstractAtoms)
   X = positions(at) |> mat
   X[[3,1],:] = X[[1,3],:]
   set_positions!(at, X)
   return at
end

"""
swap y and z position coordinates
"""
function swapyz!(at::AbstractAtoms)
   X = positions(at) |> mat
   X[[3,2],:] = X[[2,3],:]
   set_positions!(at, X)
   return at
end

"""
`dist(at, X1, X2, p = Inf)`
`dist(at1, at2, p = Inf)`

Returns the maximum distance (p = Inf) or alternatively a p-norm of
distances between the two configurations `X1, X2` or `at1, at2`.
This implementation accounts for periodic boundary conditions (in those
coordinate directions where they are set to `true`)
"""
function dist{T}(at::AbstractAtoms,
                 X1::Vector{JVec{T}}, X2::Vector{JVec{T}},
                 p = Inf)
   @assert length(X1) == length(X2)
   F = cell(at)'
   Finv = inv(F)
   bcrem = [ p ? 1.0 : Inf for p in pbc(at) ]
   d = [ norm(_project_pbc_(F, Finv, bcrem, x1 - x2))
         for (x1, x2) in zip(X1, X2) ]
   return norm(d, p)
end

dist(at::AbstractAtoms, X::Vector) = dist(at, positions(at), X)

function dist(at1::AbstractAtoms, at2::AbstractAtoms, p = Inf)
   @assert vecnorm(cell(at1) - cell(at2), Inf) < 1e-14
   return dist(at1, positions(at1), positions(at2), p)
end

function _project_pbc_(F, Finv, bcrem, x)
   λ = Finv * x     # convex coordinates
   # convex coords projected to the unit
   λp = JVecF(rem(λ[1], bcrem[1]), rem(λ[2], bcrem[2]), rem(λ[3], bcrem[3]))
   return F * λp     # convert back to real coordinates
end

function _project_coord_min_(λ, p)
   if !p
      return λ
   end
   λ = mod(λ, 1.0)   # project to cell
   if λ > 0.5        # periodic image with minimal length
      λ = λ - 1.0
   end
   return λ
end

function _project_pbc_min_(F, Finv, p, x)
   λ = Finv * x     # convex coordinates
   # convex coords projected to the unit
   λp = _project_coord_min_.(λ, JVec{Bool}(p...))
   return F * λp     # convert back to real coordinates
end

function displacement{T}(at::AbstractAtoms, X1::Vector{JVec{T}}, X2::Vector{JVec{T}})
   @assert length(X1) == length(X2)
   F = defm(at)
   Finv = inv(F)
   p = pbc(at)
   U = [ _project_pbc_min_(F, Finv, p, x2-x1)
         for (x1, x2) in zip(X1, X2) ]
   return U
end

project_min(at, u) = _project_pbc_min_(defm(at), inv(defm(at)), pbc(at), u)

function rnn(at::AbstractAtoms)
   syms = unique(chemical_symbols(at))
   return minimum(rnn.(syms))
end

function wrap_pbc!(at)
   X = positions(at)
   F = defm(at)
   X = [_project_pbc_min_(F, inv(F), pbc(at), x) for x in X]
   set_positions!(at, X)
end

function rmin(at::AbstractAtoms)
   at2 = at * 2
   X = positions(at2)
   r = norm(X[1]-X[2])
   for n = 1:length(X)-1, m = n+1:length(X)
      r = min(r, norm(X[n]-X[m]))
   end
   return r
end


text:
# rmsd100
# root mean-squared difference between two arrays in %

export rmsd100

using LinearAlgebra: norm

"""
    rmsd = rmsd100(x, y ; mask)

Compute 100 * RMSD (root mean squared difference) between `x` and `y` within domain mask.

in
- `x` : array
- `y` : another array of same size

option:
- `mask::Array{Bool}` : domain over which to compute the RMSE; default `trues(size(x))`

out
- rmsd : rmsd of `x` vs `y` within `mask` in %
"""
function rmsd100(x::AbstractArray{<:Number}, y::AbstractArray{<:Number}
        ; mask::AbstractArray{Bool} = trues(size(x)))
    return 100 * norm(x[mask] - y[mask]) / sqrt(sum(mask))
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %% [markdown]
# https://github.com/haampie/ArnoldiMethod.jl

# %%
]st

# %%
using ArnoldiMethod, LinearAlgebra, SparseArrays

# %%
A = spdiagm(
    -1 => fill(-1.0, 99),
    0 => fill(2.0, 100),
    1 => fill(-1.0, 99)
)

# %%
decomp, history = partialschur(A, nev=10, tol=1e-6, which=SR());

# %%
history

# %%
decomp

# %%
"""
    Rev0

* This revision is almost equivalent to the original.
* Stop using constants.
* Always pass global variables to functions as arguments.
* Swap the order of the for loop.
* Add @inbounds macro.
* Revise sum([f(x) for x in X]) to sum(f(x) for x in X).
"""
module Rev0

using LinearAlgebra, LinearMaps
import Arpack

function prepareDiag(L)
    diag = zeros(2^L)
    for state = 1:2^L
        for i = 1:L
            j = i==L ? 1 : i+1
            @inbounds diag[state] -= (((state >> (i-1))&1) == ((state >> (j-1))&1)) ? 1 : -1
        end
    end
    diag
end

function Hfunc!(C, B, diag, L)
    for state = 1:2^L
        @inbounds C[state] = diag[state] * B[state]
    end
    for i = 1:L
        for state = 1:2^L
            newstate = (state&(~(2^L))) ⊻ (1<<(i-1))
            if newstate == 0
                newstate = 2^L
            end
            @inbounds C[newstate] -= B[state]
        end
    end
end
prepareHfunc!(diag, L) = (C, B) -> Hfunc!(C, B, diag, L)

L = 20
diag_ = prepareDiag(L)
H = LinearMap(prepareHfunc!(diag_, L), 2^L, ismutating=true, issymmetric=true, isposdef=false)
sol_exact = -2sum(abs(sin((n-1/2) * pi/L)) for n in 1:L)

end;

# %%
"""
    Rev1

* Use Threads.@threads macro.
"""
module Rev1

using LinearAlgebra, LinearMaps
import Arpack

function prepareDiag(L)
    diag = zeros(2^L)
    for state = 1:2^L
        for i = 1:L
            j = i==L ? 1 : i+1
            @inbounds diag[state] -= (((state >> (i-1))&1) == ((state >> (j-1))&1)) ? 1 : -1
        end
    end
    diag
end

function Hfunc!(C, B, diag, L)
    Threads.@threads for state = 1:2^L
        @inbounds C[state] = diag[state] * B[state]
    end
    for i = 1:L
        Threads.@threads for state = 1:2^L
            newstate = (state&(~(2^L))) ⊻ (1<<(i-1))
            if newstate == 0
                newstate = 2^L
            end
            @inbounds C[newstate] -= B[state]
        end
    end
end
prepareHfunc!(diag, L) = (C, B) -> Hfunc!(C, B, diag, L)

L = 20
diag_ = prepareDiag(L)
H = LinearMap(prepareHfunc!(diag_, L), 2^L, ismutating=true, issymmetric=true, isposdef=false)
sol_exact = -2sum(abs(sin((n-1/2) * pi/L)) for n in 1:L)

end;

# %%
"""
    Rev2

* Use LoopVectorization.@tturbo macro.
"""
module Rev2

using LinearAlgebra, LinearMaps
import Arpack
using LoopVectorization

function prepareDiag(L)
    diag = zeros(2^L)
    for state = 1:2^L
        for i = 1:L
            j = i==L ? 1 : i+1
            @inbounds diag[state] -= (((state >> (i-1))&1) == ((state >> (j-1))&1)) ? 1 : -1
        end
    end
    diag
end

function Hfunc!(C, B, diag, L)
    N = length(diag)
    @tturbo for state = 1:N
        C[state] = diag[state] * B[state]
    end
    for i = 1:L
        @tturbo for state = 1:N
            newstate = (state&(~(2^L))) ⊻ (1<<(i-1))
            c = newstate == 0
            newstate = !c*newstate + c*N # remove if statement
            C[newstate] -= B[state]
        end
    end
end
prepareHfunc!(diag, L) = (C, B) -> Hfunc!(C, B, diag, L)

L = 20
diag_ = prepareDiag(L)
H = LinearMap(prepareHfunc!(diag_, L), 2^L, ismutating=true, issymmetric=true, isposdef=false)
sol_exact = -2sum(abs(sin((n-1/2) * pi/L)) for n in 1:L)

end;

# %%
Rev0.sol_exact

# %%
@time decomp0, history0 = partialschur(Rev0.H, nev=1, which=SR())
@time decomp0, history0 = partialschur(Rev0.H, nev=1, which=SR())
decomp0.eigenvalues[1]

# %%
@time decomp1, history1 = partialschur(Rev1.H, nev=1, which=SR())
@time decomp1, history1 = partialschur(Rev1.H, nev=1, which=SR())
decomp1.eigenvalues[1]

# %%
@time decomp2, history2 = partialschur(Rev2.H, nev=1, which=SR())
@time decomp2, history2 = partialschur(Rev2.H, nev=1, which=SR())
decomp2.eigenvalues[1]

# %%

max_stars_repo_path: 0015/Arpack/ArnoldiMethod.jl example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0015/Arpack/ArnoldiMethod.jl example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0015/Arpack/ArnoldiMethod.jl example.jl

max_forks_repo_name: genkuroki/public


text:
using LoopVectorization
using Base: OneTo
using LinearAlgebra: stride1
using Test

add_1_dim(x::AbstractArray) = reshape(x, size(x)..., 1)
check_finite(x::AbstractArray) = all(isfinite.(x)) || throw(error("x not finite!"))

"Given k-dimensional array `x` where `n=size(x,k)`, compute multinomial logistic Pr(i ∈ 1:n | x[:, ..., :, k] )"
function softmax3!(q::AA, lse::A, tmpmax::A, x::AA, maxk=size(q, ndims(q)) ) where {T<:Real, A<:Array{T}, AA<:AbstractArray{T}}
    ndims(q) == 1+ndims(lse) || throw(DimensionMismatch())
    xsizes = size(x)
    xsizes == size(q) || throw(DimensionMismatch("size(x) = $(size(x)) but size(q) = $(size(q))"))
    nk = last(xsizes)
    for i = OneTo(ndims(lse))
        size(q,i) == size(lse,i) == size(tmpmax,i) || throw(DimensionMismatch("size(x) = $(size(x)),  size(lse) = $(size(lse)), and size(tmpmax) = $(size(tmpmax))"))
    end
    0 < maxk <= nk || throw(DomainError(maxk))
    1 == stride1(q) == stride1(x) || throw(error("Arrays not strided"))

    isempty(x) && throw(error("x empty"))
    check_finite(x)

    maximum!(add_1_dim(tmpmax), x)
    fill!(lse, zero(T))

    xx = reshape(x, :, nk)
    qq = reshape(q, :, nk)

    for k in OneTo(nk)
        for i in eachindex(lse)
            tmp = exp(xx[i,k] - tmpmax[i])
            lse[i] += tmp
            k <= maxk && (qq[i,k] = tmp)
        end
    end

    qq[:,OneTo(maxk)] ./= vec(lse)
end

function softmax3_avx_broken!(q::AA, lse::A, tmpmax::A, x::AA, maxk=size(q, ndims(q)) ) where {T<:Real, A<:Array{T}, AA<:AbstractArray{T}}
        ndims(q) == 1+ndims(lse) || throw(DimensionMismatch())
        xsizes = size(x)
        xsizes == size(q) || throw(DimensionMismatch("size(x) = $(size(x)) but size(q) = $(size(q))"))
        nk = last(xsizes)
        for i = OneTo(ndims(lse))
            size(q,i) == size(lse,i) == size(tmpmax,i) || throw(DimensionMismatch("size(x) = $(size(x)),  size(lse) = $(size(lse)), and size(tmpmax) = $(size(tmpmax))"))
        end
        0 < maxk <= nk || throw(DomainError(maxk))
        1 == stride1(q) == stride1(x) || throw(error("Arrays not strided"))

        isempty(x) && throw(error("x empty"))
        check_finite(x)

        maximum!(add_1_dim(tmpmax), x)
        fill!(lse, zero(T))

        xx = reshape(x, :, nk)
        qq = reshape(q, :, nk)

        for k in OneTo(maxk)
            @avx for i in eachindex(lse)
                tmp = exp(xx[i,k] - tmpmax[i])
                lse[i] += tmp
                # FIXME - would prefer to replace 2nd loop w/ if stmt
                # k <= maxk && (qq[i,k] = tmp)
                qq[i,k] = tmp
            end
        end

        for k in maxk+1:nk
            @avx for i in eachindex(lse)
                tmp = exp(xx[i,k] - tmpmax[i])
                lse[i] += tmp
            end
        end

        qq[:,OneTo(maxk)] ./= vec(lse)
end

function softmax3_avx_fixed!(q::AA, lse::A, tmpmax::A, x::AA, maxk=size(q, ndims(q)) ) where {T<:Real, A<:Array{T}, AA<:AbstractArray{T}}
        ndims(q) == 1+ndims(lse) || throw(DimensionMismatch())
        xsizes = size(x)
        xsizes == size(q) || throw(DimensionMismatch("size(x) = $(size(x)) but size(q) = $(size(q))"))
        nk = last(xsizes)
        for i = OneTo(ndims(lse))
            size(q,i) == size(lse,i) == size(tmpmax,i) || throw(DimensionMismatch("size(x) = $(size(x)),  size(lse) = $(size(lse)), and size(tmpmax) = $(size(tmpmax))"))
        end
        0 < maxk <= nk || throw(DomainError(maxk))
        1 == stride1(q) == stride1(x) || throw(error("Arrays not strided"))

        isempty(x) && throw(error("x empty"))
        check_finite(x)

        maximum!(add_1_dim(tmpmax), x)
        fill!(lse, zero(T))

        xx = reshape(x, :, nk)
        qq = reshape(q, :, nk)
        tmpmaxvec = vec(tmpmax)  # Always required??
        lsevec = vec(lse)        # Always required??

        for k in OneTo(maxk)
            qk = view(qq, :, k) # required if using a View
            xk = view(xx, :, k) # required if using a View
            @avx for i in eachindex(lsevec)
                tmp = exp(xk[i] - tmpmaxvec[i])
                lsevec[i] += tmp
                qk[i] = tmp
            end
        end

        for k in maxk+1:nk
            qk = view(qq, :, k)
            xk = view(xx, :, k)
            @avx for i in eachindex(lsevec)
                tmp = exp(xk[i] - tmpmaxvec[i])
                lsevec[i] += tmp
            end
        end

        qq[:,OneTo(maxk)] ./= vec(lse)
end

ni, nj, nk = (100, 100, 10)
x = rand(ni, nj, nk)
q = similar(x)
tmpmax = zeros(ni,nj)
lse = similar(tmpmax)

for f! in (softmax3!, softmax3_avx_fixed!, softmax3_avx_broken!)
    @testset "$(f!) with arrays" begin
        fill!(q,0)
        fill!(lse,0)
        @show @code_warntype f!(q, lse, tmpmax, x, 1)
        f!(q, lse, tmpmax, x, 1)
        @test all(sum(q; dims=3) .<= 1)
        fill!(q,0)
        fill!(lse,0)
        f!(q, lse, tmpmax, x)
        @test sum(q; dims=3) ≈ ones(ni,nj)
    end

    @testset "$(f!) with views" begin
        nkm1 = nk-1
        @views qvw = q[:,:,1:nkm1]
        @views xvw = x[:,:,1:nkm1]
        fill!(q,0)
        fill!(lse,0)
        @show @code_warntype f!(qvw, lse, tmpmax, xvw, 1)
        f!(qvw, lse, tmpmax, xvw, 1)
        @test all(sum(qvw; dims=3) .<= 1)
        fill!(q,0)
        fill!(lse,0)
        f!(qvw, lse, tmpmax, xvw)
        @test sum(qvw; dims=3) ≈ ones(ni,nj)
    end
end


text:
"""
    rescale_estim_params!(coordinates_params_rescaling::Nothing, estim::EstimSettings, std_presample::FloatMatrix)

Default implementation for a generic `estim` does nothing.
"""
rescale_estim_params!(coordinates_params_rescaling::Nothing, estim::EstimSettings, std_presample::FloatMatrix) = nothing;

"""
    compute_loss_weights(custom_weights::FloatVector, n::Int64)

Use `custom_weights`.

    compute_loss_weights(custom_weights::Nothing, n::Int64)

Use equal weights.
"""
compute_loss_weights(custom_weights::FloatVector, n::Int64) = custom_weights;
compute_loss_weights(custom_weights::Nothing, n::Int64) = ones(n) ./ n;

"""
    compute_loss(weighted_se::AbstractArray{Float64})
    compute_loss(weighted_se::AbstractArray{Missing})
    compute_loss(weighted_se::AbstractArray{Union{Float64, Missing}})

Compute loss function. Return `[error inactive_sample]`.

`inactive_sample` is 1 when the data is entirely missing (0 otherwise).
"""
compute_loss(weighted_se::AbstractArray{Float64}) = [sum(mean(weighted_se, dims=2)), 0.0];
compute_loss(weighted_se::AbstractArray{Missing}) = [0.0, 1.0];

function compute_loss(weighted_se::AbstractArray{Union{Float64, Missing}})

    # Setup
    n, T = size(weighted_se);

    # Stop if the weighted squared errors are entirely missing
    if count(ismissing, weighted_se) == n*T
        return [0.0, 1.0];

    # Standard run
    else
        return [sum(sum_skipmissing(weighted_se))/T, 0.0];
    end
end

"""
    fc_err(validation_settings::ValidationSettings, p::Int64, λ::Number, α::Number, β::Number)

Compute the in-sample / out-of-sample error associated with the candidate hyperparameters

# Arguments
- `validation_settings`: ValidationSettings struct
- `p`: (candidate) number of lags in the vector autoregression
- `λ`: (candidate) overall shrinkage hyper-parameter for the elastic-net penalty
- `α`: (candidate) weight associated to the LASSO component of the elastic-net penalty
- `β`: (candidate) additional shrinkage for distant lags (p>1)
- `jth_jackknife_data`: j-th jackknife sample (default: nothing)

# References
Pellegrino (2022)
"""
function fc_err(validation_settings::ValidationSettings, p::Int64, λ::Number, α::Number, β::Number; jackknife_data::Union{JArray{Float64}, Nothing}=nothing, jackknife_sample_id::Union{Int64, Nothing}=nothing);

    # In-sample error
    if validation_settings.err_type == 1
        t0 = validation_settings.T;

    # Out-of-sample based error
    else
        t0 = validation_settings.t0;
    end

    # Jackknife out-of-sample
    if validation_settings.err_type > 2
        data = @view jackknife_data[:, :, jackknife_sample_id];

    # Standard in-sample or out-of-sample
    else
        data = validation_settings.Y;
    end

    # Standardise data
    data_presample = @view data[:, 1:t0];

    # Stop if the estimation sample is entirely missing
    if count(ismissing, data_presample) == validation_settings.n*t0
        return [0.0, 1.0];

    # Standard run
    else
        if validation_settings.is_stationary
            mean_presample = mean_skipmissing(data_presample);
            std_presample = std_skipmissing(data_presample);
            Y = (data .- mean_presample) ./ std_presample;
        else
            std_presample = std_skipmissing(diff(data_presample, dims=2));
            Y = data ./ std_presample;
        end

        estim = validation_settings.model_struct(Y[:, 1:t0], validation_settings.model_args, validation_settings.model_kwargs, p, λ, α, β);
        rescale_estim_params!(validation_settings.coordinates_params_rescaling, estim, std_presample);

        # Estimate model
        sspace = ecm(estim, output_sspace_data=Y);

        # Run Kalman filter
        status = kfilter_full_sample(sspace);

        # Forecast
        X_prior = mapreduce(Xt -> sspace.B*Xt, hcat, status.history_X_prior);

        # Apply original scaling
        if validation_settings.is_stationary
            forecast = (X_prior .* std_presample) .+ mean_presample;
        else
            forecast = X_prior .* std_presample;
        end

        # Compute weights
        w = compute_loss_weights(validation_settings.weights, validation_settings.n);

        # In-sample error
        if validation_settings.err_type == 1
            weighted_se = w .* (data .- forecast).^2;

        # Out-of-sample error
        else
            weighted_se = @views w .* (data[:, t0+1:end] .- forecast[:, t0+1:end]).^2;
        end

        # Return output
        return compute_loss(weighted_se);
    end
end

"""
    jackknife_err(validation_settings::ValidationSettings, jackknife_data::JArray{Float64, 3}, p::Int64, λ::Number, α::Number, β::Number)

Return the jackknife out-of-sample error.

# Arguments
- `validation_settings`: ValidationSettings struct
- `jackknife_data`: jackknife partitions
- `p`: (candidate) number of lags in the vector autoregression
- `λ`: (candidate) overall shrinkage hyper-parameter for the elastic-net penalty
- `α`: (candidate) weight associated to the LASSO component of the elastic-net penalty
- `β`: (candidate) additional shrinkage for distant lags (p>1)

# References
Pellegrino (2022)
"""
function jackknife_err(validation_settings::ValidationSettings, jackknife_data::JArray{Float64, 3}, p::Int64, λ::Number, α::Number, β::Number)

    # Error management
    if validation_settings.err_type <= 2
        error("Wrong err_type for jackknife_err!");
    end

    # Number of jackknife samples
    samples = size(jackknife_data, 3);

    # Compute jackknife loss
    verb_message(validation_settings.verb_estim, "jackknife_err > running $samples iterations on $(nworkers()) workers");

    output_fc_err = @sync @distributed (+) for j=1:samples
        fc_err(validation_settings, p, λ, α, β, jackknife_data=jackknife_data, jackknife_sample_id=j);
    end

    # Compute average jackknife loss
    loss, inactive_samples = output_fc_err;
    if samples == inactive_samples
        error("All samples are inactive! Check the initial settings or try a different randomisation method.");
    end
    loss *= 1/(samples-inactive_samples);

    # Return output
    return loss;
end

"""
    select_hyperparameters(validation_settings::ValidationSettings, γ_grid::HyperGrid)

Select the tuning hyper-parameters for the elastic-net vector autoregression.

# Arguments
- `validation_settings`: ValidationSettings struct
- `γ_grid`: HyperGrid struct

# References
Pellegrino (2022)
"""
function select_hyperparameters(validation_settings::ValidationSettings, γ_grid::HyperGrid, seed::Int64=1)

    # Check inputs
    check_bounds(validation_settings.n, 2); # It supports only multivariate models (for now ...)

    if length(γ_grid.p) != 2 || length(γ_grid.λ) != 2 || length(γ_grid.α) != 2 || length(γ_grid.β) != 2
        error("The grids include more or less than two entries. They must include only the lower and upper bounds for the grids!")
    end

    check_bounds(γ_grid.p[2], γ_grid.p[1]);
    check_bounds(γ_grid.λ[2], γ_grid.λ[1]);
    check_bounds(γ_grid.α[2], γ_grid.α[1]);
    check_bounds(γ_grid.β[2], γ_grid.β[1]);
    check_bounds(γ_grid.p[1], 1);
    check_bounds(γ_grid.λ[1], 0);
    check_bounds(γ_grid.α[1], 0, 1);
    check_bounds(γ_grid.α[2], 0, 1);
    check_bounds(γ_grid.β[1], 1);

    # Memory pre-allocation for random search
    rng = StableRNG(seed);
    errors = zeros(γ_grid.draws);
    candidates = zeros(4, γ_grid.draws);

    for draw=1:γ_grid.draws
        candidates[:,draw] = [rand(rng, γ_grid.p[1]:γ_grid.p[2]),
                              rand(rng, Uniform(γ_grid.λ[1], γ_grid.λ[2])),
                              rand(rng, Uniform(γ_grid.α[1], γ_grid.α[2])),
                              rand(rng, Uniform(γ_grid.β[1], γ_grid.β[2]))];
    end

    # Setup logger path when `validation_settings.log_folder_path` is specified
    if (validation_settings.verb == true) && ~isnothing(validation_settings.log_folder_path)
        io = open("$(validation_settings.log_folder_path)/status_err_type_$(validation_settings.err_type).txt", "w+");
        global_logger(ConsoleLogger(io));
    end

    # Generate partitions for the block jackknife out-of-sample
    if validation_settings.err_type == 3
        jackknife_data = block_jackknife(validation_settings.Y, validation_settings.subsample);

    # Generate partitions for the artificial jackknife
    elseif validation_settings.err_type == 4
        jackknife_data = artificial_jackknife(validation_settings.Y, validation_settings.subsample, validation_settings.max_samples);
    end

    for iter=1:γ_grid.draws

        # Retrieve candidate hyperparameters
        p, λ, α, β = candidates[:,iter];
        p = Int64(p);

        # Update log
        if validation_settings.verb == true
            @info "$(round(now(), Dates.Second(1))) select_hyperparameters (error estimator $(validation_settings.err_type)) > running iteration $iter (out of $(γ_grid.draws)), γ=($(round(p,digits=3)), $(round(λ,digits=3)), $(round(α,digits=3)), $(round(β,digits=3)))";
            if ~isnothing(validation_settings.log_folder_path)
                flush(io);
            end
        end

        #=
        Some candidate hyperparameters might lead to numerical instabilities. This generally happens when:
        - the candidates are extreme with respect to the model of interest
        - validation_settings.subsample is high relative to the sample size and there is not enough shrinkage

        In simulation, this issue was mostly observed with the block-jackknife.

        A priori, there is not an obvious way to construct a grid of candidates that does not result in errors in the ECM.
        The fc_err function handles these numerical issues on a case-by-case basis.
        =#

        # In-sample or standard out-of-sample
        if validation_settings.err_type <= 2
            errors[iter], inactive_sample = fc_err(validation_settings, p, λ, α, β);
            if inactive_sample == 1
                error("The estimation / validation sample is a matrix of missings!");
            end

        # Jackknife out-of-sample
        else
            errors[iter] = jackknife_err(validation_settings, jackknife_data, p, λ, α, β);
        end
    end

    # Close log
    if (validation_settings.verb == true) && ~isnothing(validation_settings.log_folder_path)
        close(io);
    end

    # Return output
    return candidates, errors;
end


text:
using JuMP
struct DummyOptimizer <: MOI.AbstractOptimizer end
MOI.is_empty(::DummyOptimizer) = true
MOI.supports_constraint(::DummyOptimizer, ::Type{<:MOI.AbstractFunction}, ::Type{<:MOI.AbstractSet}) = true
function MOI.optimize!(::DummyOptimizer, model::MOI.ModelLike)
    return MOI.Utilities.identity_index_map(model), false
end
MOI.get(::DummyOptimizer, ::MOI.TerminationStatus) = MOI.OTHER_ERROR
MOI.get(::DummyOptimizer, ::MOI.RawStatusString) = "Dummy status"

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
# Tests for score filter

using ForwardDiff

δ = 1.0
φ = 0.95
ση = 1.0
σε = 1.0

θ0 = [δ;atanh(φ);ση;σε]

β = [0.0;1.0]

N = 100

@testset "Score filter" begin
    arh = ARH(β)

    X,Y = Filtering.simulate(arh,θ0,N,fill([1.0],N))

    S = ForwardDiff.gradient(θ->-filter_likelihood(kalman_filter(arh,θ,Y,fill([1.0],N))),θ0)

    sf1 = score_filter(arh,θ0,Y,100,fill([1.0],N))
    sf2 = score_filter(arh,θ0,Y,1000,fill([1.0],N))
    sf3 = score_filter(arh,θ0,Y,10000,fill([1.0],N))

    @test loglikelihood(sf1) < 0.0
    @test size(score(sf1),1) == length(θ0)

    # This might not always succeed, but should be true
    @test mean(abs2,S .- score(sf1)) > mean(abs2,S .- score(sf2)) > mean(abs2,S .- score(sf3))
end


text:
module ContinuousBenchmarks

__precompile__(false)

using JSON

export @benchmarkd,
    set_benchmark_config_file,
    get_benchmark_files

include("config.jl")

const BENCH_DIR = abspath(joinpath(@__DIR__, "..", "benchmarks"))
const BENCHMARK_CONFIG_FILE = Base.RefValue{Union{Nothing, String}}(nothing)
CONFIG = Dict()

function set_benchmark_config_file(source::String)
    BENCHMARK_CONFIG_FILE[] = source
    try
        source = BENCHMARK_CONFIG_FILE[]
        eval(Meta.parse("""include("$source")"""))
    catch ex
        @warn "Failed to load config file: $(BENCHMARK_CONFIG_FILE[])"
    end
    for item in CONFIG
        Config.set_config(item.first, item.second)
    end
end

function get_benchmark_files()
    if Config.get_config("benchmark_files") != nothing
        return Config.get_config("benchmark_files")
    else
        bm_files = ["example.jl"]
        return map(bm_files) do x joinpath(BENCH_DIR, x) end
    end
end

## Run benchmark
macro benchmarkd(name, expr)
    quote
        result, t_elapsed, mem, gctime, memallocs  = @timed $(esc(expr))
        $(string(name)), result, t_elapsed, mem, gctime, memallocs
    end
end

include("utils.jl")
include("turing-tools.jl")
include("reporter.jl")
include("runner.jl")
include("AppServer.jl")

end #module


text:
@theme DefaultTheme Dict(
    :name => "Default",
    :description => "Based on the official Julia colours.",
    :style => S"bg: fbfbfb; fg: 444",
    :tokens => Dict(
        TEXT => S"",
        KEYWORD => S"fg: 945bb0; bold",
        KEYWORD_CONSTANT => S"fg: 3b972e; italic",
        KEYWORD_DECLARATION => S"fg: d66661; italic",
        STRING => S"fg: c93d39",
        STRING_ESCAPE => S"fg: 3b972e",
        STRING_INTERPOL => S"",
        LITERAL_DATE => S"fg: 945bb0; italic",
        COMMENT => S"fg: 999977; italic",
        NUMBER => S"fg: 3b972e",
        NAME => S"",
        NAME_FUNCTION => S"fg: 4266d5",
        NAME_DECORATOR => S"fg: d66661",
        OPERATOR => S"bold; fg: 666666",
    )
)


text:
using AlphaZero
using Test

include("../games/tictactoe/main.jl")
include("../games/connect-four/main.jl")
include("../scripts/lib/dummy_run.jl")

include("test_game.jl")

@testset "Testing Games" begin
  test_game(Tictactoe.Game)
  test_game(ConnectFour.Game)
  @test true
end

@testset "Dummy Runs" begin
  dir = "sessions/test-tictactoe"
  @test dummy_run(Tictactoe, nostdout=false, session_dir=dir)
  #@test dummy_run(ConnectFour) # Takes a bit too long for Travis
end


text:
ENV["PLOTS_DEFAULT_BACKEND"] = "Plotly"
data = readdlm("../data/rvs_0001.txt");

want_to_plot_pgram = false
if want_to_plot_pgram
  if !isdefined(:Plots) using Plots end
  if !isdefined(:LombScargle) using LombScargle end
  plot(periodpower(pgram),xscale=:log)
end

include("../src/find_peaks.jl")
num_planets = 2
#=
periods_to_try = find_n_peaks(data[:,1],data[:,2],data[:,3],num_peaks=num_planets)

try_extra_period_near_rotation = false
if try_extra_period_near_rotation
  rotation_period = 20.0
  push!(periods_to_try,rotation_period)
  num_planets += 1
end
=#

num_planets = 0
include("../src/rv_model.jl")    # provides logtarget
import RvModelCircularIndepNoise
RvModelCircularIndepNoise.set_times(data[:,1]);
RvModelCircularIndepNoise.set_obs(data[:,2]);
RvModelCircularIndepNoise.set_sigma_obs(data[:,3]);
num_offsets = 1
param_init_circ_indep = zeros(num_planets*RvModelCircularIndepNoise.num_param_per_planet+num_offsets+RvModelCircularIndepNoise.num_noise_param())
for p in 1:num_planets
  amp = p==1 ? std(data[:,2]) : 0.
  RvModelCircularIndepNoise.set_PKM0(param_init_circ_indep,periods_to_try[p], amp, 1*pi, plid=p)
end
RvModelCircularIndepNoise.set_rvoffset(param_init_circ_indep,0.)
RvModelCircularIndepNoise.set_jitter(param_init_circ_indep, 1.0)

RvModelCircularIndepNoise.logtarget(param_init_circ_indep)

if !isdefined(:QuadGK) using QuadGK end
f(x) = RvModelCircularIndepNoise.logtarget([x])
QuadGK.quadgk(f,0.0,99.0)

if !isdefined(:Optim) using Optim end

# First optimize just phase
planet_to_optimize = 1
function negative_logtarget_fixed_PK(p::Float64)
  RvModelCircularIndepNoise.set_mean_anomaly_at_t0(param_init_circ_indep,p,plid=planet_to_optimize)
  RvModelCircularIndepNoise.negative_logtarget(param_init_circ_indep)
end

result_circ_phase_indep = optimize(negative_logtarget_fixed_PK, 0, 2pi)
result_circ_phase_indep.minimizer
result_circ_phase_indep.minimum/size(data,1)
RvModelCircularIndepNoise.set_mean_anomaly_at_t0(param_init_circ_indep,mod2pi(result_circ_phase_indep.minimizer[1]))

for planet_to_optimize in 2:num_planets
  amp = std(map(i->data[i,2]- RvModelCircularIndepNoise.calc_model_rv(param_init_circ_indep,data[i,1]), 1:size(data,1)))
  RvModelCircularIndepNoise.set_PKM0(param_init_circ_indep,periods_to_try[planet_to_optimize], amp, 1*pi, plid=planet_to_optimize)

  result_circ_phase_indep = optimize(negative_logtarget_fixed_PK, 0, 2pi)
  result_circ_phase_indep.minimizer
  result_circ_phase_indep.minimum/size(data,1)
  RvModelCircularIndepNoise.set_mean_anomaly_at_t0(param_init_circ_indep,mod2pi(result_circ_phase_indep.minimizer[1]),plid=planet_to_optimize)
end

# Now setup model with GP noise
import RvModelCircularGPNoise
RvModelCircularGPNoise.set_times(data[:,1]);
RvModelCircularGPNoise.set_obs(data[:,2]);
RvModelCircularGPNoise.set_sigma_obs(data[:,3]);
num_offsets = 1
param_init_circ_gp = vcat(param_init_circ_indep[1:num_planets*RvModelCircularIndepNoise.num_param_per_planet+num_offsets],zeros(RvModelCircularGPNoise.num_noise_param()) )
RvModelCircularGPNoise.set_noise_param(param_init_circ_gp, [3.0, 50., 0.5, 20.])
RvModelCircularGPNoise.logtarget(param_init_circ_gp)

if !isdefined(:QuadGK) using QuadGK end
f(x) = RvModelCircularGPNoise.logtarget([x])
QuadGK.quadgk(f,0.0,99.0)

# Refitting with a GP model isn't really necessary, but I wrote it to make sure things were working
function negative_logtarget_fixed_PK(p::Float64)
  RvModelCircularGPNoise.set_mean_anomaly_at_t0(param_init_circ_gp,p)
  RvModelCircularGPNoise.negative_logtarget(param_init_circ_gp)
end

if !isdefined(:Optim) using Optim end
result_circ_phase_gp = optimize(negative_logtarget_fixed_PK, result_circ_phase_indep.minimizer-pi, result_circ_phase_indep.minimizer+pi)
result_circ_phase_gp.minimizer
result_circ_phase_gp.minimum/size(data,1)
RvModelCircularGPNoise.set_mean_anomaly_at_t0(param_init_circ_gp,mod2pi(result_circ_phase_gp.minimizer[1]))

# This optimization is slow and still may not complete
# result_circ = optimize(RvModelCircularGPNoise.negative_logtarget, param_init_circ_gp, ConjugateGradient())

plot_model = false
if plot_model
  pred_obs = map(t->RvModelCircularGPNoise.calc_model_rv(param_init_circ_gp,t),data[:,1]);
  times_plot = linspace(minimum(data[:,1]),maximum(data[:,1]),800);
  pred_plot = map(t->RvModelCircularGPNoise.calc_model_rv(param_init_circ_gp,t),times_plot);

  if !isdefined(:Plots) using Plots end
  plot(times_plot,pred_plot)
  scatter!(data[:,1],pred_obs)
end

plot_residuals_periodogram = false
if plot_residuals_periodogram
  residuals = map(i->data[i,2]- RvModelCircularGPNoise.calc_model_rv(param_init_circ_gp,data[i,1]), 1:size(data,1));
  if !isdefined(:LombScargle) using LombScargle end
  resid_pgram = lombscargle(data[:,1],residuals,data[:,3])
  if !isdefined(:Plots) using Plots end
  plot(periodpower(resid_pgram),xscale=:log)
end

include("../src/demcmc.jl")

perturb_scale = fill(0.01,length(param_init_circ_gp))
for p in 1:num_planets
  perturb_scale[1+RvModelCircularGPNoise.pl_offset(p)] = 1e-3
end
pop_init_circ = param_init_circ_gp .+ perturb_scale.*randn(length(param_init_circ_gp),3*length(param_init_circ_gp))
burnin_circ = run_demcmc( pop_init_circ, RvModelCircularGPNoise.logtarget, num_gen = 100, epsilon = 0.01 );

#result_circ = copy(burnin_circ)
result_circ = run_demcmc( burnin_circ["theta_last"], RvModelCircularGPNoise.logtarget, num_gen = 1000, epsilon = 0.01 );

for p in 1:num_planets
  tmp = RvModelCircularGPNoise.inv_transform_period.(result_circ["theta_all"][1+RvModelCircularGPNoise.pl_offset(p),:,:])
  println("param P_",p,": ",mean(tmp), " +/- ", std(tmp))
  tmp = RvModelCircularGPNoise.inv_transform_amplitude.(result_circ["theta_all"][2+RvModelCircularGPNoise.pl_offset(p),:,:])
  println("param K_",p,": ",mean(tmp), " +/- ", std(tmp))
  tmp = result_circ["theta_all"][3+RvModelCircularGPNoise.pl_offset(p),:,:]
  println("param M_",p,": ",mean(tmp), " +/- ", std(tmp))
end
for i in (num_planets*RvModelCircularGPNoise.num_param_per_planet+1):length(param_init_circ_gp)
  println("param #",i,": ",mean(result_circ["theta_all"][i,:,:]), " +/- ", std(result_circ["theta_all"][i,:,:]))
end

# Now with Keplerian model & GP noise
import RvModelKeplerianGPNoise
RvModelKeplerianGPNoise.set_times(data[:,1]);
RvModelKeplerianGPNoise.set_obs(data[:,2]);
RvModelKeplerianGPNoise.set_sigma_obs(data[:,3]);
num_offsets = 1
param_init_kepler_gp = zeros(num_planets*RvModelKeplerianGPNoise.num_param_per_planet+num_offsets+RvModelKeplerianGPNoise.num_noise_param())
for p in 1:num_planets
  (P,K,M0) = RvModelCircularGPNoise.extract_PKM0(param_init_circ_gp, plid=p)
  e = 0.01
  w = 0.0
  RvModelKeplerianGPNoise.set_PKewM0(param_init_kepler_gp, P, K, e, w, M0; plid=p)
end
RvModelKeplerianGPNoise.set_rvoffset(param_init_kepler_gp, RvModelCircularGPNoise.extract_rvoffset(param_init_circ_gp) )
RvModelKeplerianGPNoise.set_noise_param(param_init_kepler_gp, [3.0, 50., 0.5, 20.])
RvModelKeplerianGPNoise.logtarget(param_init_kepler_gp)

if !isdefined(:QuadGK) using QuadGK end
f(x) = RvModelKeplerianGPNoise.logtarget([x])
QuadGK.quadgk(f,0.0,99.0)

perturb_scale = fill(0.01,length(param_init_kepler_gp))
for p in 1:num_planets
  perturb_scale[1+RvModelCircularGPNoise.pl_offset(p)] = 1e-3
end
pop_init_kepler = param_init_kepler_gp .+ perturb_scale.*randn(length(param_init_kepler_gp),3*length(param_init_kepler_gp))
burnin_kepler = run_demcmc( pop_init_kepler, RvModelKeplerianGPNoise.logtarget, num_gen = 100, epsilon = 0.01 );

result_kepler = copy(burnin_kepler)
#result_kepler = run_demcmc( burnin_kepler["theta_last"], RvModelKeplerianGPNoise.logtarget, num_gen = 1000, epsilon = 0.01 );

for p in 1:num_planets
  tmp = RvModelKeplerianGPNoise.inv_transform_period.(result_kepler["theta_all"][1+RvModelKeplerianGPNoise.pl_offset(p),:,:])
  println("param P_",p,": ",mean(tmp), " +/- ", std(tmp))
  tmp = RvModelKeplerianGPNoise.inv_transform_amplitude.(result_kepler["theta_all"][2+RvModelKeplerianGPNoise.pl_offset(p),:,:])
  println("param K_",p,": ",mean(tmp), " +/- ", std(tmp))
  tmp = result_kepler["theta_all"][3+RvModelKeplerianGPNoise.pl_offset(p),:,:]
  println("param ecosw_",p,": ",mean(tmp), " +/- ", std(tmp))
  tmp = result_kepler["theta_all"][4+RvModelKeplerianGPNoise.pl_offset(p),:,:]
  println("param esinw_",p,": ",mean(tmp), " +/- ", std(tmp))
  tmp = result_kepler["theta_all"][5+RvModelKeplerianGPNoise.pl_offset(p),:,:]
  println("param M_",p,": ",mean(tmp), " +/- ", std(tmp))
end
for i in (num_planets*RvModelKeplerianGPNoise.num_param_per_planet+1):length(param_init_kepler_gp)
  println("param #",i,": ",mean(result_kepler["theta_all"][i,:,:]), " +/- ", std(result_kepler["theta_all"][i,:,:]))
end


text:
@testset "Test count word occurrence in string" begin
    @test count_word_occurrence_in_string("AAA BBB", "AAA") == 1
    @test count_word_occurrence_in_string("AAA AAA", "AAA") == 2
    @test count_word_occurrence_in_string("AAAAA", "AAA") == 1
end


text:
function internal_wave_solution(; L, background_stratification=false)
    # Internal wave parameters
     ν = κ = 1e-9
    z₀ = -L/3
     δ = L/20
    a₀ = 1e-3
     m = 16
     k = 1
     f = 0.2
     ℕ = 1.0
     σ = sqrt( (ℕ^2*k^2 + f^2*m^2) / (k^2 + m^2) )

    # Numerical parameters
    Δt = 0.01 * 1/σ

    cᵍ = m * σ / (k^2 + m^2) * (f^2/σ^2 - 1)
     U = a₀ * k * σ   / (σ^2 - f^2)
     V = a₀ * k * f   / (σ^2 - f^2)
     W = a₀ * m * σ   / (σ^2 - ℕ^2)
     B = a₀ * m * ℕ^2 / (σ^2 - ℕ^2)

    a(x, y, z, t) = exp( -(z - cᵍ*t - z₀)^2 / (2*δ)^2 )

    u(x, y, z, t) = a(x, y, z, t) * U * cos(k*x + m*z - σ*t)
    v(x, y, z, t) = a(x, y, z, t) * V * sin(k*x + m*z - σ*t)
    w(x, y, z, t) = a(x, y, z, t) * W * cos(k*x + m*z - σ*t)

    # Buoyancy field depends on whether we use background fields or not
    wavy_b(x, y, z, t) = a(x, y, z, t) * B * sin(k*x + m*z - σ*t)
    background_b(x, y, z, t) = ℕ^2 * z

    if background_stratification # Move stratification to a background field
        background_fields = (; b=background_b)
        b = wavy_b
    else
        background_fields = NamedTuple()
        b(x, y, z, t) = wavy_b(x, y, z, t) + background_b(x, y, z, t)
    end

    u₀(x, y, z) = u(x, y, z, 0)
    v₀(x, y, z) = v(x, y, z, 0)
    w₀(x, y, z) = w(x, y, z, 0)
    b₀(x, y, z) = b(x, y, z, 0)

    solution = (; u, v, w, b)

    # Form model keyword arguments
    closure = IsotropicDiffusivity(ν=ν, κ=κ)
    buoyancy = BuoyancyTracer()
    tracers = :b
    coriolis = FPlane(f=f)
    model_kwargs = (; closure, buoyancy, tracers, coriolis)

    return solution, model_kwargs, background_fields, Δt, σ
end

function internal_wave_dynamics_test(model, solution, Δt)
    # Make initial conditions
    u₀(x, y, z) = solution.u(x, y, z, 0)
    v₀(x, y, z) = solution.v(x, y, z, 0)
    w₀(x, y, z) = solution.w(x, y, z, 0)
    b₀(x, y, z) = solution.b(x, y, z, 0)

    set!(model, u=u₀, v=v₀, w=w₀, b=b₀)

    simulation = Simulation(model, stop_iteration=10, Δt=Δt)
    try run!(simulation); catch; end # so the test continues to execute if there's a NaN

    # Tolerance was found by trial and error...
    @test relative_error(model.velocities.u, solution.u, model.clock.time) < 1e-4

    return nothing
end


text:
using TestImages, Images
function voronoi_img!(img, n_centroids)
    n,m = size(img)
    w   = minimum([n,m])
    dist = (point,vector) -> sqrt.((point[1].-vector[:,1]).^2 .+ (point[2].-vector[:,2]).^2)
    dots = [rand(1:n, n_centroids) rand(1:m, n_centroids)]
    c = []
    for i in 1:size(dots,1)
        p = dots[i,:]
        append!(c, [img[p[1],p[2]]])
    end
    dots = [dots c]

    for x in 1:n, y in 1:m
        distances = dist([x,y],dots) # distance
        nn = findmin(distances)[2]
        img[x,y]  = dots[nn,:][3]
    end
end
img = testimage("mandrill")
voronoi_img!(img, 300)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# test CommonConvWrapper

using Test
using NLsolve
using IncrementalInference
using Statistics

import IncrementalInference: getSample

mutable struct FunctorWorks
  a::Array{Float64,2}
end

@testset "FunctorWorks" begin

function (fw::FunctorWorks)(x)
  fw.a[1,1] = -1.0
  nothing
end

A = rand(2,3)
At = deepcopy(A)
At[1,1] = -1.0

fvar = FunctorWorks(A)
fvar(0.0)
@test At == A

end

mutable struct FunctorArray{T}
  # This was a tuple, but array will like work better in the long term
  fnc!::Function
  a::Array{T, 1}
end

##

@testset "FunctorArray" begin

function testarray!(a1::Array{Float64, 2}, a2::Array{Float64,2})
  a1[1,1] = -1.0
  @show a1
  nothing
end

function (fi::FunctorArray)(x)
  fi.fnc!(fi.a...)
end

A = rand(2,3)
B = rand(2,3)
t = Array{Array{Float64,2},1}()
push!(t,A)
push!(t,B)
At = deepcopy(A)
At[1,1] = -1.0

fvar = FunctorArray(testarray!, t)
fvar([0.0])
@test At == A

end

##

# abstract Nonparametric <: Function
# This is what the internmediate user would be contributing
mutable struct Pose1Pose1Test{T} <: AbstractRelativeRoots
  Dx::T
end

getSample(cf::CalcFactor{<:Pose1Pose1Test}, N::Int=1) = (reshape(rand(cf.factor.Dx,N),1,N),)

#proposed standardized parameter list, does not have to be functor
function (cf::CalcFactor{<:Pose1Pose1Test})(Dx,
                                            p1,
                                            p2 )
  #
  return Dx[1] - (p2[1] - p1[1])
end

##

@testset "Test with CommonConvWrapper for un-permuted root finding..." begin

##

N = 110
p1 = rand(1,N)
p2 = rand(1,N)
t = Array{Array{Float64,2},1}()
push!(t,p1)
push!(t,p2)

odo = Pose1Pose1Test(Normal(100.0,1.0))

fg = initfg()
X0 = addVariable!(fg, :x0, ContinuousEuclid{1})
initManual!(fg, :x0, zeros(1,100))
X1 = addVariable!(fg, :x1, ContinuousEuclid{1})
addFactor!(fg, [:x0;:x1], odo, graphinit=false)

pts = approxConv(fg, getFactor(fg, :x0x1f1), :x1)

## Now check the contents of internal CCW

ccw = IIF._getCCW(fg, :x0x1f1)

@test 90.0 < Statistics.mean(ccw.params[ccw.varidx]) < 110.0
@test -10.0 < Statistics.mean(ccw.params[1]) < 10.0

##

println("and in the reverse direction")

initManual!(fg, :x1, 100*ones(1,100))

pts = approxConv(fg, getFactor(fg, :x0x1f1), :x0)

@test -10.0 < Statistics.mean(ccw.params[1]) < 10.0
@test 90.0 < Statistics.mean(ccw.params[2]) < 110.0

##

end

# use the range only example, should give a circle with nothing in the middle

@testset "Generic convolution testing in factor graph context..." begin

##

N=100
p1 = randn(1,N)
d1 = kde!(p1)
p2 = randn(1,N)
t = Array{Array{Float64,2},1}()
push!(t,p1)
push!(t,p2)

fg = initfg()
# fg.registeredModuleFunctions[:Main] = getSample

v1=addVariable!(fg, :x1, ContinuousScalar, N=N)
v2=addVariable!(fg, :x2, ContinuousScalar, N=N)
bws = getBW(d1)[:,1]
f1 = addFactor!(fg, [v1], Prior(kde!(p1, bws)) )

odo = Pose1Pose1Test(Normal(100.0,1.0))
f2 = addFactor!(fg, [v1;v2], odo)

tree = buildTreeReset!(fg)

pts = getVal(getVariable(fg,:x1))
@test abs(Statistics.mean(pts)-0.0) < 10.0
pts = getVal(getVariable(fg,:x2))
@test abs(Statistics.mean(pts)-0.0) < 10.0

# inferOverTreeR!(fg, tree, N=N)
# inferOverTreeR!(fg, tree)
# @time [inferOverTreeR!(fg, tree, N=N) for i in 1:3];
tree, smt, hist = solveTree!(fg)

# using Gadfly
# plot(y=rand(10))
# plotKDE(getBelief(fg,:x2))

pts = getVal(getVariable(fg,:x1))
@test abs(Statistics.mean(pts)-0.0) < 10.0

pts = getVal(getVariable(fg,:x2))
@test abs(Statistics.mean(pts)-100.0) < 10.0

##

end