text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

@doc raw"""
    AcousticPerturbationEquations2D(v_mean_global, c_mean_global, rho_mean_global)

Acoustic perturbation equations (APE) in two space dimensions. The equations are given by
```math
\begin{aligned}
  \frac{\partial\mathbf{v'}}{\partial t} + \nabla (\bar{\mathbf{v}}\cdot\mathbf{v'})
    + \nabla\left( \frac{\bar{c}^2 \tilde{p}'}{\bar{\rho}} \right) &= 0 \\
  \frac{\partial \tilde{p}'}{\partial t} +
    \nabla\cdot (\bar{\rho} \mathbf{v'} + \bar{\mathbf{v}} \tilde{p}') &= 0.
\end{aligned}
```
The bar ``\bar{(\cdot)}`` indicates time-averaged quantities. The unknowns of the APE are the
perturbed velocities ``\mathbf{v'} = (v_1', v_2')^T`` and the scaled perturbed pressure
``\tilde{p}' = \frac{p'}{\bar{c}^2}``, where ``p'`` denotes the perturbed pressure and the
perturbed variables are defined by ``\phi' = \phi - \bar{\phi}``.

In addition to the unknowns, Trixi currently stores the mean values in the state vector,
i.e. the state vector used internally is given by
```math
\mathbf{u} =
  \begin{pmatrix}
    v_1' \\ v_2' \\ \tilde{p}' \\ \bar{v}_1 \\ \bar{v}_2 \\ \bar{c} \\ \bar{\rho}
  \end{pmatrix}.
```
This affects the implementation and use of these equations in various ways:
* The flux values corresponding to the mean values must be zero.
* The mean values have to be considered when defining initial conditions, boundary conditions or
  source terms.
* [`AnalysisCallback`](@ref) analyzes these variables too.
* Trixi's visualization tools will visualize the mean values by default.

The constructor accepts a 2-tuple `v_mean_global` and scalars `c_mean_global` and `rho_mean_global`
which can be used to make the definition of initial conditions for problems with constant mean flow
more flexible. These values are ignored if the mean values are defined internally in an initial
condition.

The equations are based on the APE-4 system introduced in the following paper:
- Roland Ewert and Wolfgang Schröder (2003)
  Acoustic perturbation equations based on flow decomposition via source filtering
  [DOI: 10.1016/S0021-9991(03)00168-2](https://doi.org/10.1016/S0021-9991(03)00168-2)
"""
struct AcousticPerturbationEquations2D{RealT<:Real} <: AbstractAcousticPerturbationEquations{2, 7}
  v_mean_global::SVector{2, RealT}
  c_mean_global::RealT
  rho_mean_global::RealT
end

function AcousticPerturbationEquations2D(v_mean_global::NTuple{2,<:Real}, c_mean_global::Real,
                                         rho_mean_global::Real)
  return AcousticPerturbationEquations2D(SVector(v_mean_global), c_mean_global, rho_mean_global)
end

function AcousticPerturbationEquations2D(; v_mean_global::NTuple{2,<:Real}, c_mean_global::Real,
                                         rho_mean_global::Real)
  return AcousticPerturbationEquations2D(SVector(v_mean_global), c_mean_global, rho_mean_global)
end

varnames(::typeof(cons2cons), ::AcousticPerturbationEquations2D) = ("v1_prime", "v2_prime", "p_prime_scaled",
                                                                    "v1_mean", "v2_mean", "c_mean", "rho_mean")
varnames(::typeof(cons2prim), ::AcousticPerturbationEquations2D) = ("v1_prime", "v2_prime", "p_prime",
                                                                    "v1_mean", "v2_mean", "c_mean", "rho_mean")

# Convenience functions for retrieving state variables and mean variables
function cons2state(u, equations::AcousticPerturbationEquations2D)
  return SVector(u[1], u[2], u[3])
end

function cons2mean(u, equations::AcousticPerturbationEquations2D)
  return SVector(u[4], u[5], u[6], u[7])
end

varnames(::typeof(cons2state), ::AcousticPerturbationEquations2D) = ("v1_prime", "v2_prime", "p_prime_scaled")
varnames(::typeof(cons2mean), ::AcousticPerturbationEquations2D) = ("v1_mean", "v2_mean", "c_mean", "rho_mean")

"""
    global_mean_vars(equations::AcousticPerturbationEquations2D)

Returns the global mean variables stored in `equations`. This makes it easier to define flexible
initial conditions for problems with constant mean flow.
"""
function global_mean_vars(equations::AcousticPerturbationEquations2D)
  return equations.v_mean_global[1], equations.v_mean_global[2], equations.c_mean_global,
         equations.rho_mean_global
end

"""
    initial_condition_constant(x, t, equations::AcousticPerturbationEquations2D)

A constant initial condition where the state variables are zero and the mean flow is constant.
Uses the global mean values from `equations`.
"""
function initial_condition_constant(x, t, equations::AcousticPerturbationEquations2D)
  v1_prime = 0.0
  v2_prime = 0.0
  p_prime_scaled = 0.0

  return SVector(v1_prime, v2_prime, p_prime_scaled, global_mean_vars(equations)...)
end

"""
    initial_condition_convergence_test(x, t, equations::AcousticPerturbationEquations2D)

A smooth initial condition used for convergence tests in combination with
[`source_terms_convergence_test`](@ref). Uses the global mean values from `equations`.
"""
function initial_condition_convergence_test(x, t, equations::AcousticPerturbationEquations2D)
  c = 2.0
  A = 0.2
  L = 2.0
  f = 2.0 / L
  a = 1.0
  omega = 2 * pi * f
  init = c + A * sin(omega * (x[1] + x[2] - a*t))

  v1_prime = init
  v2_prime = init
  p_prime = init^2

  prim = SVector(v1_prime, v2_prime, p_prime, global_mean_vars(equations)...)

  return prim2cons(prim, equations)
end

"""
  source_terms_convergence_test(u, x, t, equations::AcousticPerturbationEquations2D)

Source terms used for convergence tests in combination with
[`initial_condition_convergence_test`](@ref).
"""
function source_terms_convergence_test(u, x, t, equations::AcousticPerturbationEquations2D)
  v1_mean, v2_mean, c_mean, rho_mean = cons2mean(u, equations)

  c = 2.0
  A = 0.2
  L = 2.0
  f = 2.0 / L
  a = 1.0
  omega = 2 * pi * f

  si, co = sincos(omega * (x[1] + x[2] - a * t))
  tmp = v1_mean + v2_mean - a

  du1 = du2 = A * omega * co * (2 * c/rho_mean + tmp + 2/rho_mean * A * si)
  du3 = A * omega * co * (2 * c_mean^2 * rho_mean + 2 * c * tmp + 2 * A * tmp * si) / c_mean^2

  du4 = du5 = du6 = du7 = 0.0

  return SVector(du1, du2, du3, du4, du5, du6, du7)
end

"""
    initial_condition_gauss(x, t, equations::AcousticPerturbationEquations2D)

A Gaussian pulse in a constant mean flow. Uses the global mean values from `equations`.
"""
function initial_condition_gauss(x, t, equations::AcousticPerturbationEquations2D)
  v1_prime = 0.0
  v2_prime = 0.0
  p_prime = exp(-4*(x[1]^2 + x[2]^2))

  prim = SVector(v1_prime, v2_prime, p_prime, global_mean_vars(equations)...)

  return prim2cons(prim, equations)
end

"""
    initial_condition_gauss_wall(x, t, equations::AcousticPerturbationEquations2D)

A Gaussian pulse, used in the `gauss_wall` example elixir in combination with
[`boundary_condition_wall`](@ref). Uses the global mean values from `equations`.
"""
function initial_condition_gauss_wall(x, t, equations::AcousticPerturbationEquations2D)
  v1_prime = 0.0
  v2_prime = 0.0
  p_prime = exp(-log(2) * (x[1]^2 + (x[2] - 25)^2) / 25)

  prim = SVector(v1_prime, v2_prime, p_prime, global_mean_vars(equations)...)

  return prim2cons(prim, equations)
end

"""
    boundary_condition_wall(u_inner, orientation, direction, x, t, surface_flux_function,
                            equations::AcousticPerturbationEquations2D)

Boundary conditions for a solid wall.
"""
function boundary_condition_wall(u_inner, orientation, direction, x, t, surface_flux_function,
                                 equations::AcousticPerturbationEquations2D)
  # Boundary state is equal to the inner state except for the perturbed velocity. For boundaries
  # in the -x/+x direction, we multiply the perturbed velocity in the x direction by -1.
  # Similarly, for boundaries in the -y/+y direction, we multiply the perturbed velocity in the
  # y direction by -1
  if direction in (1, 2) # x direction
    u_boundary = SVector(-u_inner[1], u_inner[2], u_inner[3], cons2mean(u_inner, equations)...)
  else # y direction
    u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3], cons2mean(u_inner, equations)...)
  end

  # Calculate boundary flux
  if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
  end

  return flux
end

"""
  initial_condition_monopole(x, t, equations::AcousticPerturbationEquations2D)

Initial condition for the monopole in a boundary layer setup, used in combination with
[`boundary_condition_monopole`](@ref).
"""
function initial_condition_monopole(x, t, equations::AcousticPerturbationEquations2D)
  m = 0.3 # Mach number

  v1_prime = 0.0
  v2_prime = 0.0
  p_prime = 0.0

  v1_mean = x[2] > 1 ? m : m * (2*x[2] - 2*x[2]^2 + x[2]^4)
  v2_mean = 0.0
  c_mean = 1.0
  rho_mean = 1.0

  prim = SVector(v1_prime, v2_prime, p_prime, v1_mean, v2_mean, c_mean, rho_mean)

  return prim2cons(prim, equations)
end

"""
  boundary_condition_monopole(u_inner, orientation, direction, x, t, surface_flux_function,
                              equations::AcousticPerturbationEquations2D)

Boundary condition for a monopole in a boundary layer at the -y boundary, i.e. `direction = 3`.
This will return an error for any other direction. This boundary condition is used in combination
with [`initial_condition_monopole`](@ref).
"""
function boundary_condition_monopole(u_inner, orientation, direction, x, t, surface_flux_function,
                                     equations::AcousticPerturbationEquations2D)
  if direction != 3
    error("expected direction = 3, got $direction instead")
  end

  # Wall at the boundary in -y direction with a monopole at -0.05 <= x <= 0.05. In the monopole area
  # we use a sinusoidal boundary state for the perturbed variables. For the rest of the -y boundary
  # we set the boundary state to the inner state and multiply the perturbed velocity in the
  # y-direction by -1.
  if -0.05 <= x[1] <= 0.05 # Monopole
    v1_prime = 0.0
    v2_prime = p_prime = sin(2 * pi * t)

    prim_boundary = SVector(v1_prime, v2_prime, p_prime, u_inner[4], u_inner[5], u_inner[6], u_inner[7])

    u_boundary = prim2cons(prim_boundary, equations)
  else # Wall
    u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3], u_inner[4], u_inner[5], u_inner[6],
                         u_inner[7])
  end

  # Calculate boundary flux
  flux = surface_flux_function(u_boundary, u_inner, orientation, equations)

  return flux
end

"""
    boundary_condition_zero(u_inner, orientation, direction, x, t, surface_flux_function,
                            equations::AcousticPerturbationEquations2D)

Boundary condition that uses a boundary state where the state variables are zero and the mean
variables are the same as in `u_inner`.
"""
function boundary_condition_zero(u_inner, orientation, direction, x, t, surface_flux_function,
                                 equations::AcousticPerturbationEquations2D)
  value = zero(eltype(u_inner))
  u_boundary = SVector(value, value, value, cons2mean(u_inner, equations)...)

  # Calculate boundary flux
  if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
  end

  return flux
end

"""
    boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function,
                                 equations::AcousticPerturbationEquations2D)

Use an orthogonal projection of the perturbed velocities to zero out the normal velocity
while retaining the possibility of a tangential velocity in the boundary state.
Further details are available in the paper:
- Marcus Bauer, Jürgen Dierke and Roland Ewert (2011)
  Application of a discontinuous Galerkin method to discretize acoustic perturbation equations
  [DOI: 10.2514/1.J050333](https://doi.org/10.2514/1.J050333)
"""
function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, x, t,
                                      surface_flux_function, equations::AcousticPerturbationEquations2D)
  # normalize the outward pointing direction
  normal = normal_direction / norm(normal_direction)

  # compute the normal perturbed velocity
  u_normal = normal[1] * u_inner[1] + normal[2] * u_inner[2]

  # create the "external" boundary solution state
  u_boundary = SVector(u_inner[1] - 2.0 * u_normal * normal[1],
                       u_inner[2] - 2.0 * u_normal * normal[2],
                       u_inner[3], cons2mean(u_inner, equations)...)

  # calculate the boundary flux
  flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)

  return flux
end

# Calculate 1D flux for a single point
@inline function flux(u, orientation::Integer, equations::AcousticPerturbationEquations2D)
  v1_prime, v2_prime, p_prime_scaled = cons2state(u, equations)
  v1_mean, v2_mean, c_mean, rho_mean = cons2mean(u, equations)

  # Calculate flux for conservative state variables
  if orientation == 1
    f1 = v1_mean * v1_prime + v2_mean * v2_prime + c_mean^2 * p_prime_scaled / rho_mean
    f2 = zero(eltype(u))
    f3 = rho_mean * v1_prime + v1_mean * p_prime_scaled
  else
    f1 = zero(eltype(u))
    f2 = v1_mean * v1_prime + v2_mean * v2_prime + c_mean^2 * p_prime_scaled / rho_mean
    f3 = rho_mean * v2_prime + v2_mean * p_prime_scaled
  end

  # The rest of the state variables are actually variable coefficients, hence the flux should be
  # zero. See https://github.com/trixi-framework/Trixi.jl/issues/358#issuecomment-784828762
  # for details.
  f4 = f5 = f6 = f7 = zero(eltype(u))

  return SVector(f1, f2, f3, f4, f5, f6, f7)
end

# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::AcousticPerturbationEquations2D)
  # Calculate v = v_prime + v_mean
  v_prime_ll = u_ll[orientation]
  v_prime_rr = u_rr[orientation]
  v_mean_ll = u_ll[orientation + 3]
  v_mean_rr = u_rr[orientation + 3]

  v_ll = v_prime_ll + v_mean_ll
  v_rr = v_prime_rr + v_mean_rr

  c_mean_ll = u_ll[6]
  c_mean_rr = u_rr[6]

  λ_max = max(abs(v_ll), abs(v_rr)) + max(c_mean_ll, c_mean_rr)
end

# Calculate 1D flux for a single point in the normal direction
# Note, this directional vector is not normalized
@inline function flux(u, normal_direction::AbstractVector, equations::AcousticPerturbationEquations2D)
  v1_prime, v2_prime, p_prime_scaled = cons2state(u, equations)
  v1_mean, v2_mean, c_mean, rho_mean = cons2mean(u, equations)

  f1 = normal_direction[1] * (v1_mean * v1_prime + v2_mean * v2_prime + c_mean^2 * p_prime_scaled / rho_mean)
  f2 = normal_direction[2] * (v1_mean * v1_prime + v2_mean * v2_prime + c_mean^2 * p_prime_scaled / rho_mean)
  f3 = ( normal_direction[1] * (rho_mean * v1_prime + v1_mean * p_prime_scaled)
       + normal_direction[2] * (rho_mean * v2_prime + v2_mean * p_prime_scaled) )

  # The rest of the state variables are actually variable coefficients, hence the flux should be
  # zero. See https://github.com/trixi-framework/Trixi.jl/issues/358#issuecomment-784828762
  # for details.
  f4 = f5 = f6 = f7 = zero(eltype(u))

  return SVector(f1, f2, f3, f4, f5, f6, f7)
end

# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, equations::AcousticPerturbationEquations2D)
  # Calculate v = v_prime + v_mean
  v_prime_ll = normal_direction[1]*u_ll[1] + normal_direction[2]*u_ll[2]
  v_prime_rr = normal_direction[1]*u_rr[1] + normal_direction[2]*u_rr[2]
  v_mean_ll = normal_direction[1]*u_ll[4] + normal_direction[2]*u_ll[5]
  v_mean_rr = normal_direction[1]*u_rr[4] + normal_direction[2]*u_rr[5]

  v_ll = v_prime_ll + v_mean_ll
  v_rr = v_prime_rr + v_mean_rr

  c_mean_ll = u_ll[6]
  c_mean_rr = u_rr[6]

  # The v_normals are already scaled by the norm
  λ_max = max(abs(v_ll), abs(v_rr)) + max(c_mean_ll, c_mean_rr) * norm(normal_direction)
end

# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the mean values
@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction,
                                                              equations::AcousticPerturbationEquations2D)
  λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations)
  diss = -0.5 * λ * (u_rr - u_ll)
  z = zero(eltype(u_ll))
  return SVector(diss[1], diss[2], diss[3], z, z, z, z)
end

@inline have_constant_speed(::AcousticPerturbationEquations2D) = Val(false)

@inline function max_abs_speeds(u, equations::AcousticPerturbationEquations2D)
  v1_mean = u[4]
  v2_mean = u[5]
  c_mean = u[6]

  return abs(v1_mean) + c_mean, abs(v2_mean) + c_mean
end

# Convert conservative variables to primitive
@inline function cons2prim(u, equations::AcousticPerturbationEquations2D)
  p_prime_scaled = u[3]
  c_mean = u[6]
  p_prime = p_prime_scaled * c_mean^2

  return SVector(u[1], u[2], p_prime, u[4], u[5], u[6], u[7])
end

# Convert primitive variables to conservative
@inline function prim2cons(u, equations::AcousticPerturbationEquations2D)
  p_prime = u[3]
  c_mean = u[6]
  p_prime_scaled = p_prime / c_mean^2

  return SVector(u[1], u[2], p_prime_scaled, u[4], u[5], u[6], u[7])
end

# Convert conservative variables to entropy variables
@inline cons2entropy(u, equations::AcousticPerturbationEquations2D) = u

end # @muladd


text:
using Pkg; Pkg.activate(first(splitdir(@__DIR__)))

using CUDA
using DLPack
using PyCall

@pyimport hoomd
@pyimport hoomd.md as md
@pyimport hoomd.dump as dump
@pyimport hoomd.group as group
@pyimport hoomd.benchmark as benchmark

@pyimport numpy as np

N = 10000
mode_string = "gpu"

if mode_string ∉ ("cpu", "gpu")
    throw("""Execution mode argument must be either "cpu" or "gpu".""")
end

# %%
hoomd.context.initialize("--mode=" * mode_string)
rcut = 3.0
sqrt_N = int(sqrt(N))

system = hoomd.init.create_lattice(
    unitcell = hoomd.lattice.sq(a = 2.0),
    n = [sqrt_N, sqrt_N]
)
nlist = md.nlist.cell(check_period = 1)
# basic LJ forces from Hoomd
lj = md.pair.lj(rcut, nlist)
lj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0)
md.integrate.mode_standard(dt=0.005)
md.integrate.langevin(group = group.all(), kT=1.0, seed=42)
# equilibrate for 4k steps first
hoomd.run(4000)

# %%
import pysages

pysages.set_backend("hoomd")
simulation = hoomd.context.current
simulation_view = pysages.wrap(simulation)

# %%
ξ1 = pysages.cvs.distance(np.array([0, 1]))
ξ = pysages.collective_variable(ξ1)

grid = pysages.Grid(
    lower = (0,),
    upper = (2.0,),
    shape = (1024,),
    periodicity = (True,)
)

sampling_method = pysages.methods.abf(simulation_view, grid, ξ, N = 100)

# %%
sampler = pysages.bind(simulation, sampling_method)

# %%
hoomd.run(5000)

# %%
benchmark_results = hoomd.benchmark.series(
    warmup = 6000,
    repeat = 5,
    steps = 50000,
    limit_hours = 2
)

ext: jl


text:
using ADVUMPS
using ADVUMPS: model_tensor, tensorfromclassical
using Test
using Zygote, OMEinsum

@testset "exampletensor" begin
    β = rand()
    @test model_tensor(Ising(),β) ≈ tensorfromclassical([β -β; -β β])
end


text:
# Define TDistLocationScale as the location scale transformation of the Distributions.jl TDist
export TDistLocationScale
TDistLocationScale = LocationScale{Float64,TDist{Float64}}

export BetaFourParameters
BetaFourParameters = LocationScale{Float64,Beta{Float64}}


text:
module LinearInterp
"""
linear interpolation without the hot air
I really need this function, and honestly I don't care that much
about its speed or flexibility. I just need sth that works for my
BL formulation
"""
struct Interp1D
    x::Vector{Real}
    y::Vector{Real}
    m::Vector{Real}
    b::Vector{Real}
end

"""
interp1D(x::Array{T, 1}, y::Array{T, 1}) where T<:Real
Linear interpollation function which returns a Interp1D struct
    ```
    x = collect(sort(rand(10)))
    y = sin.(10*x)
    xfine = collect(x[1]:0.01:x[end])
    intp1 = interp1D(x, y)
    yfine = intp1(xfine)
    ```
"""
function interp1D(x::Array{T, 1}, y::Array{T, 1}) where T<:Real
    m = (y[2:end]-y[1:end-1])./(x[2:end]-x[1:end-1])
    b = -m.*x[1:end-1]+y[1:end-1]
    return Interp1D(x, y, m, b)
end

function evaluate(intrp::Interp1D, xval::Real)
    if xval<intrp.x[1]
        return intrp.y[1]
    elseif xval>=intrp.x[end]
        return intrp.y[end]
    end
    for i in 1:length(intrp.x)-1
        if intrp.x[i]<=xval<intrp.x[i+1]
            return intrp.m[i]*xval+intrp.b[i]
        end
    end
end

function evaluate(intrp::Interp1D, xvec::Array{T, 1}) where T<:Real
    [evaluate(intrp, xvec[i]) for i in 1:length(xvec)]
end

(intrp1d::Interp1D)(x::Real) = evaluate(intrp1d, x)
(intrp1d::Interp1D)(x::Array{T, 1}) where T<:Real = evaluate(intrp1d, x)

end # module

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
module Persistence1D

include("../deps/deps.jl")

struct Persistence{T}
    data::Vector{T}
    minima::Vector{Int}
    maxima::Vector{Int}
    gminindex::Int
    gminvalue::T
    persistence::Vector{T}
end

function find_persistence(v::Array{Float32},
                          threshold = 0,
                          persistence = true)
    min_p = Ref{Ptr{Cint}}(0)
    max_p = Ref{Ptr{Cint}}(0)
    nmin = Ref{Csize_t}(0)
    nmax = Ref{Csize_t}(0)
    gminindex,gminvalue = Ref{Cint}(0),Ref{Cfloat}(0)
    pers_p,npers = persistence ? (Ref{Ptr{Cfloat}}(0),Ref{Csize_t}(1)) : (C_NULL,Ref{Csize_t}(0))

    ccall((:find_extrema, _jl_libpersistence1d),
          Cvoid,
          (Ptr{Cfloat}, Csize_t,
           Ref{Ptr{Cint}}, Ref{Csize_t},
           Ref{Ptr{Cint}}, Ref{Csize_t},
           Ref{Cint}, Ref{Cfloat},
           Ref{Ptr{Cfloat}}, Ref{Csize_t},
           Cfloat),
          v, length(v),
          min_p, nmin,
          max_p, nmax,
          gminindex, gminvalue,
          pers_p, npers,
          threshold)

    nmin = convert(Int64,nmin[])
    nmax = convert(Int64,nmax[])
    min_v = unsafe_wrap(Array,min_p[], nmin, own=true)
    max_v = unsafe_wrap(Array,max_p[], nmax, own=true)

    gminindex,gminvalue = convert(Int32,gminindex[]),convert(Float32,gminvalue[])

    persistence || return [min_v;gminindex],max_v
    npers = convert(Int64,npers[])
    pers_v = unsafe_wrap(Array,pers_p[], npers, own=true)
    Persistence(v, min_v, max_v, gminindex, gminvalue, pers_v)
end

find_persistence(v::Array{T}, threshold = 0, persistence = true) where {T<:Real} =
    find_persistence(convert(Array{Float32}, v), threshold, persistence)

find_extrema(v::Array{T}, threshold = 0) where {T<:Real} =
    find_persistence(v, threshold, false)

import Base.filter
function filter(p::Persistence, threshold::Real)
    pis = find(p.persistence .> threshold)
    Persistence(p.data,
                p.minima[pis],
                p.maxima[pis],
                p.gminindex,
                p.gminvalue,
                p.persistence[pis])
end

export find_persistence, find_extrema, filter

include("reconstruct.jl")

end


text:
"""
    write_data(vars_results::Dict{Symbol, DataFrames.DataFrame}, save_path::AbstractString; kwargs...)

Receives the variables dictionary from the operation model results and writes each variable dataframe
to a file. The default file type is feather.

# Arguments
-`vars_results::Dict{Symbol, DataFrames.DataFrame} = res.variables`: dictionary of DataFrames
-`save_path::AbstractString`: the file path that the variable files should be written to populate.

# Example
```julia
res = solve_op_problem!(op_problem)
write_data(res.variables, "Users/downloads")
```
# Accepted Key Words
-`file_type::String = CSV`: default filetype is Feather, but this key word can be used to make it CSV.
if a different file type is desired the code will have to be changed to accept it.

"""
# writing a dictionary of dataframes to files

function write_data(vars_results::Dict{Symbol, DataFrames.DataFrame}, save_path::AbstractString; kwargs...)
    file_type = get(kwargs, :file_type, Feather)
    if file_type == Feather || file_type == CSV
        for (k,v) in vars_results
            file_path = joinpath(save_path,"$(k).$(lowercase("$file_type"))")
            file_type.write(file_path, vars_results[k])
        end
    else
        error("unsupported file type: $file_type")
    end

    return
end

# writing a dictionary of dataframes to files and appending the time

function write_data(vars_results::Dict{Symbol, DataFrames.DataFrame}, time::DataFrames.DataFrame, save_path::AbstractString; kwargs...)
    file_type = get(kwargs, :file_type, Feather)
    if file_type == Feather || file_type == CSV
        for (k,v) in vars_results
            if size(time,2) == size(v,2)
                var = hcat(time, v)
            else
                var = v
            end
            file_path = joinpath(save_path,"$(k).$(lowercase("$file_type"))")
            file_type.write(file_path, var)
        end
    else
        error("unsupported file type: $file_type")
    end
    return
end

# taking the variable DataFrame and writing them to a featherfile
## TO DO ## is this one being used??
#=
function write_data(results::Dict, save_path::AbstractString, variable::String; kwargs...)
    file_type = get(kwargs, :file_type, Feather)
    if file_type == Feather || file_type == CSV

        file_path = joinpath(save_path,"$(variable).$(lowercase("$file_type"))")
        file_type.write(file_path)

    else
        error("unsupported file type: $file_type")
    end
    return
end
=#
# taking the variable DataFrame and writing them to a featherfile

function write_data(data::DataFrames.DataFrame, save_path::AbstractString, file_name::String; kwargs...)
    if isfile(save_path)
        save_path = dirname(save_path)
    end
    file_type = get(kwargs, :file_type, Feather)
    if file_type == Feather || file_type == CSV
        file_path = joinpath(save_path,"$(file_name).$(lowercase("$file_type"))")
        file_type.write(file_path, data)
    else
        error("unsupported file type: $file_type")
    end
    return
end

function _write_optimizer_log(optimizer_log::Dict, save_path::AbstractString)

    JSON.write(joinpath(save_path, "optimizer_log.json"), JSON.json(optimizer_log))

end

function _write_data(canonical::Canonical, save_path::AbstractString; kwargs...)
    file_type = get(kwargs, :file_type, Feather)
    if file_type == Feather || file_type == CSV
        for (k, v) in get_variables(canonical)
            file_path = joinpath(save_path,"$(k).$(lowercase("$file_type"))")
            file_type.write(file_path, _result_dataframe_variables(v))
        end
    else
        error("unsupported file type: $file_type")
    end
    return
end

function write_data(canonical::Canonical, save_path::AbstractString, dual_con::Vector{Symbol}; kwargs...)
    file_type = get(kwargs, :file_type, Feather)
    if file_type == Feather || file_type == CSV
        duals = get_model_duals(canonical, dual_con)
        for (k, v) in duals
            file_path = joinpath(save_path,"$(k)_dual.$(lowercase("$file_type"))")
            file_type.write(file_path, _result_dataframe_vars(v))
        end
    else
        error("unsupported file type: $file_type")
    end
    return
end

function write_data(op_problem::OperationsProblem, save_path::AbstractString; kwargs...)
    _write_data(op_problem.canonical, save_path; kwargs...)
    return
end

function write_data(stage::_Stage, save_path::AbstractString; kwargs...)
    _write_data(stage.canonical, save_path; kwargs...)
    return
end

# These functions are writing directly to the feather file and skipping printing to memory.
function _export_model_result(stage::_Stage, start_time::Dates.DateTime, save_path::String)
    write_data(stage, save_path)
    write_data(get_time_stamps(stage, start_time), save_path, "time_stamp")
    return
end

function _export_model_result(stage::_Stage, start_time::Dates.DateTime, save_path::String, dual_con::Vector{Symbol})
    write_data(stage, save_path)
    write_data(stage, save_path, dual_con)
    write_data(get_time_stamp(stage, start_time), save_path, "time_stamp")
    return
end

function _export_optimizer_log(optimizer_log::Dict{Symbol, Any},
                               canonical::Canonical,
                               path::String)

    optimizer_log[:obj_value] = JuMP.objective_value(canonical.JuMPmodel)
    optimizer_log[:termination_status] = Int(JuMP.termination_status(canonical.JuMPmodel))
    optimizer_log[:primal_status] = Int(JuMP.primal_status(canonical.JuMPmodel))
    optimizer_log[:dual_status] = Int(JuMP.dual_status(canonical.JuMPmodel))
    try
        optimizer_log[:solve_time] = MOI.get(canonical.JuMPmodel, MOI.SolveTime())
    catch
        @warn("SolveTime() property not supported by the Solver")
        optimizer_log[:solve_time] = NaN #"Not Supported by solver"
    end
    _write_optimizer_log(optimizer_log, path)
    return
end

""" Exports Operational Model Results to a path"""
function write_model_results(results::Results, save_path::String)
    if !isdir(save_path)
        error("Specified path is not valid. Run write_results to save results.")
    end
    new_folder_path = replace_chars("$save_path/$(round(Dates.now(),Dates.Minute))", ":", "-")
    folder_path = mkdir(new_folder_path)
    write_data(results.variables, folder_path)
    _write_optimizer_log(results.optimizer_log, folder_path)
    write_data(results.time_stamp, folder_path, "time_stamp")
    println("Files written to $folder_path folder.")
    return
end
""" Exports Simulation Model Results to a path"""
function write_model_results(res::Results, path::String, results::String; kwargs...)
    if !isdir(path)
        error("Specified path is not valid. Run write_results to save results.")
    end
    folder_path = joinpath(path,results)
    write_data(res.variables, res.time_stamp, folder_path; kwargs...)
    write_data(res.duals, folder_path; kwargs...)
    _write_optimizer_log(res.optimizer_log, folder_path)
    write_data(res.time_stamp, folder_path, "time_stamp"; kwargs...)
    println("Files written to $folder_path folder.")
    return
end

""" Exports the OpModel JuMP object in MathOptFormat"""
function write_op_problem(op_problem::OperationsProblem, save_path::String)
    _write_canonical(op_problem.canonical, save_path)
    return
end

""" Exports the OpModel JuMP object in MathOptFormat"""
function _write_canonical(canonical::Canonical, save_path::String)
    MOF_model = MOPFM()
    MOI.copy_to(MOF_model, JuMP.backend(canonical.JuMPmodel))
    MOI.write_to_file(MOF_model, save_path)
    return
end


text:
module GenerativeModels

    using Reexport
    using Random
    using Flux

    using Distributions
    using DistributionsAD
    @reexport using ConditionalDists
    @reexport using IPMeasures

    using IPMeasures: AbstractKernel
    using ConditionalDists: AbstractConditionalDistribution
    const ACD = AbstractConditionalDistribution
    const TuMvNormal = Union{DistributionsAD.TuringDenseMvNormal,
                             DistributionsAD.TuringDiagMvNormal,
                             DistributionsAD.TuringScalMvNormal}

    export VAE, GAN, VAMP
    export elbo, mmd_mean, mmd_rand, generator_loss, discriminator_loss
    export train!, softplus_safe, save_checkpoint, load_checkpoint
    export mmd_mean_vamp, init_vamp_mean, init_vamp_sample

    export NeuralStatistician

    abstract type AbstractGM end
    abstract type AbstractVAE <: AbstractGM end
    abstract type AbstractGAN <: AbstractGM end

    # to make Flux.gpu work on VAE/GAN/etc priors we need:
    Flux.@functor DistributionsAD.TuringScalMvNormal
    Flux.@functor DistributionsAD.TuringDiagMvNormal
    Flux.@functor DistributionsAD.TuringDenseMvNormal

    include("utils.jl")
    include("vae.jl")
    include("gan.jl")
    include("vamp.jl")
    include("statistician.jl")

end # module


text:
using PaddedViews

"""
    scale(cube::AbstractArray{T,3}, scales[, size])

Linearly stretch each frame in `cube` by the corresponding scale in `scales`. Uses bilinear interpolaiton internally. The output size can be specified or else we will choose the smallest size that contains the largest stretch.
"""
function scale(cube::AbstractArray{T, 3}, scales, size) where T
    out = similar(cube, Base.size(cube, 1), size...)
    Threads.@threads for i in axes(cube, 1)
        frame = @view cube[i, :, :]
        out[i, :, :] .= scale(frame, scales[i], size)
    end
    return out
end

scale(cube::AbstractArray{T, 3}, scales) where {T} = scale(cube, scales, stretch_size(cube, scales))

stretch_size(cube, scales) = stretch_size((size(cube, 2), size(cube, 3)), scales)

function stretch_size(size::Tuple, scales)
    max_size = ceil.(Int, maximum(scales) .* size)
    return check_size(size, max_size)
end

scale(frame::AbstractMatrix, scale) = imresize(frame; ratio=scale)

"""
    scale(frame::AbstractMatrix, scale[, size])

Linearly stretch `frame` with the ratio `scale`, padding symmetrically with zeros to match `out_size`. Uses bilinear interpolation internally. If `size` is specified, will return a padded view.
"""
function scale(frame::AbstractMatrix{T}, scale, out_size) where T
    res_frame = imresize(frame; ratio=scale)
    init_pos = (out_size .- size(res_frame)) .÷ 2 .+ 1
    return PaddedView(zero(T), res_frame, out_size, init_pos)
end

"""
    invscale(cube::AbstractArray{T,3}, scales[, size])

Linearly contract each frame in `cube` by the corresponding scale in `scales`. Uses bilinear interpolaiton internally. The output size can be specified or else we will choose the smallest size that contains the largest stretch.
"""
function invscale(cube::AbstractArray{T, 3}, scales, size) where T
    out = similar(cube, Base.size(cube, 1), size...)
    Threads.@threads for i in axes(cube, 1)
        out[i, :, :] .=  invscale(cube[i, :, :], scales[i], size)
    end
    return out
end

invscale(cube::AbstractArray{T, 3}, scales) where {T} = invscale(cube, scales, shrink_size(cube, scales))

shrink_size(cube, scales) = shrink_size((size(cube, 2), size(cube, 3)), scales)

function shrink_size(size::Tuple, scales)
    min_size = ceil.(Int, size ./ maximum(scales))
    return check_size(size, min_size)
end

"""
    invscale(frame::AbstractMatrix, scale[, size])

Linearly contract `frame` with the ratio `scale`. Uses bilinear interpolation internally. If `size` is specified, will return a crop view.
"""
invscale(frame::AbstractMatrix, scale) = imresize(frame; ratio=inv(scale))
invscale(frame::AbstractMatrix, scale, size) = cropview(imresize(frame; ratio=inv(scale)), size; force=true, verbose=false)

"""
    scale_and_stack(spcube::AbstractArray{T,4}, scales)

Given a 4-D spectral ADI (SDI) tensor, scales each, temporal slice according to `scales` and then concatenates into a cube with `nλ * nf` frames.
"""
function scale_and_stack(spcube::AbstractArray{T, 4}, scales) where T
    nλ, nf, ny, nx = size(spcube)
    frame_size = (ny, nx)
    max_scale = maximum(scales)
    out_size = check_size(frame_size, ceil.(Int, max_scale .* frame_size))
    out = similar(spcube, nλ * nf, out_size...)
    Threads.@threads for n in axes(spcube, 2)
        slice = (n - 1) * nλ + firstindex(spcube, 2):n * nλ
        frame = @view spcube[:, n, :, :]
        out[slice, :, :] .= scale(frame, scales)
    end

    return out
end

"""
    invscale_and_collapse(stack_cube::AbstractArray{T,3}, scales, size; kwargs...)

Given an SDI tensor that has been stacked into a cube, invscales each spectral slice and combines with [`collapse`](@ref). The output will be cropped to `size`.
"""
function invscale_and_collapse(stack_cube::AbstractArray{T, 3}, scales, size; kwargs...) where T
    nλ = length(scales)
    n = Base.size(stack_cube, 1) ÷ nλ
    out = similar(stack_cube, n, size...)
    Threads.@threads for n in axes(out, 1)
        slice = (n - 1) * nλ + firstindex(out, 1):n * nλ
        cube = @view stack_cube[slice, :, :]
        out[n, :, :] .= collapse(invscale(cube, scales, size); kwargs...)
    end
    return out
end

"""
    scale_list(wavelengths)

Returns a list of scaling factors for aligning SDI tensors from a list of wavelengths.
"""
scale_list(wavelengths) = maximum(wavelengths) ./ wavelengths


text:
#Function to calculate Clebsch-Gordan coeficients
#22/01/2016

function ClebschGord(jj1,jj2,JJ,mm1,mm2,MM)
return CG =
    (-1)^(jj1-jj2+MM)*sqrt(2*JJ+1)*Wigner3j(jj1,jj2,JJ,mm1,mm2,-MM)
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
#=
  SEND+MOST=MONEY in Julia + ConstraintSolver.jl

  Alphametic problem were we maximize MONEY.

  This version do two things:
    - find the maximum of MONEY
    - and then find all solutions for the maximum value of MONEY.

  Problem from the lecture notes:
  http://www.ict.kth.se/courses/ID2204/notes/L01.pdf

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
const CS = ConstraintSolver
include("constraints_utils.jl")

function send_most_money()
    model = Model(optimizer_with_attributes(CS.Optimizer,
                                                        # "all_solutions"=>true,
                                                        "all_optimal_solutions"=>true,
                                                        "logging"=>[]))
    @variable(model, 0 <= x[1:8] <= 9, Int)
    s,e,n,d,m,o,t,y = x

    @variable(model, 10000 <= MONEY <= 99999, Int)

    @constraint(model, x in CS.AllDifferentSet())
    @constraint(model, s != 0)
    @constraint(model, m != 0)
    @constraint(model, MONEY == 10000*m + 1000*o + 100*n + 10*e + y)
    @constraint(model,
                              1000*s + 100*e + 10*n + d
                           +  1000*m + 100*o + 10*s + t
                 == MONEY
        )

    @objective(model,Max,MONEY)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        for sol in 1:num_sols
            println("solution #$sol")
            xx = convert.(Integer,JuMP.value.(x; result=sol))
            MONEYx = convert.(Integer,JuMP.value.(MONEY; result=sol))
            println("x:$xx MONEY:$MONEYx\n")
        end

    end
end

send_most_money()

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
@testset "rendering.jl" begin
    img = WordCloud.rendertext("test", 88.3, color="blue", angle=20, border=2)
    mat, svg = WordCloud.rendertext(Random.randstring(rand(1:8)), rand(5:50), angle=rand(0:180), type=:both)
    @test all(WordCloud.tobitmap(svg) .≈ mat)

    for i in 1:20
        @test 0.9 < WordCloud.Render.gamma(i+1) / prod(1:i) < 1.1
    end

    h = 300+300rand()
    w = 300+300rand()
    sh = WordCloud.Render.shape(ellipse, h, w,color=0)
    true_area = WordCloud.occupying(WordCloud.imagemask(WordCloud.tobitmap(sh), (0, 0, 0, 0)), false)
    @test 0.8 < WordCloud.ellipse_area(h, w) / true_area < 1.2

    h = 300+300rand()
    w = 300+300rand()
    r = 10 + 140rand()
    sh = WordCloud.Render.shape(box, h, w, cornerradius=r,color=0)
    true_area = WordCloud.occupying(WordCloud.imagemask(WordCloud.tobitmap(sh), (0, 0, 0, 0)), false)
    @test 0.8 < WordCloud.box_area(h, w, cornerradius=r) / true_area < 1.2

    h = 300+300rand()
    w = 300+300rand()
    rt = 3rand()
    sh = WordCloud.Render.shape(squircle, h, w, rt=rt,color=0)
    true_area = WordCloud.occupying(WordCloud.imagemask(WordCloud.tobitmap(sh), (0, 0, 0, 0)), false)
    @test 0.8 < WordCloud.squircle_area(h, w, rt=rt) / true_area < 1.2

    h = 300+300rand()
    w = 300+300rand()
    npoints = rand(3:10)
    sh = WordCloud.Render.shape(ngon, h, w, npoints=npoints, color=0)
    true_area = WordCloud.occupying(WordCloud.imagemask(WordCloud.tobitmap(sh), (0, 0, 0, 0)), false)
    @test 0.8 < WordCloud.ngon_area(h, w, npoints=npoints) / true_area < 1.2

    h = 300+300rand()
    w = 300+300rand()
    npoints = rand(3:10)
    starratio = 0.3 + 0.7rand()
    sh = WordCloud.Render.shape(star, h, w, npoints=npoints, starratio=starratio, color=0)
    true_area = WordCloud.occupying(WordCloud.imagemask(WordCloud.tobitmap(sh), (0, 0, 0, 0)), false)
    WordCloud.star_area(h, w, npoints=npoints, starratio=starratio) / true_area
end


text:
module OmegaMH

using InferenceBase
# using TransformVariables

include("mhcore.jl")
# include("proposals.jl")
include("ssproposal.jl")
include("interface.jl")

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text: @cdef :PyInstanceMethod_New PyPtr (PyPtr,)


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Johannes Stoljar
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

################################## INSTALLATION ####################################################
# (1) Enter Package Manager via      ]
# (2) Install FMI via                add FMI       or   add "https://github.com/ThummeTo/FMI.jl"
# (3) Install FMIFlux via            add FMIFlux   or   add "https://github.com/ThummeTo/FMIFlux.jl"
# (4) Install FMIZoo via             add FMIZoo    or   add "https://github.com/ThummeTo/FMIZoo.jl"
# (5) Install Flux via               add Flux
# (6) Install DifferentialEquations  add DifferentialEquations
# (7) Install Plots via              add Plots
# (8) Install Random via             add Random
################################ END INSTALLATION ##################################################

# this example covers creation and training of a ME-NeuralFMU
using FMI
using FMI.FMIImport: fmi2StringToValueReference, fmi2ValueReference, fmi2Real
using FMIFlux
using FMIZoo
using Flux
using DifferentialEquations: Tsit5
import Plots

# set seed
import Random
Random.seed!(1234);

tStart = 0.0
tStep = 0.1
tStop = 5.0
tSave = collect(tStart:tStep:tStop)

realFMU = fmiLoad("SpringFrictionPendulum1D", "Dymola", "2022x")
fmiInstantiate!(realFMU; loggingOn=false)
fmiSetupExperiment(realFMU, tStart, tStop)

fmiEnterInitializationMode(realFMU)
fmiExitInitializationMode(realFMU)

x₀ = fmiGetContinuousStates(realFMU)

vrs = ["mass.s", "mass.v", "mass.a", "mass.f"]
solution = fmiSimulate(realFMU, tStart, tStop; recordValues=vrs, saveat=tSave, reset=false)
fmiPlot(solution)

realSimData = solution.values.saveval
posReal = collect(data[1] for data in realSimData)
velReal = collect(data[2] for data in realSimData)
fmiUnload(realFMU)

simpleFMU = fmiLoad("SpringPendulum1D", "Dymola", "2022x")
fmiInstantiate!(simpleFMU; loggingOn=false)
solution = fmiSimulate(simpleFMU, tStart, tStop; recordValues=vrs[1:end-1], saveat=tSave, reset=false)
simpleSimData = solution.values.saveval
posSimple = collect(data[1] for data in simpleSimData)

# loss function for training
function lossSum()
    global neuralFMU

    solution = neuralFMU(x₀)

    posNet = collect(data[1] for data in solution.states.u)
    # velNet = collect(data[2] for data in solution.states.u)

    Flux.mse(posReal, posNet)
end

# callback function for training
global counter = 0
function callb()
    global counter += 1

    if counter % 30 == 1
        avgLoss = lossSum()
        @info "   Loss [$counter]: $(round(avgLoss, digits=5))
        Avg displacement in data: $(round(sqrt(avgLoss), digits=5))"

        fig = plotResults()
        println("Fig. update.")
        display(fig)
    end
end

function plotResults()
    solutionAfter = neuralFMU(x₀, tStart)
    fig = Plots.plot(xlabel="t [s]", ylabel="mass position [m]", linewidth=2,
        xtickfontsize=12, ytickfontsize=12,
        xguidefontsize=12, yguidefontsize=12,
        legendfontsize=12, legend=:topright)

    t = solutionAfter.states.t
    posNeuralFMU = collect(data[1] for data in solutionAfter.states.u)

    Plots.plot!(fig, tSave, posSimple, label="SimpleFMU", linewidth=2)
    Plots.plot!(fig, tSave, posReal, label="RealFMU", linewidth=2)
    Plots.plot!(fig, t, posNeuralFMU, label="NeuralFMU", linewidth=2)
    fig
end

# NeuralFMU setup
numStates = fmiGetNumberOfStates(simpleFMU)
additionalVRs = [fmi2StringToValueReference(simpleFMU, "mass.m")]
numAdditionalVRs = length(additionalVRs)

net = Chain(inputs -> fmiEvaluateME(simpleFMU, inputs, -1.0, zeros(fmi2ValueReference, 0), zeros(fmi2Real, 0), additionalVRs),
            Dense(numStates+numAdditionalVRs, 16, tanh),
            Dense(16, 16, tanh),
            Dense(16, numStates))

neuralFMU = ME_NeuralFMU(simpleFMU, net, (tStart, tStop), Tsit5(); saveat=tSave)
solutionBefore = neuralFMU(x₀, tStart)
fmiPlot(solutionBefore)

# train it ...
paramsNet = Flux.params(neuralFMU)

optim = ADAM()
# Feel free to increase training steps or epochs for better results
Flux.train!(lossSum, paramsNet, Iterators.repeated((), 1000), optim; cb=callb)

solutionAfter = neuralFMU(x₀, tStart)

fig = Plots.plot(xlabel="t [s]", ylabel="mass position [m]", linewidth=2,
                 xtickfontsize=12, ytickfontsize=12,
                 xguidefontsize=12, yguidefontsize=12,
                 legendfontsize=8, legend=:topright)

posNeuralFMU = collect(data[1] for data in solutionAfter.states.u)

Plots.plot!(fig, tSave, posSimple, label="SimpleFMU", linewidth=2)
Plots.plot!(fig, tSave, posReal, label="RealFMU", linewidth=2)
Plots.plot!(fig, tSave, posNeuralFMU, label="NeuralFMU (1000 epochs)", linewidth=2)
fig

fmiUnload(simpleFMU)


text:
"""
    single_reduce!(x, y)

Inplace reduction of y into x.
"""
@inline single_reduce!(x, y) = begin x .+= y; x end

for T in [judiVector, judiWeights, judiWavefield]
    @eval @inline single_reduce!(x::$T, y) = begin push!(x, y); x end
end

@inline single_reduce!(x::Ref, y::Ref) = begin x[] += y[]; x end

@inline function single_reduce!(x::Tuple, y::Tuple)
    nx = length(x)
    ny = length(y)
    (nx == ny) || throw(ArgumentError("Incompatible lengths ($nx, $ny)"))
    for (xi, yi) in zip(x, y)
        single_reduce!(xi, yi)
    end
    x
end

"""
    safe_gc()

Generic GC, compatible with different julia versions of it.
"""
safe_gc() = try Base.GC.gc(); catch; gc() end

"""
    local_reduce!(future, other)

Reduce `other` future into local `future`. This is perform remotely on the julia worker
`future.where`.
"""
function local_reduce!(my_future::_TFuture, other_future::_TFuture)
    y = fetch(other_future)
    x = fetch(my_future)
    single_reduce!(x, y)
    y = []
    nothing
end

"""
    reduce_level!(futures, nleaf)

Reduce one level of the reduction tree consisting of `nleaf` futures. Each `leaf`
reduces into itself `futures[i]`.
"""
function reduce_level!(futures::Vector{_TFuture}, nleaf::Int)
    nleaf == 0 && return
    @sync for i = 1:2:2*nleaf
        @async remotecall_wait(local_reduce!, futures[i].where, futures[i], futures[i+1])
    end
    deleteat!(futures, 2:2:2*nleaf)
end

"""
    reduce!(x)
binary-tree reduction of a Vector of Futures `x`

Modified from `DistributedOperations.jl` (Striped from custom types in it) and extended to multiple outputs
with different reduction functions.

Copyright (c) 2020: Chevron U.S.A Inc.

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
"""
function reduce!(futures::Vector{_TFuture})
    # Get length a next power of two for binary reduction
    M = length(futures)
    L = round(Int,log2(prevpow(2,M)))
    m = 2^L
    # remainder
    R = M - m

    # Reduce the remainder
    reduce_level!(futures, R)

    # Binary tree reduction
    for l = L:-1:1
        reduce_level!(futures, 2^(l-1))
    end
    fetch(futures[myid()])
end


text:
"""
    TypeY(var_1::TypeX)

Lorem ipsum dolor sit amet.
"""
struct TypeY

  var_1::TypeX

end


text:
using Test, MPI
using ProxAL

testdir = @__DIR__

@testset "Integration tests" begin
    include("blockmodel.jl")
end

# We can finalize here as now we launch external processes

# Testing using 1 process
@testset "Sequential tests" begin
    include("single.jl")
end

# Testing using 2 processes

@testset "Parallel tests" begin
    mpiexec() do cmd
        run(`$cmd -n 2 $(Base.julia_cmd()) --project=$testdir/.. $testdir/distributed.jl`)
    end
    @test true

    mpiexec() do cmd
        run(`$cmd -n 1 $(Base.julia_cmd()) --project=$testdir/.. $testdir/distributed.jl`)
    end
    @test true
end


text:
export create_operator, NoiseSampler, mix_data!

function create_operator(; nrec = 204, nsrc = 102, sigma = 2000.0f0, sim_src = false)

    n = (256, 256)
    v0 = ones(Float32, n) * 2.5f0
    v0 = (v0 .* range(1.0f0, 4.5f0 / maximum(v0), length = n[1]))'
    s = 1.0f0 ./ v0
    s[:, 1:10] .= 1 / 1.5f0
    m0 = s .^ 2.0f0
    d = (20.0f0, 12.5f0)
    o = (0.0f0, 0.0f0)

    # Setup info and model structure.
    model0 = Model(n, d, o, m0)

    # Setup receiver geometry.
    xrec = range(0.0f0, stop = model0.d[1] * (model0.n[1] - 1.0f0), length = nrec)
    yrec = 0.0f0
    zrec = range(2.0f0 * model0.d[2], stop = 2.0f0 * model0.d[2], length = nrec)

    # Receiver sampling and recording time.
    timeR = 2000.0f0
    dtR = 1.378f0

    # Setup receiver structure.
    recGeometry = Geometry(xrec, yrec, zrec; dt = dtR, t = timeR, nsrc = nsrc)

    # Setup source geometry (cell array with source locations for each shot).
    if sim_src == false
        xsrc = convertToCell(
            range(0.0f0, stop = model0.d[1] * (model0.n[1] - 1.0f0), length = nsrc),
        )
        ysrc = convertToCell(range(0.0f0, stop = 0.0f0, length = nsrc))
        zsrc = convertToCell(
            range(2.0f0 * model0.d[2], stop = 2.0f0 * model0.d[2], length = nsrc),
        )
    else
        xsrc = Array{Any}(undef, nsrc)
        ysrc = Array{Any}(undef, nsrc)
        zsrc = Array{Any}(undef, nsrc)
        for j = 1:nsrc
            xsrc[j] = collect(
                range(0.0f0, stop = model0.d[1] * (model0.n[1] - 1.0f0), length = nsrc),
            )
            ysrc[j] = [0.0f0]
            zsrc[j] = collect(
                range(2.0f0 * model0.d[2], stop = 2.0f0 * model0.d[2], length = nsrc),
            )
        end
        xsrc = convert(Array{Array{Float32,1},1}, xsrc)
        ysrc = convert(Array{Array{Float32,1},1}, ysrc)
        zsrc = convert(Array{Array{Float32,1},1}, zsrc)
    end

    # Source sampling and number of time steps.
    timeS = 2000.0f0
    dtS = 1.378f0

    # Setup source structure.
    srcGeometry = Geometry(xsrc, ysrc, zsrc; dt = dtS, t = timeS)

    # Setup wavelet.
    f0 = 0.03f0
    if sim_src == false
        wavelet = ricker_wavelet(timeS, dtS, f0)
        q = judiVector(srcGeometry, wavelet)
    else
        # Create array with different time shifts of the wavelet.
        src_enc = randn(Float32, nsrc, nsrc) ./ sqrt(nsrc)
        wavelet = zeros(Float32, srcGeometry.nt[1], nsrc)
        for j = 1:nsrc
            wavelet[:, j] = ricker_wavelet(timeS, dtS, f0)
        end

        q = judiVector(srcGeometry, wavelet)
        for j = 1:nsrc
            q.data[j] .= q.data[j] .* src_enc[j:j, :]
        end
    end

    # Setup options for forward modeling.
    ntComp = get_computational_nt(srcGeometry, recGeometry, model0, dt = dtS)
    info = Info(prod(n), nsrc, ntComp)
    opt = Options(
        optimal_checkpointing = false,
        isic = false,
        space_order = 16,
        free_surface = false,
        save_data_to_disk = false,
        return_array = true,
    )

    # Setup operators.
    Pr = judiProjection(info, recGeometry)
    F0 = judiModeling(info, model0; options = opt)
    Ps = judiProjection(info, srcGeometry)
    J = judiJacobian(Pr * F0 * adjoint(Ps), q)

    noise_dist =
        NoiseSampler(ricker_wavelet(timeS, dtS, f0)[1:50], sigma, Int(size(J)[1] / nsrc))

    if sim_src == false
        return J, noise_dist
    else
        return J, noise_dist, src_enc
    end
end

# Band-limited noise sampler.
mutable struct NoiseSampler
    op::Conv
    kernel::Array{Float32,1}
    sigma::Float32
    n::Int64
end

function NoiseSampler(kernel::Array{Float32,1}, sigma::Float32, n::Int64)
    op = Conv(reshape(kernel, :, 1, 1), zeros(Float32, 1), identity; pad = SamePad())
    return NoiseSampler(op, kernel, sigma, n)
end

function rand(NS::NoiseSampler)
    e = randn(Float32, (NS.n, 1, 1))
    e = NS.op(e)
    e *= 1.0f0 / norm(e)
    e *= sqrt(NS.n) * NS.sigma
    return e[:, 1, 1]
end

# Mix data according to the source encodings.
function mix_data!(d, src_enc)
    @info "Mixing shots to create simultaneous source data"
    @assert size(src_enc, 2) == size(d, 1)
    d_tmp = zeros(Float32, size(d))
    p = Progress(size(d, 1))
    for i = 1:size(d, 1)
        Base.flush(Base.stdout)
        for j = 1:size(d, 1)
            d_tmp[i, :] = d_tmp[i, :] + d[j, :] * src_enc[i, j]
        end
        ProgressMeter.next!(p; showvalues = [(:Shot, i)])
    end
    d .= d_tmp
    d_tmp = 0.0f0
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((5, 2, 2, true, true, true),),
    ((5, 1, 2, false, true, true),),
    ((5, 1, 2, false, true, false),),
    ((5, 1, 2, true, false, false), :ExpPSD),
    ((5, 1, 2, true, false, false), :SOCExpPSD),
    ((:iris, 2, true, false, true),),
    ]
insts["fast"] = [
    ((10, 1, 10, true, false, false), :ExpPSD),
    ((10, 1, 10, true, false, false), :SOCExpPSD),
    ((100, 1, 250, true, true, true),),
    ((100, 2, 5, false, true, true),),
    ((200, 2, 20, true, true, true),),
    ((50, 3, 2, true, true, true),),
    ((50, 3, 2, false, false, true),),
    ((50, 3, 2, true, false, true), :ExpPSD),
    ((50, 3, 2, true, false, true), :SOCExpPSD),
    ((50, 3, 4, true, true, true),),
    ((500, 3, 14, true, true, true),),
    ((100, 8, 2, true, true, true),),
    ((250, 4, 4, true, true, true),),
    ((250, 4, 4, false, false, true),),
    ((200, 32, 2, true, false, true),),
    ((:iris, 4, false, false, false),),
    ((:iris, 5, true, false, false),),
    ((:iris, 6, true, true, true),),
    ((:cancer, 4, true, true, true),),
    ]
insts["various"] = [
    ((50, 2, 16, true, false, false),),
    ((50, 2, 16, true, false, false), :ExpPSD),
    ((50, 2, 16, true, false, false), :SOCExpPSD),
    ((50, 2, 16, false, false, false),),
    ((50, 2, 16, false, false, true),),
    ((50, 8, 4, true, false, false),),
    ((50, 8, 4, true, false, false), :ExpPSD),
    ((50, 8, 4, true, false, false), :SOCExpPSD),
    ((50, 8, 4, false, false, false),),
    ((50, 8, 4, false, false, true),),
    ((50, 32, 2, true, false, false),),
    ((50, 32, 2, true, false, false), :ExpPSD),
    ((50, 32, 2, true, false, false), :SOCExpPSD),
    ((:iris, 6, true, true, true),),
    ((:cancer, 4, true, true, true),),
    ]
return (DensityEstJuMP, insts)


text:
using Documenter, Expectations

# Compile the raw documentation.
makedocs(sitename = "Expectations.jl")

# Push the documentation to the server.
deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/QuantEcon/Expectations.jl.git",
    forcepush = true)


text:
struct ArVar{Ti <: Integer}
    N::Int
    M::Int
    q::Int
    q2::Int
    lambdaJ::Float64
    lambdaH::Float64
    Z::Array{Ti,2}
    W::Array{Float64,1}
    IdxZ::Array{Int,2} # partial index computation to speed up energy calculation
    idxperm::Array{Int,1}

    function ArVar(N, M, q, lambdaJ, lambdaH, Z::Array{Ti,2}, W::Array{Float64,1}, permorder::Union{Symbol,Vector{Int}}) where Ti <: Integer
        idxperm = if typeof(permorder) == Symbol
            S = entropy(Z,W)
            if permorder === :ENTROPIC
                sortperm(S)
            elseif permorder === :REV_ENTROPIC
                sortperm(S,rev=true)
            elseif permorder === :RANDOM
                randperm(N)
            elseif permorder === :NATURAL
                collect(1:N)
            else
                error("the end of the world has come")
            end
        elseif typeof(permorder) <: Vector
            (length(permorder) != N) && error("length permorder ≠ $N")
            isperm(permorder) && (permorder)
        else
            error("permorder can only be a Symbol or a Vector")
        end
        permuterow!(Z,idxperm)
        IdxZ = Array{Int,2}(undef, N, M)
        q2 = q * q
        for i in 1:M
            for j in 1:N
                IdxZ[j,i] = (j - 1) * q2 + q * (Z[j,i] - 1)
            end
        end
        new{Ti}(N, M, q, q^2, lambdaJ, lambdaH, Z, W, IdxZ,idxperm)
    end
end

function Base.show(io::IO, arvar::ArVar)
    @extract arvar : N M q lambdaJ lambdaH Z
    print(io,"ArVar [N=$N M=$M q=$q λJ = $lambdaJ λH = $lambdaH Z::$(eltype(Z))]")
end
struct ArAlg
    method::Symbol
    verbose::Bool
    epsconv::Float64
    maxit::Int
end
struct ArNet
    idxperm::Array{Int,1}
    p0::Array{Float64,1}
    J::Array{Array{Float64,3},1}
    H::Array{Array{Float64,1},1}
end

ArNet(θ,var::ArVar) = ArNet(var.idxperm, unpack_params(θ, var)...)

function Base.show(io::IO, arnet::ArNet)
    N = length(arnet.idxperm)
    q = length(arnet.H[1])
    print(io,"ArNet [N=$N q=$q]")
end

function (A::ArNet)(x::Vector{T}) where T <: Integer
    @extract A:J H p0 idxperm
    backorder = sortperm(idxperm)
    N = length(x)
    length(H) == N - 1 || throw(DimensionMismatch("incompatible size between input and fields"))
    q = length(p0)
    permuterow!(x,idxperm)
    permuterow!(_outputarnet(x, J, H, p0, N, q),backorder)
end

function (A::ArNet)(x::Matrix{T}) where T <: Integer
    @extract A : J H p0 idxperm
    backorder = sortperm(idxperm)
    N, M = size(x)
    length(H) == N - 1 || throw(DimensionMismatch("incompatible size between input and fields"))
    q = length(p0)
    output = Array{eltype(p0),2}(undef,N,M)
    #Threads.@threads
    for i in 1:M
        output[:,i] .= _outputarnet(view(x,:,i), J, H, p0, N, q)
    end
    permuterow!(output,backorder)
end

(A::ArNet)(arvar::ArVar) = A(arvar.Z)

function _outputarnet( xs, J, H, p0, N, q)
    dest = Vector{Float64}(undef,N)
    _outputarnet!(dest, xs, J, H, p0, N, q)
end

let DtotH = Dict{Tuple{Int,Int},Vector{Float64}}()
    global _outputarnet!
    function _outputarnet!(dest, x, J, H, p0, N, q)
        dest[1] = p0[x[1]]
        totH = Base.get!(DtotH,(q,Threads.threadid()),Vector{Float64}(undef,q))
        #fill!(totH,0.0)
        @inbounds for site in 1:N-1
            Js = J[site]
            h = H[site]
            copy!(totH,h)
            @avx for i in 1:site
                for a in 1:q
                    totH[a] += Js[a,x[i],i]
                end
            end
            softmax!(totH)
            dest[site+1]=totH[x[site+1]]
        end
        dest
    end
end

max_issues_repo_name: pagnani/ArDCA

max_forks_repo_name: pagnani/ArDCA


text:
"""
    nearneighbor(cmd0::String="", arg1=nothing; kwargs...)

Reads arbitrarily located (x,y,z[,w]) triples [quadruplets] and uses a nearest neighbor algorithm
to assign an average value to each node that have one or more points within a radius centered on the node.
The average value is computed as a weighted mean of the nearest point from each sector inside the search
radius. The weighting function used is w(r) = 1 / (1 + d ^ 2), where d = 3 * r / search_radius and r is
distance from the node. This weight is modulated by the weights of the observation points [if supplied].

Full option list at [`nearneighbor`]($(GMTdoc)nearneighbor.html)

Parameters
----------

- **I** | **inc** :: [Type => Str | Number]

    *x_inc* [and optionally *y_inc*] is the grid spacing.
    ($(GMTdoc)nearneighbor.html#i)
- **N** | **sectors** | **nn** | **nearest** :: [Type => Number | Str | Bool (for nn or nearest)]

    The circular area centered on each node is divided into `sectors` sectors.
    ($(GMTdoc)nearneighbor.html#n)
- $(GMT.opt_R)
- **S** | **search_radius** :: [Type => Number]

    Sets the search_radius that determines which data points are considered close to a node.
    ($(GMTdoc)nearneighbor.html#s)

- **E** | **empty** :: [Type => Bool]

    Set the value assigned to empty nodes when G is set [NaN].
    ($(GMTdoc)nearneighbor.html#e)
- **G** | **outgrid** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = nearneighbor(....) form.
    ($(GMTdoc)nearneighbor.html#g)
- $(GMT.opt_V)
- **W** | **weights** :: [Type => Bool]

    Input data have a 4th column containing observation point weights.
    ($(GMTdoc)nearneighbor.html#w)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_n)
- $(GMT.opt_r)
- $(GMT.opt_swap_xy)
"""
function nearneighbor(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("nearneighbor", cmd0, arg1)

	d = KW(kwargs)
	help_show_options(d)		# Check if user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:R :I :V_params :bi :di :e :f :h :i :n :r :yx])
	cmd  = parse_these_opts(cmd, d, [[:E :empty], [:G :outgrid], [:S :search_radius], [:Z :weights], [:A]])
	cmd  = add_opt(d, cmd, 'N', [:N :sectors], (n="", min_sectors="+m"), true)
	opt  = add_opt(d, "", 'N', [:N :nn :nearest])
	if (opt != "")  cmd *= " -Nn"  end

	if (isa(arg1, Matrix{<:Number}))  arg1 = GMTdataset(arg1)  end      # Must find why need this
	common_grd(d, cmd0, cmd, "nearneighbor ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
nearneighbor(arg1, cmd0::String=""; kw...) = nearneighbor(cmd0, arg1; kw...)


text:
struct Wrapper
  tickPrice::Function
  tickSize::Function
  tickOptionComputation::Function
  tickGeneric::Function
  tickString::Function
  tickEFP::Function
  orderStatus::Function
  openOrder::Function
  openOrderEnd::Function
#  winError::Function
#  connectionClosed::Function
  updateAccountValue::Function
  updatePortfolio::Function
  updateAccountTime::Function
  accountDownloadEnd::Function
  nextValidId::Function
  contractDetails::Function
  bondContractDetails::Function
  contractDetailsEnd::Function
  execDetails::Function
  execDetailsEnd::Function
  error::Function
  updateMktDepth::Function
  updateMktDepthL2::Function
  updateNewsBulletin::Function
  managedAccounts::Function
  receiveFA::Function
  historicalData::Function
  scannerParameters::Function
  scannerData::Function
  realtimeBar::Function
  currentTime::Function
  fundamentalData::Function
  deltaNeutralValidation::Function
  tickSnapshotEnd::Function
  marketDataType::Function
  commissionReport::Function
  position::Function
  positionEnd::Function
  accountSummary::Function
  accountSummaryEnd::Function
  verifyMessageAPI::Function
  verifyCompleted::Function
  displayGroupList::Function
  displayGroupUpdated::Function
  verifyAndAuthMessageAPI::Function
  verifyAndAuthCompleted::Function
#  connectAck::Function
  positionMulti::Function
  positionMultiEnd::Function
  accountUpdateMulti::Function
  accountUpdateMultiEnd::Function
  securityDefinitionOptionalParameter::Function
  securityDefinitionOptionalParameterEnd::Function
  softDollarTiers::Function
  familyCodes::Function
  symbolSamples::Function
  mktDepthExchanges::Function
  tickNews::Function
  smartComponents::Function
  tickReqParams::Function
  newsProviders::Function
  newsArticle::Function
  historicalNews::Function
  historicalNewsEnd::Function
  headTimestamp::Function
  histogramData::Function
  historicalDataUpdate::Function
  rerouteMktDataReq::Function
  rerouteMktDepthReq::Function
  marketRule::Function
  pnl::Function
  pnlSingle::Function
  historicalTicks::Function
  historicalTicksBidAsk::Function
  historicalTicksLast::Function
  tickByTickAllLast::Function
  tickByTickBidAsk::Function
  tickByTickMidPoint::Function
  orderBound::Function
  completedOrder::Function
  completedOrdersEnd::Function
end
function Wrapper(; kw...)

  default = (_...) -> println("Default Implementation")

  args = (get(kw, a, default) for a ∈ fieldnames(Wrapper))

  Wrapper(args...)
end

function simple_wrap()

  d = Dict{Symbol,Any}()

  w = Wrapper(

    tickPrice= (tickerId::Int, field::String, price::Float64, size::Int, attrib::TickAttrib) ->
                 println("Price: $tickerId $field $price $size $attrib"),

    tickSize= (tickerId::Int, field::String, size::Int) -> println("Size: $tickerId $field $size"),

    tickOptionComputation= function(tickerId::Int, tickType::String, impliedVol::Union{Float64,Nothing}, delta::Union{Float64,Nothing}, optPrice::Union{Float64,Nothing}, pvDividend::Union{Float64,Nothing}, gamma::Union{Float64,Nothing}, vega::Union{Float64,Nothing}, theta::Union{Float64,Nothing}, undPrice::Union{Float64,Nothing})
                             d[:option] = (tickType, impliedVol, delta, optPrice, pvDividend, gamma, vega, theta, undPrice)
                             println("Option: $tickerId $tickType")
                           end,

    tickGeneric= (tickerId::Int, tickType::String, value::Float64) ->
                   println("Generic: $tickerId $tickType $value"),

    tickString= (tickerId::Int, tickType::String, value::String) ->
                  println("String: $tickerId $tickType $value"),

    orderStatus= function(orderId::Int, status::String, filled::Float64, remaining::Float64, avgFillPrice::Float64, permId::Int, parentId::Int, lastFillPrice::Float64, clientId::Int, whyHeld::String, mktCapPrice::Union{Float64,Nothing})
                   d[:orderstatus] = (orderId=orderId,
                                      status=status,
                                      filled=filled,
                                      remaining=remaining,
                                      avgFillPrice=avgFillPrice,
                                      permId=permId,
                                      lastFillPrice=lastFillPrice,
                                      clientId=clientId,
                                      whyHeld=whyHeld,
                                      mktCapPrice=mktCapPrice)
                   println("OrderStatus: $orderId $status")
                 end,

    openOrder= function(orderId::Int, contract::Contract, order::Order, orderstate::OrderState)
                 d[:order_contract] = contract
                 d[:order] = order
                 d[:orderstate] = orderstate
                 println("OpenOrder: $orderId $(orderstate.status)")
              end,

    openOrderEnd= () -> println("OpenOrderEnd."),

    updateAccountValue= (key::String, val::String, currency::String, accountName::String) ->
                          println("AccountValue: $key $val $currency $accountName"),

    updatePortfolio= (contract::Contract, position::Float64, marketPrice::Float64, marketValue::Float64, averageCost::Float64, unrealizedPNL::Float64, realizedPNL::Float64, accountName::String) ->
                       println("Portfolio:  $(contract.symbol) $position $marketPrice $marketValue $averageCost $unrealizedPNL $realizedPNL $accountName"),

    updateAccountTime= (timeStamp::String) -> println("AccountTime: $timeStamp"),

    accountDownloadEnd= (accountName::String) -> println("AccountDownloadEnd: $accountName"),

    nextValidId= (orderId::Int) -> d[:nextId] = orderId,

    contractDetails = function(reqId::Int, contractDetails::ContractDetails)
                        d[:cd] = contractDetails
                        println("ContractDetails: $reqId ", contractDetails.contract.conId)
                      end,

    contractDetailsEnd= (reqId::Int) -> println("ContractDetailsEnd: $reqId"),

    execDetails= function(reqId::Int, contract::Contract, execution::Execution)
                    d[:ex_con] = contract
                    d[:execution] = execution
                    println("ExecDetails: $reqId")
                 end,

    execDetailsEnd= (reqId::Int) -> println("ExecDetailsEnd: $reqId"),

    error= (id::Union{Int,Nothing}, errorCode::Int, errorString::String) ->
                     println("Error: $(something(id, "NA")) $errorCode $errorString"),

    updateMktDepth= (id::Int, position::Int, operation::Int, side::Int, price::Float64, size::Int) ->
                      println("MktDepth: $id $position $operation $side $price $size"),

    updateMktDepthL2= (id::Int, position::Int, marketMaker::String, operation::Int, side::Int, price::Float64, size::Int, isSmartDepth::Bool) ->
                        println("MktDepthL2: $id $position $marketMaker $operation $side $price $size $isSmartDepth"),

    updateNewsBulletin= function(msgId::Int, msgType::Int, newsMessage::String, originExch::String)
                          d[:newsbulletin] = newsMessage
                          println("NewsBulletin: $msgId $msgType $originExch")
                        end,

    managedAccounts= (accountsList::String) ->d[:accounts] = accountsList,

    receiveFA= function(pFaDataType::faDataType, cxml::String)
                 d[:fa] = cxml
                 println("ReceiveFA: $pFaDataType")
               end,

    historicalData= function(reqId::Int, bar::DataFrame)
                      d[:history] = bar
                      println("HistoricalData: $reqId $(size(bar))")
                    end,

    scannerParameters= function(xml::String)
                        d[:scannerparam] = xml
                        println("ScannerParameters.")
                       end,

    scannerData= function(reqId::Int, rank::Vector{Int}, contractDetails::Vector{ContractDetails}, distance::Vector{String}, benchmark::Vector{String}, projection::Vector{String}, legsStr::Vector{String})
                   d[:scannerdata] = (rank=rank,
                                      cd=contractDetails,
                                      distance=distance,
                                      benchmark=benchmark,
                                      projection=projection,
                                      legsStr=legsStr)
                   println("ScannerData: $reqId")
                 end,

    realtimeBar= (reqId::Int, time::Int, open::Float64, high::Float64, low::Float64, close::Float64, volume::Int, wap::Float64, count::Int) ->
                   println("RealTimeBar: $reqId $time $open $high, $low $close $volume $wap $count"),

    currentTime= (time::Int) -> println("CurrentTime: $time"),

    fundamentalData= function(reqId::Int, data::String)
                       d[:fundamental] = data
                       println("FundamentalData: $reqId")
                     end,

    tickSnapshotEnd= (reqId::Int) -> println("TickSnapshotEnd: $reqId"),

    marketDataType= (reqId::Int, marketDataType::MarketDataType) -> println("MarketDataType: $reqId $marketDataType"),

    commissionReport= function(commissionReport::CommissionReport)
                        d[:commission] = commissionReport
                        println("CommissionReport.")
                      end,

    position= (account::String, contract::Contract, position::Float64, avgCost::Float64) ->
                println("Position: $account $(contract.symbol) $position $avgCost"),

    positionEnd= () -> println("PositionEnd."),

    accountSummary= (reqId::Int, account::String, tag::String, value::String, currency::String) -> println("AccountSummary: $reqId $account $tag $value $currency"),

    accountSummaryEnd= (reqId::Int) -> println("AccountSummaryEnd: $reqId"),

    displayGroupList= (reqId::Int, groups::String) -> println("DisplayGroupList: $reqId $groups"),

    displayGroupUpdated= (reqId::Int, contractInfo::String) -> println("DisplayGroupUpdated: $reqId $contractInfo"),

    positionMulti= (reqId::Int, account::String, modelCode::String, contract::Contract, position::Float64, avgCost::Float64) ->
                     println("PositionMulti: $reqId $account $modelCode $(contract.symbol) $position $avgCost"),

    positionMultiEnd= (reqId::Int) -> println("PositionMultiEnd $reqId"),

    accountUpdateMulti= (reqId::Int, account::String, modelCode::String, key::String, value::String, currency::String) -> println("AccountUpdateMulti: $reqId $account, $modelCode $key=$value $currency"),

    accountUpdateMultiEnd= (reqId::Int) -> println("AccountUpdateMultiEnd: $reqId"),

    securityDefinitionOptionalParameter= function(reqId::Int, exchange::String, underlyingConId::Int, tradingClass::String, multiplier::String, expirations::Vector{String}, strikes::Vector{Float64})
                                           d[:sdop] = (e=expirations, s=strikes)
                                           println("SDOP: $reqId $exchange $underlyingConId $tradingClass $multiplier $(length(expirations)) $(length(strikes))")
                                         end,

    securityDefinitionOptionalParameterEnd= (reqId::Int) -> println("SDOPEnd: $reqId"),

    softDollarTiers= function(reqId::Int, tiers::Vector{SoftDollarTier})
                       d[:softdollar] = tiers
                       println("SoftDollarTiers: $reqId $(length(tiers))")
                     end,

    familyCodes= function(familyCodes::Vector{FamilyCode})
                   d[:familycodes] = familyCodes
                   println("FamilyCodes: $(length(familyCodes))")
                 end,

    symbolSamples= function(reqId::Int, contractDescriptions::Vector{ContractDescription})
                     d[:symbolsamples] = contractDescriptions
                     println("SymbolSamples: $reqId")
                   end,

    mktDepthExchanges= function(depthMktDataDescriptions::DataFrame)
                        d[:mktdepthexchanges] = depthMktDataDescriptions
                        println("MktDepthExchanges.")
                      end,

    tickNews= (tickerId::Int, timeStamp::Int, providerCode::String, articleId::String, headline::String, extraData::String) ->
                println("TickNews: $tickerId $timeStamp $providerCode $articleId $headline"),

    smartComponents= function(reqId::Int, theMap::DataFrame)
                       d[:smartcomponents] = theMap
                       println("SmartComponents: $reqId $(size(theMap))")
                     end,

    tickReqParams= (tickerId::Int, minTick::Union{Float64,Nothing}, bboExchange::String, snapshotPermissions::Int) ->
                     println("TickReqParams: $tickerId ",
                             something(minTick, "NA"),
                             " $bboExchange $snapshotPermissions"),

    newsProviders= function(newsProviders::DataFrame)
                     d[:newsproviders] = newsProviders
                     println("NewsProviders.")
                   end,

    newsArticle= function(requestId::Int, articleType::Int, articleText::String)
                   d[:newsarticle] = articleText
                   println("NewsArticle: $requestId $articleType")
                 end,

    headTimestamp= (reqId::Int, headTimestamp::String) -> println("HeadTimestamp: $reqId $headTimestamp"),

    histogramData= function(reqId::Int, data::DataFrame)
                    d[:histogram] = data
                    println("HistogramData: $reqId")
                  end,

    historicalDataUpdate= (reqId::Int, bar::NamedTuple) ->
                            println("HistoricalDataUpdate: $reqId $bar"),

    marketRule= function(marketRuleId::Int, priceIncrements::DataFrame)
                  d[:marketrule] = priceIncrements
                  println("MarketRule: $marketRuleId")
                end,

    pnl= (reqId::Int, dailyPnL::Float64, unrealizedPnL::Union{Float64,Nothing}, realizedPnL::Union{Float64,Nothing}) ->
           println("PnL: $reqId $dailyPnL ",
                   something(unrealizedPnL, "NA"), " ",
                   something(realizedPnL, "NA")),

    pnlSingle= (reqId::Int, pos::Int, dailyPnL::Float64, unrealizedPnL::Union{Float64,Nothing}, realizedPnL::Union{Float64,Nothing}, value::Float64) ->
                 println("PnLSingle: $reqId $pos $dailyPnL ",
                         something(unrealizedPnL, "NA"), " ",
                         something(realizedPnL, "NA"),   " ",
                         value),

    orderBound= (orderId::Int, apiClientId::Int, apiOrderId::Int) ->
                  println("OrderBound $orderId $apiClientId $apiOrderId"),

    completedOrder= function(contract::Contract, order::Order, orderState::OrderState)
                      d[:completed_contract] = contract
                      d[:completed] = order
                      d[:completed_state] = orderState
                 println("CompletedOrder: $(contract.symbol) $(orderState.status)")
              end,

    completedOrdersEnd= () -> println("CompletedOrdersEnd."),
  )

  d, w
end


text:
@inline function check_bind_bounds(stmt::Stmt, pos::Integer)
    @assert pos > 0 && pos <= stmt.bind_count "Bind position $pos out of bounds."
end

@inline function check_bind_bounds(stmt::Stmt, name::AbstractString)
    name_upper = uppercase(name)
    @assert haskey(stmt.bind_names_index, name_upper) "Bind name $name_upper not found in statement."
end

@inline check_bind_bounds(stmt::Stmt, name::Symbol) = check_bind_bounds(stmt, string(name))

const NameOrPositionTypes = Union{Integer, String, Symbol}

#
# Bind Variable to Stmt
#

@inline function Base.setindex!(stmt::Stmt, value::Variable, name_or_position::NameOrPositionTypes)
    bind!(stmt, value, name_or_position)
end

@generated function bind!(stmt::Stmt, variable::Variable, name_or_position::K) where {K<:NameOrPositionTypes}

    if name_or_position <: Integer
        name_or_position_exp = :(UInt32(name_or_position))
        bind_function_name = :dpiStmt_bindByPos
    elseif name_or_position <: Symbol
        name_or_position_exp = :(string(name_or_position))
        bind_function_name = :dpiStmt_bindByName
    elseif name_or_position <: String
        name_or_position_exp = :name_or_position
        bind_function_name = :dpiStmt_bindByName
    else
        error("Unsupported type for argument `name_or_position`: $name_or_position.")
    end

    return quote
        key = $name_or_position_exp
        check_bind_bounds(stmt, key)
        result = $(bind_function_name)(stmt.handle, key, variable.handle)
        error_check(context(stmt), result)
        nothing
    end
end

#
# Bind Value to Stmt
#

@generated function bind!(stmt::Stmt, value::JuliaOracleValue{O,N}, name_or_position::K) where {K<:NameOrPositionTypes,O,N}

    # https://github.com/oracle/odpi/issues/99
    if O == ORA_ORACLE_TYPE_TIMESTAMP_TZ || O == ORA_ORACLE_TYPE_TIMESTAMP_LTZ
        error("Can't bind Timestamp with TimeZone directly to Statement. Use a Variable instead.")
    end

    if O == ORA_ORACLE_TYPE_RAW || O == ORA_ORACLE_TYPE_LONG_RAW
        error("Can't bind RAW data type directly to Statement. Use a Variable instead.")
    end

    if name_or_position <: Integer
        name_or_position_exp = :(UInt32(name_or_position))
        bind_function_name = :dpiStmt_bindValueByPos
    elseif name_or_position <: Symbol
        name_or_position_exp = :(string(name_or_position))
        bind_function_name = :dpiStmt_bindValueByName
    elseif name_or_position <: String
        name_or_position_exp = :name_or_position
        bind_function_name = :dpiStmt_bindValueByName
    else
        error("Unsupported type for argument `name_or_position`: $name_or_position.")
    end

    return quote
        key = $name_or_position_exp
        check_bind_bounds(stmt, key)
        data_ref = Ref{OraData}()
        set_oracle_value_at!(value, value[], data_ref)
        result = $(bind_function_name)(stmt.handle, key, N, data_ref)
        error_check(context(stmt), result)
        nothing
    end
end

@inline function Base.setindex!(stmt::Stmt, value::JuliaOracleValue, name_or_position::N) where {N<:NameOrPositionTypes}
    bind!(stmt, value, name_or_position)
end

@inline function Base.setindex!(stmt::Stmt, value::T, name_or_position::N) where {T, N<:NameOrPositionTypes}
    bind!(stmt, JuliaOracleValue(value), name_or_position)
end

@inline function Base.setindex!(stmt::Stmt, ::Missing, name_or_position::N) where {N<:NameOrPositionTypes}
    error("Cannot bind missing value to statement without type information. Use `stmt[pos, julia_type] = value` or `stmt[pos, oracle_type, native_type] = value`.")
end

@inline function Base.setindex!(stmt::Stmt, ::Missing, name_or_position::N, oracle_type::OraOracleTypeNum, native_type::OraNativeTypeNum) where {N<:NameOrPositionTypes}
    val = JuliaOracleValue(oracle_type, native_type, Missing)
    val[] = missing
    bind!(stmt, val, name_or_position)
end

@inline function Base.setindex!(stmt::Stmt, m::Missing, name_or_position::N, ::Type{T}) where {T,N<:NameOrPositionTypes}
    ott = infer_oracle_type_tuple(T)
    setindex!(stmt, m, name_or_position, ott.oracle_type, ott.native_type)
end


text:
using Turf

shape = GeoJSON.parsefile(joinpath(dirname(pathof(Turf)), "..", "docs") * "/examples/simplify/geometry.geojson")

simplify!(shape, 0.4)

result = GeoJSON.geojson(shape)

open(joinpath(dirname(pathof(Turf)), "..", "docs") * "/examples/simplify/geometry.result.geojson", "w") do file
    write(file, result)
end


text:
# power analysis functions for mean for normal distribution with known variance

"""
    powerRTest(r::Real = 0, n::Real = 0, alpha::Float64 = 0.05, alternative = "two.sided")

Compute power of a sample of `n` observations to test the effect size `r`
at type I error = `alpha` (default: 0.05).
The effect size `r` is a linear correlation coefficient. The option
`alternative` is `two.sided` for H₀: r = 0 vs H₁: r ≠ 0, `less` for H₁: r < 0, and `greater`
for H₁: r > 0.
"""
function powerRTest(;
    r::Real = 0.0,
    n::Real = 0,
    alpha::Real = 0.05,
    alternative::String = "two")

    check_args(r=r,n=n,alpha=alpha)

    if alternative == "less"
        tside = 1
        r = -r
    elseif alternative in ("two","two.sided","two-sided","two sided")
        tside = 2
        r = abs(r)
    elseif alternative == "greater"
        tside = 1
    end

    if tside == 1
        ttt = cquantile(TDist(n-2),alpha)
        rc = sqrt(ttt^2 / (ttt^2 + n - 2))
        zr = atanh(r) + r / (2 * (n - 1))
        zrc = atanh(rc)
        return cdf(Normal(),(zr - zrc) * sqrt(n-3))
    end
    ttt = cquantile(TDist(n-2),alpha/2)
    rc = sqrt(ttt^2 / (ttt^2 + n - 2))
    zr = atanh(r) + r / (2 * (n - 1))
    zrc = atanh(rc)
    return cdf(Normal(),(zr - zrc) * sqrt(n-3)) + cdf(Normal(),(-zr - zrc) * sqrt(n-3))
end

"""
    samplesizeRTest(r::Real = 0, alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute the sample size to detect the effect size `r` with power > `power`
at type I error = `alpha` (default: 0.05).
The effect size `r` is a linear correlation coefficient. The option
`alternative` is `two.sided` for H₀: r = 0 vs H₁: r ≠ 0, `less` for H₁: r < 0, and `greater`
for H₁: r > 0.
"""
function samplesizeRTest(;
    r::Real = 0.0,
    alpha::Real = 0.05,
    power::Real = 0.0,
    alternative::String = "two")

    check_args(r=r,alpha=alpha,power=power)

    return ceil(Int64,fzero(x->powerRTest(n = x, r = r, alpha = alpha, alternative = alternative) - power, 4.0+1e-10, 1e+09))
end

"""
    effectsizeRTest(n::Real = 0, alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute the effect size that a sample size `n` can detect with power > `power`
at type I error = `alpha` (default: 0.05).
The effect size `r` is a linear correlation coefficient. The option
`alternative` is `two.sided` for H₀: r = 0 vs H₁: r ≠ 0, `less` for H₁: r < 0, and `greater`
for H₁: r > 0.
"""
function effectsizeRTest(;
    n::Real = 0,
    alpha::Real = 0.05,
    power::Real = 0.0,
    alternative::String = "two"
    )

    check_args(n=n,alpha=alpha,power=power)

    return fzero(x -> powerRTest(n = n, r = x, alpha = alpha, alternative = alternative) - power,.001,100)
end

"""
    alphaRTest(r::Real = 0, n::Real = 0, power::Float64 = 0.8, alternative = "two.sided")

Compute the probability of type I error that the sample size `n` has in detecting the
effect size `r` with power > `power` at type I error = `alpha` (default: 0.05).
The effect size `r` is a linear correlation coefficient. The option
`alternative` is `two.sided` for H₀: r = 0 vs H₁: r ≠ 0, `less` for H₁: r < 0, and `greater`
for H₁: r > 0.
"""
function alphaRTest(;
    n::Real = 0,
    r::Real = 0.0,
    power::Real = 0.0,
    alternative::String = "two"
    )

    check_args(r=r,n=n,power=power)

    return fzero(x->powerRTest(n = n, r = r, alpha = x, alternative = alternative) - power, 1e-10, 1 - 1e-10)
end

"""
    pwr.RTest(r::Real = 0, n::Real = 0, alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute one of the test parameters such as sample size (`n`),
effect size (`r`), type I error (`alpha`), and power (`power`).
The parameter to be estimated must be set to zero (default).
The effect size `r` is a linear correlation coefficient. The option
`alternative` is `two.sided` for H₀: r = 0 vs H₁: r ≠ 0, `less` for H₁: r < 0, and `greater`
for H₁: r > 0.
"""
function RTest(;
    n::Real = 0,
    r::Real = 0.0,
    alpha::Float64 = 0.0,
    power::Float64 = 0.0,
    alternative::String = "two")

    if sum([x == 0 for x in (n,r,alpha,power)]) != 1
        error("exactly one of n, r, power, and alpha must be zero")
    end

    if power == 0.0
        power = powerRTest(n = n, r = r, alpha = alpha, alternative = alternative)
    elseif alpha == 0.0
        alpha = alphaRTest(n = n, r = r, power = power, alternative = alternative)
    elseif r == 0.0
        r = effectsizeRTest(n = n, alpha = alpha, power = power, alternative = alternative)
    elseif n == 0
        n = samplesizeRTest(r = r, alpha = alpha, power = power, alternative = alternative)
    end

    alt = Dict("two" => "two-sided","two.sided" => "two-sided", "two sided" => "two-sided", "less" => "less", "greater" => "greater")

    return htest(
        string("Approximate correlation power calculation (arctangh transformation)"),
        OrderedDict(
            "n" => n,
            "r" => r,
            "alpha" => alpha,
            "power" => power,
            "alternative" => alt[alternative])
        )
end

max_stars_repo_name: mwsohn/pwr

max_issues_repo_name: mwsohn/pwr

max_forks_repo_name: mwsohn/pwr


text:
# ------------------------------------------------------------------------------------------
# # Julia for Data Science - Plotting
#
# ### Data visualization: generating nice looking plots in Julia is straight forward
# In what's next, we will see some of the tools that Julia plotting provides to produce high
# quality figures for your data. In particular we'll look at
#
# 1. Plotting mathematical functions
# 1. Visualizing statistics
# 1. Subplotting
#
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Part 1: plot math functions (specifically latex equations) in our plots
# ------------------------------------------------------------------------------------------

using LaTeXStrings
using Plots
pyplot(leg=false)
x = 1:0.2:4

# ------------------------------------------------------------------------------------------
# Create three functions and plot them all!
# ------------------------------------------------------------------------------------------

y1 = sqrt.(x)
y2 = log.(x)
y3 = x.^2

f1 = plot(x,y1)
plot!(f1,x,y2) # "plot!" means "plot on the same canvas we just plot on"
plot!(f1,x,y3)

# ------------------------------------------------------------------------------------------
# Now we can annotate each of these plots! using either text, or latex strings
# ------------------------------------------------------------------------------------------

annotate!(f1,[(x[6],y1[6],text(L"\sqrt{x}",16,:center)),
          (x[11],y2[11],text(L"log(x)",:right,16)),
          (x[6],y3[6],text(L"x^2",16))])

# ------------------------------------------------------------------------------------------
# ## Part 2: Visualizing statistics
#
# **2D histograms** are really easy!
# ------------------------------------------------------------------------------------------

n = 1000
set1 = randn(n)
set2 = randn(n)
histogram2d(set1,set2,nbins=20,colorbar=true)

# ------------------------------------------------------------------------------------------
# **Let's go back to our houses dataset and learn even more things about it!**
# ------------------------------------------------------------------------------------------

using DataFrames
houses = readtable("houses.csv")
filter_houses = houses[houses[!, :sq__ft].>0,:]
x = filter_houses[!, :sq__ft]
y = filter_houses[!, :price]

gh = histogram2d(x,y,nbins=20,colorbar=true)
xaxis!(gh,"square feet")
yaxis!(gh,"price")

# ------------------------------------------------------------------------------------------
# Interesting!
#
# Most houses sold are in the range 1000-1500 and they cost approximately 150,000 dollars
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# *Let's see more stats plots.*
#
# We can convince ourselves that random distrubutions are indeed very similar.
#
# Let's do that through **box plots** and **violin plots**.
# ------------------------------------------------------------------------------------------

using StatsPlots
y = rand(10000,6) # generate 6 random samples of size 1000 each
f2 = violin(["Series 1" "Series 2" "Series 3" "Series 4" "Series 5"],y,leg=false,color=:red)

boxplot!(["Series 1" "Series 2" "Series 3" "Series 4" "Series 5"],y,leg=false,color=:green)

# ------------------------------------------------------------------------------------------
# These plots look almost identical, so we do have the same distribution indeed.
#
# Let's study the price distributions in different cities in the houses dataset.
# ------------------------------------------------------------------------------------------

some_cities = ["SACRAMENTO","RANCHO CORDOVA","RIO LINDA","CITRUS HEIGHTS","NORTH HIGHLANDS","ANTELOPE","ELK GROVE","ELVERTA" ] # try picking pther cities!

fh = plot(xrotation=90)
for ucity in some_cities
    subs = filter_houses[filter_houses[!, :city].==ucity,:]
    city_prices = subs[!, :price]
    violin!(fh,[ucity],city_prices,leg=false)
end
display(fh)

# ------------------------------------------------------------------------------------------
# ## Part 3: Subplots are very easy!
#
# To create a plot with subplots, all we need to do is throw the variables bound to
# individual plots inside another call to `plot`!
# ------------------------------------------------------------------------------------------

x = -10:.1:10
p1 = plot(x, x)
p2 = plot(x, x.^2)
p3 = plot(x, x.^3)
p4 = plot(x, x.^4)
plot(p1,p2,p3,p4,layout=(2,2),legend=false)

# ------------------------------------------------------------------------------------------
# You can create your own layout as follows.
# ------------------------------------------------------------------------------------------

mylayout = @layout([a{0.5h};[b{0.7w} c]])
plot(fh,f2,gh,layout=mylayout,legend=false)

# this layout:
# 1
# 2 3

# ------------------------------------------------------------------------------------------
# ### Please let us know how we're doing!
#
# https://tinyurl.com/JuliaDataScience
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# # Wiring diagrams in TikZ
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/graphics/tikz_wiring_diagrams.ipynb)
#
# Catlab can draw morphism expressions as TikZ pictures. To use this feature,
# you must import the [TikzPictures.jl](https://github.com/sisl/TikzPictures.jl)
# package before loading `Catlab.Graphics`.

import TikzPictures
using Catlab.Graphics

# ### Symmetric monoidal category

using Catlab.Doctrines

A, B = Ob(FreeSymmetricMonoidalCategory, :A, :B)
f = Hom(:f, A, B)
g = Hom(:g, B, A)
h = Hom(:h, otimes(A,B), otimes(A,B));
#-
to_tikz(f; arrowtip="Stealth", labels=true)
#-
to_tikz(compose(f,g); arrowtip="Stealth", labels=true)
#-
comp1 = compose(otimes(g,f), h, otimes(f,g))
#-
to_tikz(comp1)
#-
comp2 = compose(otimes(id(A),f), h, otimes(id(A),g))
#-
to_tikz(comp2)
#-
to_tikz(Hom(:H, otimes(A,A), otimes(B,B,B)))
#-
twist = compose(braid(A,A), otimes(f,f), braid(B,B))
#-
to_tikz(twist)

# ### Biproduct category

A, B = Ob(FreeBiproductCategory, :A, :B)
f = Hom(:f, A, B)
g = Hom(:g, B, A);
#-
split = compose(mcopy(A), otimes(f,f))
#-
to_tikz(split)
#-
combine = compose(split, mmerge(B))
#-
to_tikz(combine, labels=true)
#-
to_tikz(compose(create(A), f, g, delete(A)))

# ### Compact closed category

A, B = Ob(FreeCompactClosedCategory, :A, :B)
f = Hom(:f, A, B)
g = Hom(:g, B, A)
h = Hom(:h, otimes(A,B), otimes(A,B));
#-
to_tikz(dcounit(A); arrowtip="Stealth", labels=true)
#-
to_tikz(dunit(A); arrowtip="Stealth", labels=true)
#-
transpose = compose(
  otimes(dunit(A), id(dual(B))),
  otimes(id(dual(A)), f, id(dual(B))),
  otimes(id(dual(A)), dcounit(B)),
)
#-
to_tikz(transpose)
#-
trace = compose(
  f,
  otimes(dunit(A), id(B)),
  otimes(id(dual(A)), h),
  otimes(braid(dual(A),A), id(B)),
  otimes(dcounit(A), id(B)),
  g
)
#-
to_tikz(trace)


text:
# ----------------- logit -----------------

# pdf
function logisticpdf(z::Real)
    isfinite(z) || return zero(z)
    lz = logistic(-abs(z))
    return lz * (1 - lz)
end

function logisticcdf(z::Real)
    isfinite(z)     && return logistic(z)
    z == typemin(z) && return zero(z)
    z == typemax(z) && return one(z)
    throw(DomainError())
end

# cdf
cdf(   z::Real, ::Type{Val{:logit}}) = logisticcdf(z)
ccdf(  z::Real, ::Type{Val{:logit}}) = logisticcdf(-z)
invcdf(z::Real, ::Type{Val{:logit}}) = logit(z)

@inline pdf(z::Real, ::Type{Val{:logit}}) = logisticpdf(z)

function dpdf(z::Real, ::Type{Val{:logit}})
    isfinite(z) || return zero(z)
    F = logistic(z)
    cF = one(z) - F
    return F * cF * (one(z)-2*F)
end

# logcdf
logcdf(   z::Real, ::Type{Val{:logit}}) = -log1pexp(-z)  # z - log1pexp(z)
logccdf(  z::Real, ::Type{Val{:logit}}) = -log1pexp( z)  #   - log1pexp(z)
dlogcdf(  z::Real, ::Type{Val{:logit}}) = one(z) - logistic(z)
dlogccdf( z::Real, ::Type{Val{:logit}}) =        - logistic(z)
d2logcdf( z::Real, ::Type{Val{:logit}}) = - pdf(z, Val{:logit})
d2logccdf(z::Real, ::Type{Val{:logit}}) = - pdf(z, Val{:logit})

dlogcdf_trunc(a::Real, b::Real, ::Type{Val{:logit}}) = (logisticpdf(b) - logisticpdf(a)) / (logisticcdf(b) - logisticcdf(a))

# ----------------- probit -----------------

# dlogcdf
cdf(     z::Real, ::Type{Val{:probit}}) = normcdf(z)
ccdf(    z::Real, ::Type{Val{:probit}}) = normccdf(z)
invcdf(  z::Real, ::Type{Val{:probit}}) = norminvcdf(z)
pdf(     z::Real, ::Type{Val{:probit}}) = isfinite(z) ?      normpdf(z) : zero(z)
dpdf(    z::Real, ::Type{Val{:probit}}) = isfinite(z) ? -z * normpdf(z) : zero(z)

logcdf(  z::Real, ::Type{Val{:probit}}) = normlogcdf(z)
logccdf( z::Real, ::Type{Val{:probit}}) = normlogccdf(z)
dlogcdf( z::Real, ::Type{Val{:probit}}) =   normpdf(z) / normcdf(z)
dlogccdf(z::Real, ::Type{Val{:probit}}) = - normpdf(z) / normccdf(z)

function dlogcdf_trunc(a::Real, b::Real, ::Type{Val{:probit}})
    a == typemin(a) && return   normpdf(b) / normcdf(b)
    b == typemax(b) && return - normpdf(a) / normccdf(a)
    return - sqrt(2/π) * _F1(a/sqrt2, b/sqrt2)
end

# d2logcdf
function d2logcdf(z::Real, ::Type{Val{:probit}})
    ϕΦ = normpdf(z) / normcdf(z)
    return - ϕΦ * (z + ϕΦ)
end
function d2logccdf(z::Real, ::Type{Val{:probit}})
    ϕomΦ = normpdf(z) / normccdf(z)
    return ϕomΦ * (z + ϕomΦ)
end


text:
using XLA, Test

@testset "XLA" begin

@testset "IR" begin
  include("ir.jl")
end

@testset "Compile" begin
  include("compile.jl")
end

@testset "Flux" begin
  include("flux.jl")
end

end