text:
#=##############################################################################
# DESCRIPTION
    Vehicles based on FLOWVLM (VLM and blade-element) models.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Apr 2020
  * License   : MIT
=###############################################################################

################################################################################
# ABSTRACT VLM VEHICLE TYPE
################################################################################
"""
    `AbstractVLMVehicle{N, M, R} <: AbstractVehicle{N, M, R}`

Implementations of `AbstractUVLMVehicle` are expected to have the following
fields:

# Required inputs
* `system::vlm.WingSystem`

# Optional inputs
* `tilting_systems::NTuple{N, vlm.WingSystem}`
* `rotor_systems::NTuple{M, Array{vlm.Rotor, 1}}`
* `vlm_system::vlm.WingSystem`
* `wake_system::vlm.WingSystem`
* `grids::Array{gt.GridTypes, 1}`

# Internal properties
* `V::Array{R, 1}`                          # Current vehicle velocity
* `W::Array{R, 1}`                          # Current vehicle angular velocity
* `prev_data::Array{Any, 1}`                # Information about previous step
* `grid_O::Array{Array{R, 1}, 1}`           # Origin of every grid

They also need to implement the following functions required by
`AbstractVehicle`:

* `shed_wake(...)`
* `generate_static_particle_fun(...)`
* `save_vtk(...)`

See the documentation and code of `AbstractVehicle` for more details.

In general, in order for implementation to work correctly, it is required that
the components of `wake_system` are also components of either `vlm_system` or
`rotor_systems`, and that none of the components of `rotor_systems` are also
components of `vlm_system`.
"""
abstract type AbstractVLMVehicle{N, M, R} <: AbstractVehicle{N, M, R} end

# Implementations
for header_name in ["unsteady", "quasisteady"]
  include("FLOWUnsteady_vehicle_vlm_"*header_name*".jl")
end

##### FUNCTIONS  ###############################################################
get_ntltsys(self::AbstractVLMVehicle) = typeof(self).parameters[1]

get_nrtrsys(self::AbstractVLMVehicle) = typeof(self).parameters[2]

function add_dV(self::AbstractVLMVehicle, dV)
    self.V .+= dV
    return nothing
end

function add_dW(self::AbstractVLMVehicle, dW)
    self.W .+= dW
    return nothing
end

function set_V(self::AbstractVLMVehicle, V)
    self.V .= V
    return nothing
end

function set_W(self::AbstractVLMVehicle, W)
    self.W .= W
    return nothing
end

function tilt_systems(self::AbstractVLMVehicle{N,M,R}, angles::NTuple{N, Array{R2, 1}}
                                                    ) where{N, M, R, R2<:Real}
    # Iterate over tilting system
    for i in 1:get_ntltsys(self)
        sys = self.tilting_systems[i]
        Oaxis = gt.rotation_matrix2([-a for a in angles[i]]...)
        vlm.setcoordsystem(sys, sys.O, Oaxis)
    end
    return nothing
end
function tilt_systems(self::AbstractVLMVehicle{0,M,R}, angles::Tuple{})  where{M, R}
    return nothing
end

function nextstep_kinematic(self::AbstractVLMVehicle, dt::Real)
    dX = dt*self.V                  # Translation
    dA = 180/pi * dt*self.W         # Angular rotation (degrees)

    O = self.system.O + dX          # New origin of the system (translation)
    M = gt.rotation_matrix2([-a for a in dA]...) # Rotation matrix
    Oaxis = M*self.system.Oaxis     # New axes of the system (rotation)

    # Save state of current vehicle
    _update_prev_vlm_system(self, deepcopy(self.vlm_system))
    _update_prev_wake_system(self, deepcopy(self.wake_system))
    _update_prev_rotor_systems(self, deepcopy(self.rotor_systems))

    # Translation and rotation
    vlm.setcoordsystem(self.system, O, Oaxis)

    # Translation and rotation of every grid
    for i in 1:length(self.grids)
        self.grid_O[i] .+= dX

        # Translation
        gt.lintransform!(self.grids[i], Array(1.0I, 3, 3) , dX)

        # Brings the grid back to the global origin
        gt.lintransform!(self.grids[i], Array(1.0I, 3, 3) , -self.grid_O[i])

        # Rotation and brings the grid back to its position
        gt.lintransform!(self.grids[i], M, self.grid_O[i])
    end

    return nothing
end

"""
Precalculations before calling the solver.

Calculates kinematic velocity on VLM an adds them as a solution field
"""
function precalculations(self::AbstractVLMVehicle, Vinf::Function,
                                pfield::vpm.ParticleField, t::Real, dt::Real,
                                nt::Int)

    if nt!=0
        # ---------- 1) Recalculate wake horseshoes with kinematic velocity -
        # Calculate kinematic velocity at every control point
        Vkin = _Vkinematic_wake(self, dt; t=t, targetX="CP")
        vlm._addsolution(self.wake_system, "Vkin", Vkin; t=t)

        Vkin = _Vkinematic_vlm(self, dt; t=t, targetX="CP")
        vlm._addsolution(self.vlm_system, "Vkin", Vkin; t=t)

        # Recalculate horseshoes
        vlm.getHorseshoes(self.vlm_system; t=t, extraVinf=_extraVinf1)
        vlm.getHorseshoes(self.wake_system; t=t, extraVinf=_extraVinf1)

        # ---------- 2) Paste previous Gamma solution ----------------------
        for i in 1:length(self.wake_system.wings)
            wing = vlm.get_wing(self.wake_system, i)
            prev_wing = vlm.get_wing(_get_prev_wake_system(self), i)

            if typeof(prev_wing) != vlm.Rotor
                sol = deepcopy(prev_wing.sol["Gamma"])
            else
                sol = deepcopy(prev_wing._wingsystem.sol["Gamma"])
            end

            vlm._addsolution(wing, "Gamma", sol; t=t)
        end
        for i in 1:length(self.vlm_system.wings)
            wing = vlm.get_wing(self.vlm_system, i)
            prev_wing = vlm.get_wing(_get_prev_vlm_system(self), i)

            if typeof(prev_wing) != vlm.Rotor
                sol = deepcopy(prev_wing.sol["Gamma"])
            else
                sol = deepcopy(prev_wing._wingsystem.sol["Gamma"])
            end

            vlm._addsolution(wing, "Gamma", sol; t=t)
        end
    end

    return nothing
end

function save_vtk_base(self::AbstractVLMVehicle, filename; path=nothing,
                                   num=nothing, save_wopwopin=false, optargs...)
    strn = vlm.save(self.system, filename; path=path, num=num, optargs...)

    for (i, grid) in enumerate(self.grids)
        strn *= gt.save(grid, filename*"_Grid$i"; path=path, num=num)
    end

    # Generate inputs for PSU-WOPWOP
    if save_wopwopin
        for (si, rotors) in enumerate(self.rotor_systems)
            # Compact patch for loading
            generate_vtkliftinglines(rotors, filename*"_Sys$(si)", path; num=num,
                                                                suf="_compact")
            # Loft for thickness
            # NOTE: this is not needed since FLOWNoise can read the lofted VTK
            # straight up
            # for (ri, rotors) enumerate(rotors)
            #     vlm.save(rotor, filename*"_Sys$(si)_Rotor$(ri)"; path=path, addtiproot=true,
            #                             wopwop=true, wopbin=true, wopv=1.0, num=num)
            # end

            # Loading file
            save_gammas = [[vlm.get_blade(rotor, i).sol["Gamma"]
                                      for i in 1:rotor.B] for rotor in rotors]
            # NOTE TO SELF: Ftot is a force per unit length
            save_Ftot = [[vlm.get_blade(rotor, i).sol["Ftot"]
                                      for i in 1:rotor.B] for rotor in rotors]
            loadfname = "loading_Sys$(si)"*(num!=nothing ? ".$num" : "")*".jld"
            JLD.save(joinpath(path, loadfname), "rotorgammas", save_gammas,
                                                        "Ftot", save_Ftot)
        end
    end

    return strn
end

##### INTERNAL FUNCTIONS  ######################################################

_get_prev_vlm_system(self::AbstractVLMVehicle) = self.prev_data[1]
_get_prev_wake_system(self::AbstractVLMVehicle) = self.prev_data[2]
_get_prev_rotor_systems(self::AbstractVLMVehicle) = self.prev_data[3]
function _update_prev_vlm_system(self::AbstractVLMVehicle, system)
    self.prev_data[1] = system
end
function _update_prev_wake_system(self::AbstractVLMVehicle, system)
    self.prev_data[2] = system
end
function _update_prev_rotor_systems(self::AbstractVLMVehicle, rotor_systems)
    self.prev_data[3] = rotor_systems
end

"""
Return the maximum number of static particles.
"""
function _get_m_static(self::AbstractVLMVehicle)
    m = vlm.get_m(self.vlm_system)

    for rotors in self.rotor_systems
        for rotor in rotors
            m += vlm.get_m(rotor)
        end
    end

    return m
end

"""
Returns the local translational velocity of every control point in `vlm_system`.
"""
function _Vkinematic_vlm(self::AbstractVLMVehicle, args...; optargs...)
    return _Vkinematic(self.vlm_system, _get_prev_vlm_system(self), args...;
                                                                    optargs...)
end

"""
Returns the local translational velocity of every control point in `wake_system`.
"""
function _Vkinematic_wake(self::AbstractVLMVehicle, args...; optargs...)
    return _Vkinematic(self.wake_system, _get_prev_wake_system(self), args...;
                                                                    optargs...)
end

function _Vkinematic(system, prev_system, dt::Real; t=0.0, targetX="CP")

    cur_Xs = _get_Xs(system, targetX; t=t)
    prev_Xs = _get_Xs(prev_system, targetX; t=t)

    return [-(cur_Xs[i]-prev_Xs[i])/dt for i in 1:size(cur_Xs, 1)]
end

"""
Returns the local translational velocity of every midpoint in the ri-th rotor in
the si-th system formatted as the output of `_get_midXs()`.
"""
function _Vkinematic_rotor(rotor_systems::NTuple{M, Array{vlm.Rotor, 1}},
                           prev_rotor_systems::NTuple{M, Array{vlm.Rotor, 1}},
                           si, ri, dt::Real
                          ) where{M}

    cur_Xs = _get_midXs(rotor_systems[si][ri])
    prev_Xs = _get_midXs(prev_rotor_systems[si][ri])

    Vkin = [-(cur_Xs[i]-prev_Xs[i])/dt for i in 1:size(cur_Xs, 1)]

    return Vkin

end

function _get_Xs(system::Union{vlm.Wing, vlm.WingSystem, vlm.Rotor},
                                                        targetX::String; t=0.0)
    if targetX=="CP"
        Xs = [vlm.getControlPoint(system, i
                                        ) for i in 1:vlm.get_m(system)]
    else
        Xs = [vlm.getHorseshoe(system, i; t=t)[vlm.VLMSolver.HS_hash[targetX]
                                        ] for i in 1:vlm.get_m(system)]
    end
    return Xs
end

function _get_Xs(system::Union{vlm.Wing, vlm.WingSystem, vlm.Rotor},
                                        targetXs::Array{T,1}; t=0.0) where{T}
    return vcat([_get_Xs(system, targetX; t=t) for targetX in targetXs]...)
end

"Returns the i-th bound-vortex midpoint of the j-th blade in `rotor`"
function _get_midX(rotor::vlm.Rotor, j::Int, i::Int)
    return (vlm.getHorseshoe(vlm.get_blade(rotor, j), i)[2] +
            vlm.getHorseshoe(vlm.get_blade(rotor, j), i)[3]) / 2
end

"Returns all bound-vortex midpoints of the j-th blade in `rotor`"
function _get_midXs(rotor::vlm.Rotor, j::Int)
    return [_get_midX(rotor, j, i)
                            for i in 1:vlm.get_m(vlm.get_blade(rotor, j)) ]
end
"Returns all bound-vortex midpoints in `rotor`"
_get_midXs(rotor::vlm.Rotor) = vcat([_get_midXs(rotor, j) for j in 1:rotor.B]...)

"Returns all bound-vortex midpoints in an array of rotors"
_get_midXs(rotors::Array{vlm.Rotor, 1}) = vcat([_get_midXs(rotor) for rotor in rotors]...)

"Returns all bound-vortex midpoints in multiple rotor systems"
function _get_midXs(rotor_systems::NTuple{M, Array{vlm.Rotor, 1}}) where{M}
    return vcat([_get_midXs(rotors) for rotors in rotor_systems]...)
end

"Receives an array `midXs` formatted as the output of `_get_midXs(rotor_systems)`
and returns the section of the array that corresponds to the ri-th rotor in the
si-th system (formatted linearly).
"
function _parse_midXs(rotor_systems::NTuple{M, Array{vlm.Rotor, 1}},
                        midXs::Array{Array{R, 1}}, si, ri) where{R<:Real, M}

    if si>length(rotor_systems) || si<=0
        error("Got invalid system index $si; max is $(length(rotor_systems))")
    end

    # Find lower bound of midXs belonging to rotor ri
    i_low = 1
    for (sj, rotors) in enumerate(rotor_systems[1:si])       # Iterate over rotor systems
        for (rj, rotor) in enumerate(rotors)                 # Iterative over rotors
            if sj==si && rj==ri
                break
            end
            i_low += vlm.get_m(rotor)
        end
    end

    # Find upper bound of midXs belonging to rotor ri
    rotor = rotor_systems[si][ri]
    i_up = i_low-1 + vlm.get_m(rotor)

    # # Format it as blades
    out = midXs[i_low:i_up]
    # out = [out[ ((i-1)*rotor.m+1):(i*rotor.m) ] for i in 1:rotor.B]
    return out
end

"""
    Receives the output of `parse_midXs(rotor_systems, midXs, si, ri)` and
outputs the array formatted as blades.
"""
function _format_blades(arr::Array{Array{R, 1}, 1},
                        rotor_systems::NTuple{M, Array{vlm.Rotor, 1}},
                        si, ri) where{R, M}

    rotor = rotor_systems[si][ri]
    arr = [arr[ ((i-1)*rotor.m+1):(i*rotor.m) ] for i in 1:rotor.B]

    return arr
end

"""
Function for FLOWVLM to fetch and add extra velocities to what every
control point imposes as the boundary condition
"""
function _extraVinf1(i, t; wing=nothing)
    if wing==nothing; error("Logic error!"); end;
    return wing.sol["Vkin"][i]
end
function _extraVinf2(i, t; wing=nothing)
    if wing==nothing; error("Logic error!"); end;
    return wing.sol["Vvpm"][i] + wing.sol["Vkin"][i]
end
function _extraVinf3(i, t; wing=nothing)
    if wing==nothing; error("Logic error!"); end;
    return wing.sol["Vind"][i] + wing.sol["Vkin"][i]
end

"""
Receives the FLOWVLM object `system` (Wing/WingSystem/Rotor), and adds vortex
particles to the particle field `pfield` at each trailing edge position where
an infinite vortex starts. `dt` indicates the length of the interval of time
that the vortex shedding represents.

Give it a previous system to detect differences in circulation and add
unsteady particles.
"""
function VLM2VPM(system::Union{vlm.Wing, vlm.WingSystem, vlm.Rotor}, pfield, dt,
                    Vinf;
                    t=0.0, prev_system=nothing, unsteady_shedcrit=-1.0,
                    p_per_step=1, sigmafactor=1.0, overwrite_sigma=nothing,
                    check=true, debug=false, tol=1e-6, omit_shedding=[])

  m = vlm.get_m(system)   # Number of lattices

  # Velocity at horseshoe points Ap and Bp
  if prev_system==nothing
      Vinfs_Ap = vlm.getVinfs(system; t=t, target="Ap")
      Vinfs_Bp = vlm.getVinfs(system; t=t, target="Bp")
  else
      Vinfs_Ap = []
      Vinfs_Bp = []
      for i in 1:m
          Ap = vlm.getHorseshoe(system, i)[1]
          prev_Ap = vlm.getHorseshoe(prev_system, i)[1]
          push!(Vinfs_Ap, -(Ap-prev_Ap)/dt + Vinf(Ap, t))

          Bp = vlm.getHorseshoe(system, i)[4]
          prev_Bp = vlm.getHorseshoe(prev_system, i)[4]
          push!(Vinfs_Bp, -(Bp-prev_Bp)/dt + Vinf(Bp, t))
      end
  end

  # Adds a particle at each infinite vortex
  prev_HS = [nothing for i in 1:8]
  for i in 1:m  # Iterates over lattices
    HS = vlm.getHorseshoe(system, i)

    Ap, A, B, Bp, CP, infDA, infDB, Gamma = HS
    (prev_Ap, prev_A, prev_B, prev_Bp, prev_CP,
      prev_infDA, prev_infDB, prev_Gamma) = prev_HS
    cntgs = true          # Contiguous horseshoes flag

    if nothing in HS; error("Logic error! $HS"); end;
    if true in [isnan.(elem) for elem in HS]; error("Logic error! $HS"); end;

    # ----------- Case of left wing tip ---------------------------
    if i==1
      # Adds particle at Ap
      X = Ap                                    # Particle position
      gamma = Gamma                             # Infinite vortex circulation
      V = norm(Vinfs_Ap[i])                     # Freestream at X
      infD = -infDA                             # Direction of vorticity
      sigma = sigmafactor*V*dt                  # Vortex blob radius
      vol = pi*(norm(Bp-Ap)/2)^2*V*dt           # Volume of particle
      l = -infD*V*dt                             # Distance the TE travels

      # Check if the system loops around
      if norm(Ap - vlm.getHorseshoe(system, m)[4]) / norm(Bp - Ap) <= tol
          gamma -= vlm.getHorseshoe(system, m)[8]
      end

      if unsteady_shedcrit<=0
          if !(i in omit_shedding)
              add_particle(pfield, X, gamma, dt, V, infD, sigma, vol,
                                l, p_per_step; overwrite_sigma=overwrite_sigma)
          end
      end

    # ----------- Case of wing tip on discontinuous wing --------
    elseif norm(Ap - prev_Bp) / norm(Bp - Ap) > tol
      if debug; println("Wing tip found at i=$i. $Ap!=$prev_Bp"); end;

      cntgs = false

      # Adds particle at previous Bp
      X = prev_Bp                               # Particle position
      gamma = -prev_Gamma                       # Infinite vortex circulation
      V = norm(Vinfs_Bp[i-1])                   # Freestream at X
      infD = -prev_infDB                        # Direction of vorticity
      sigma = sigmafactor*V*dt                  # Vortex blob radius
      vol = pi*(norm(prev_Bp-prev_Ap)/2)^2*V*dt # Volume of particle
      l = -infD*V*dt                             # Distance the TE travels

      if unsteady_shedcrit<=0
          if !(i-1 in omit_shedding)
              add_particle(pfield, X, gamma, dt, V, infD, sigma, vol,
                              l, p_per_step; overwrite_sigma=overwrite_sigma)
          end
      end

      # Adds particle at Ap
      X = Ap                                    # Particle position
      gamma = Gamma                             # Infinite vortex circulation
      V = norm(Vinfs_Ap[i])                     # Freestream at X
      infD = -infDA                             # Direction of vorticity
      sigma = sigmafactor*V*dt                  # Vortex blob radius
      vol = pi*(norm(Bp-Ap)/2)^2*V*dt           # Volume of particle
      l = -infD*V*dt                             # Distance the TE travels

      if unsteady_shedcrit<=0
          if !(i in omit_shedding)
              add_particle(pfield, X, gamma, dt, V, infD, sigma, vol,
                                l, p_per_step; overwrite_sigma=overwrite_sigma)
          end
      end

    # ----------- Case of contiguous horseshoes -----------------
    else

      # Verify logic
      if check
        crit1 = norm(Vinfs_Ap[i]-Vinfs_Bp[i-1]) / norm((Vinfs_Ap[i]+Vinfs_Bp[i-1])/2)
        crit2 = norm(infDA-prev_infDB) / norm((infDA+prev_infDB)/2)
        if crit1 > tol
          error("Logic error! Vinfs_Ap[i]!=Vinfs_Bp[i-1] "*
                  "( $(Vinfs_Ap[i])!=$(Vinfs_Bp[i-1]) )")
        elseif crit2 > tol
          error("Logic error! infDA!=prev_infDB "*
                  "( $(infDA)!=$(prev_infDB) )")
        elseif norm(infDA) - 1 > tol
          error("Logic error! norm(infDA)!= 1 "*
                  "( $(norm(infDA)) )")
        end
      end

      # Adds particle at Ap
      X = Ap                                    # Particle position
      gamma = Gamma - prev_Gamma                # Infinite vortex circulation
      V = norm(Vinfs_Ap[i])                     # Freestream at X
      infD = -infDA                             # Direction of vorticity
      sigma = sigmafactor*V*dt                  # Vortex blob radius
      vol = pi*(norm(Bp-Ap)/2)^2*V*dt           # Volume of particle
      l = -infD*V*dt                             # Distance the TE travels

      if unsteady_shedcrit<=0
          if !(i in omit_shedding)
              add_particle(pfield, X, gamma, dt, V, infD, sigma, vol,
                              l, p_per_step; overwrite_sigma=overwrite_sigma)
          end
      end

      # ----------- Case of discontinuous right wing tip -----------------------
      if i==m && norm(Bp - vlm.getHorseshoe(system, 1)[1]) / norm(Bp - Ap) > tol
        # Adds particle at Bp
        X = Bp                                    # Particle position
        gamma = -Gamma                            # Infinite vortex circulation
        V = norm(Vinfs_Bp[i])                     # Freestream at X
        infD = -infDB                             # Direction of vorticity
        sigma = sigmafactor*V*dt                  # Vortex blob radius
        vol = pi*(norm(Bp-Ap)/2)^2*V*dt           # Volume of particle
        l = -infD*V*dt                             # Distance the TE travels

        if unsteady_shedcrit<=0
            if !(i in omit_shedding)
                add_particle(pfield, X, gamma, dt, V, infD, sigma, vol,
                            l, p_per_step; overwrite_sigma=overwrite_sigma)
            end
        end
      end
    end

    # Adds unsteady particle
    if prev_system!=nothing && unsteady_shedcrit>0

      (p_Ap, p_A, p_B, p_Bp, p_CP,
          p_infDA, p_infDB, p_Gamma) = vlm.getHorseshoe(prev_system, i)
      X = (Ap+Bp)/2                               # LE midpoint position
      p_X = (p_Ap+p_Bp)/2                         # Previous LE midpoint position
      gamma = Gamma - p_Gamma                     # Bound circulation increase
      infD = (!cntgs || i==1 ? Ap : (prev_Ap+prev_Bp)/2) - X # Direction of circulation
      sigma = sigmafactor*norm(B-A)               # Vortex blob radius
      vol = 4/3*pi*(sigma/2)^3                    # Volume of particle
      l = -(X-p_X) + Vinf(X, t)*dt                # Distance the TE travels

      # Adds particle only if difference is greater than 1%
      if abs(gamma/p_Gamma) > unsteady_shedcrit && !(i in omit_shedding)
        add_particle(pfield, X, gamma, 1.0, 1.0, infD, sigma, vol,
                                        l, 1; overwrite_sigma=overwrite_sigma)
      end
    end

    prev_HS = HS
  end
end
##### END OF ABSTRACT VLM VEHICLE ##############################################


text:
module extras

"""
these methods are optional (coming from a C variant language)
a necssacity for those needing mental reflush
cleaner to return nothing in this case
There are 3 (equally valid) ways I might use to do this. All 3 return nothing:
"""

ext: jl

max_stars_repo_name: adamwillisXanax/SolvencyPredictor

max_issues_repo_name: adamwillisXanax/SolvencyPredictor

max_forks_repo_name: adamwillisXanax/SolvencyPredictor


text:
using Test
using ProxyInterfaces

# ProxyInterfaces.@indexing
# ----------------------

struct MyIndexing{T}
  value::T
end

ProxyInterfaces.indexing(::Type{MyIndexing{T}}) where T = T
ProxyInterfaces.indexing(myindexing::MyIndexing) = myindexing.value
ProxyInterfaces.@indexing MyIndexing

myindexing = MyIndexing([11,22,33])
# WARNING  == is intentionally not overwritten in order to enable the use of this proxy next to something else
# hence the default == applies, which compares by ref, and not by content.
@test !(myindexing == MyIndexing([11, 22, 33]))
@test myindexing[1] == 11
@test myindexing[end] == 33
@test firstindex(myindexing) == 1
@test lastindex(myindexing) == 3

# no mutable functionality is defined
@test_throws MethodError myindexing[2] = 2
@test_throws MethodError empty!(myindexing)
@test_throws MethodError fill!(myindexing, 0)

# no other methods
@test_throws MethodError 22 in myindexing
@test_throws MethodError length(myindexing) == 3
@test_throws MethodError size(myindexing) == (3,)

# ProxyInterfaces.@indexing_mutable
# ------------------------------

struct MyMutableIndexing{T}
  value::T
end

ProxyInterfaces.indexing(::Type{MyMutableIndexing{T}}) where T = T
ProxyInterfaces.indexing(myindexing::MyMutableIndexing) = myindexing.value
ProxyInterfaces.@indexing_mutable MyMutableIndexing

myindexing_mutable = MyMutableIndexing([12, 13])
myindexing_mutable[2] = 44
@test myindexing_mutable[2] == 44


text:
function all_tests(parts)
   display(parts)

   t = PArrays.PTimer(parts,verbose=true)
   PArrays.tic!(t)

   TestApp.GeometryTests.main(parts)
   PArrays.toc!(t,"Geometry")

   TestApp.CellDataTests.main(parts)
   PArrays.toc!(t,"CellData")

   TestApp.FESpacesTests.main(parts)
   PArrays.toc!(t,"FESpaces")

   TestApp.MultiFieldTests.main(parts)
   PArrays.toc!(t,"MultiField")

   TestApp.PoissonTests.main(parts)
   PArrays.toc!(t,"Poisson")

  TestApp.PLaplacianTests.main(parts)
  PArrays.toc!(t,"PLaplacian")

  TestApp.PeriodicBCsTests.main(parts)
  PArrays.toc!(t,"PeriodicBCs")

  display(t)
end


text:
abstract type Normalization end
struct RightNormalization <: Normalization end
struct LeftNormalization <: Normalization end
struct MixNormalization <: Normalization end

export normalization!,RightNormalization,LeftNormalization,MixNormalization,Normalization
"""
    using QR sweep mps from left to right or from right to left
    low level API, should not be exposed to user.
    Sweep from left to right
"""
function sweep!(mps::AbstractMPS,endpoint::Int,::LeftNormalization)
    res = 0.0
    for site =1:endpoint
        res+=single_tensor_sweep!(mps,site,LeftNormalization())
    end
    return res
end

"""
    using QR sweep mps from left to right or from right to left
    low level API, should not be exposed to user.
    Sweep from right to left
"""
function sweep!(mps::AbstractMPS,endpoint::Int,::RightNormalization)
    res = 0.0
    for site = mps.L:-1:endpoint
        res+=single_tensor_sweep!(mps,site,RightNormalization())
    end
    return res
end

"""
    using QR sweep mps from left to right or from right to left
    low level API, should not be exposed to user.
    Sweep from right to endpoint-1 and left to endpoint+1
"""
function sweep!(mps::AbstractMPS,endpoint::Int,::MixNormalization)
    res = 0.0
    res += sweep!(mps,endpoint-1,LeftNormalization())
    res += sweep!(mps,endpoint+1,RightNormalization())
    return res
end

"""
    Cut-off subprogram for a Left normalization MPS.
"""
function cutoff!(mps::AbstractMPS,Dcut::Int,::LeftNormalization)
    res = zeros(mps.L-1)
    for site = mps.L:-1:2
        res[site-1] = single_tensor_cutoff!(mps,site,Dcut,LeftNormalization())
    end
    return res
end

"""
    normalization for MPS. This function will set the norm of a mps to 1.
        return res = log(s...)
"""
function normalization!(mps::AbstractMPS,::LeftNormalization)
    res = sweep!(mps,mps.L-1,LeftNormalization())
    mps[mps.L] = mps[mps.L]/sqrt(sum(mps[mps.L].*mps[mps.L]))
    return res
end

function normalization!(mps::AbstractMPS,::RightNormalization)
    res = sweep!(mps,2,RightNormalization())
    mps[1] = mps[1]/sqrt(sum(mps[1].*mps[1]))
    return res
end

function normalization!(mps::AbstractMPS,site::Int,::MixNormalization)
    res = sweep!(mps,site,MixNormalization())
    mps[site] = mps[site]/sqrt(sum(mps[site].*mps[site]))
    return res
end

max_stars_repo_name: qiyang-ustc/JMPS

max_issues_repo_name: qiyang-ustc/JMPS

max_forks_repo_name: qiyang-ustc/JMPS


text:
export UniformPlane, Point, Rect, Circle, Hole

"""
Subtypes of `Hole` are holes in a plane.

**Subtypes**

- Point
- Rect
- Circle
"""
abstract type Hole <: Element end

"""
    Point(x, y, z)
    Point(<keyword arguments>)

`Point` objects are used for point holes in a `UniformPlane`.

Keyword arguments are `x`, `y`, `z`.
"""
struct Point <: Hole
  xyz1  :: Array{Float64,1}
end
Point(x,y,z) = Point([x,y,z,1])
Point(; x=0, y=0, z=0) = Point(x,y,z)

function Base.show(io::IO, x::Point)
  print(io,"+ hole point ")
  @printf(io,"(%.9e, %.9e, %.9e)",x.xyz1[1],x.xyz1[2],x.xyz1[3])
  println(io)
  return nothing
end

"""
    Rect(x1, y1, z1, x2, y2, z2)
    Rect(<keyword arguments>)

`Rect` objects are used for rectangular holes in a `UniformPlane`.

Keyword arguments are `x1`, `y1`, `z1`, `x2`, `y2`, `z2`.
"""
struct Rect <: Hole
  corner1  :: Array{Float64,1}
  corner2  :: Array{Float64,1}
end
Rect(x1,y1,z1,x2,y2,z2) = Rect([x1,y1,z1,1],[x2,y2,z2,1])
Rect(; x1=0, y1=0, z1=0, x2=0, y2=0, z2=0) =
  Rect(x1,y1,z1,x2,y2,z2)

function Base.show(io::IO, x::Rect)
  print(io,"+ hole rect ")
  @printf(io,"(%.9e, %.9e, %.9e, %.9e, %.9e, %.9e)",
          x.corner1[1],x.corner1[2],x.corner1[3],
          x.corner2[1],x.corner2[2],x.corner2[3])
  println(io)
  return nothing
end

"""
    Circle(x, y, z, r)
    Circle(<keyword arguments>)

`Circle` objects are used for circular holes in a `UniformPlane`.

Keyword arguments are `x`, `y`, `z`, `r`.
"""
struct Circle <: Hole
  xyz1    :: Array{Float64,1} # center
  radius  :: Float64
end
Circle(x,y,z,r) = Circle([x,y,z,1],r)
Circle(; x=0, y=0, z=0, r=0) =
  Circle(x,y,z,r)

function Base.show(io::IO, x::Circle)
  print(io,"+ hole circle ")
  @printf(io, "(%.9e, %.9e, %.9e, %.9e)",
          x.xyz1[1],x.xyz1[2],x.xyz1[3],
          x.radius)
  println(io)
  return nothing
end

"""
    UniformPlane(<keyword arguments>)

`UniformPlane` objects `show` a FastHenry uniform discretized plane command.

**Keyword Arguments**

- `name`              -- plane name
- `x1`,`y1`,`z1`      -- coordinate of first corner
- `x2`,`y2`,`z2`      -- coordinate of second corner
- `x3`,`y3`,`z3`      -- coordinate of third corner
- `thick`             -- segment thickness
- `seg1`,`seg2`       -- number of segments along each side of plane
- `segwid1`,`segwid2` -- segment width
- `sigma`,`rho`       -- specify conductivity or resistivity
- `nhinc`             -- integer number of filaments
- `rh`                -- filament ratio
- `relx`,`rely`,`relz`-- see FastHenry documentation
- `nodes`             -- Array of `Node` connection points to plane
- `holes`             -- Array of `Hole` (`Point`, `Rect`, or `Circle`)
"""
struct UniformPlane <: Element
  name :: Symbol
  corner1 :: Array{Float64,1}
  corner2 :: Array{Float64,1}
  corner3 :: Array{Float64,1}
  thick :: Float64
  seg1 :: Int
  seg2 :: Int
  segwid1 :: Float64
  segwid2 :: Float64
  sigma :: Float64
  rho :: Float64
  nhinc :: Int
  rh :: Int
  relx :: Float64
  rely :: Float64
  relz :: Float64
  nodes :: Array{Node,1}
  holes :: Array{Hole,1}
  function UniformPlane(name, corner1, corner2, corner3, thick, seg1, seg2,
                segwid1, segwid2, sigma, rho, nhinc, rh, relx, rely, relz,
                nodes::Array{Node,1}, holes)
    if isnan(thick)
      throw(ArgumentError("Must specify thick for plane"))
    end
    if seg1<1
      throw(ArgumentError("Must specify seg1>0 for plane"))
    end
    if seg2<1
      throw(ArgumentError("Must specify seg2>0 for plane"))
    end
    if nhinc<0
      throw(ArgumentError("nhinc must be a positive integer"))
    end
    if rh<0
      throw(ArgumentError("rh must be positive integer"))
    end
    if ~isnan(rho) && ~isnan(sigma)
      throw(ArgumentError("Cannot specify both rho and sigma"))
    end
    corner1 = convert(Array{Float64,1},corner1)
    corner2 = convert(Array{Float64,1},corner2)
    corner3 = convert(Array{Float64,1},corner3)
    thick = Float64(thick)
    v1 = corner1 - corner2
    v2 = corner3 - corner2
    if abs(dot(normalize(v1), normalize(v2))) > 1e-9
      throw(ArgumentError("Corners do not form a rectangle"))
    end
    for node in nodes
      plane = threepointstoabcd(corner1,corner2,corner3)
      pointtoplanedistance = distance(xyz(node),plane)
      if abs(pointtoplanedistance)>0.5*thick
        throw(ArgumentError("Node(s) not in plane"))
      end
    end
    new(Symbol(name), corner1, corner2, corner3, thick, seg1, seg2,
                segwid1, segwid2, sigma, rho, nhinc, rh, relx, rely, relz,
                nodes, holes)
  end
end
function UniformPlane(name, corner1, corner2, corner3, thick, seg1, seg2,
                segwid1, segwid2, sigma, rho, nhinc, rh, relx, rely, relz,
                nodes::Node, holes)
  UniformPlane(name, corner1, corner2, corner3, thick, seg1, seg2,
        segwid1, segwid2, sigma, rho, nhinc, rh, relx, rely, relz,
        [nodes], holes)
end
UniformPlane(;name = Symbol(""),
              x1=0, y1=0, z1=0, x2=0, y2=0, z2=0, x3=0, y3=0, z3=0,
              thick = NaN, seg1 = 0, seg2 = 0,
              segwid1 = NaN, segwid2 = NaN,
              sigma = NaN, rho = NaN,
              nhinc = 0, rh = 0,
              relx = NaN, rely = NaN, relz = NaN,
              nodes = Array{Node,1}(), holes = []) =
  UniformPlane(name,[x1,y1,z1,1],[x2,y2,z2,1],[x3,y3,z3,1], thick, seg1, seg2,
                segwid1, segwid2, sigma, rho, nhinc, rh, relx, rely, relz,
                nodes, holes)

function threepointstoabcd(c1::Array{Float64,1},
                           c2::Array{Float64,1},
                           c3::Array{Float64,1})
  n = cross(c1[1:3]-c2[1:3],c3[1:3]-c2[1:3]) # vector normal to plane
  d = dot(n,c2[1:3])
  return (n...,-d)
end

function distance(point::Array{Float64,1},
                  plane::Tuple{Float64,Float64,Float64,Float64})
  (a,b,c,d) = plane
  (a*point[1] + b*point[2] + c*point[3] + d)/(sqrt(a^2+b^2+c^2))
end


text:
immutable ChineseRestaurantProcess
	n::Int # Number of customers
	alpha::Float64
end

typealias CRP ChineseRestaurantProcess

function Distributions.rand(d::ChineseRestaurantProcess)
	s = [1]
	for t in 2:d.n
		psum = 0.0
		r = rand()
		table = 1
		n_tables = length(s)
		while table <= n_tables
			c = s[table]
			p = c / (t - 1 + d.alpha)
			psum += p
			if r < psum
				s[table] += 1
				break
			end
			table += 1
		end
		if table > n_tables
			push!(s, 1)
		end
	end
	# TODO: Always sort result to provide canonical form for partitions?
	return s
end

function Distributions.logpdf(d::ChineseRestaurantProcess, x::Vector)
	s = 0.0
	n_tables = length(x)
	n = 0
	for table in 1:n_tables
		c = x[table]
		n += 1
		p = d.alpha / (n - 1 + d.alpha)
		s += log(p)
		for i in 2:c
			n += 1
			p = (i - 1) / (n - 1 + d.alpha)
			s += log(p)
		end
	end
	return s
end

function Distributions.pdf(d::ChineseRestaurantProcess, x::Vector)
	return exp(logpdf(d, x))
end


text:
module sarscov2primers

import ArgParse

const version_info = v"0.1.0"
const prog_string = "primerpairs"
const version_string = "$(prog_string) command line interface $(version_info), by Dmitri Pavlichin"
const description_string = "Generate primer pairs for SARS-CoV-2 genomes"
const default_output_file = "out.csv"
const ref_genome_length = 29891 # SARS-CoV-2 reference genome length
const k = 25 # k-mer k

const arg_parse_settings = ArgParse.ArgParseSettings(
    prog = prog_string,
    description = join((version_string, description_string), "\n", "\n\n\n\n"),
    version = version_string,
)

ArgParse.@add_arg_table! arg_parse_settings begin
    "--coordinates", "-c"
        help = "pair of start and stop coordinates of target in SARS-CoV-2 reference genome"
        nargs = 2
        default = [k, ref_genome_length - k]
        arg_type = Int
    "--max-amplicon-length"
        help = "maximum amplicon length between forward and reverse primers (in the SARS-CoV-2 reference genome; length may be larger in other genomes if there are insertions)"
        default = 2500
        arg_type = Int
    "--min-amplicon-length"
        help = "minimum amplicon length between forward and reverse primers (in the SARS-CoV-2 reference genome; length may be shorter in other genomes if there are deletions)"
        default = 1
        arg_type = Int
    "--max-gc-diff"
        help = "maximum difference in GC content between forward and reverse primers"
        default = 2
        arg_type = Int
    "--output-file", "-o"
        help = "output file"
        default = default_output_file
        arg_type = String
    "--version", "-V"
        help = "display version number and exit"
        action = :show_version
end

# input files
#const sarscov2_4K_genomes_file = "../../genomes/sars-cov-2/sars-cov-2_4K_2020-4-9.fa"
const unique_conserved_kmers_file =          "kmer_counts/unique_conserved_25-mers_sars-cov-2_4K.bed"
const unique_conserved_specific_kmers_file = "kmer_counts/unique_conserved_specific_25-mers_sars-cov-2_4K.bed"
const input_files = (unique_conserved_kmers_file, unique_conserved_specific_kmers_file)

"""Return misisng expected input files, possibly uncompressing them"""
function missingsources()::Vector{String}
    missing_sources = String[]

    for file in input_files
        if !isfile(file)
            if isfile(file*".gz")
                @info "uncompressing $(file)"
                run(pipeline(`gunzip -k $(file)`))
            else
                push!(missing_sources, file)
            end
        end
    end

    return missing_sources
end

"""Reverse complement of DNA string `dna`"""
function rc(dna::AbstractString)
    replace(codeunits(dna), UInt8('A')=>UInt8('T'), UInt8('T')=>UInt8('A'), UInt8('C')=>UInt8('G'), UInt8('G')=>UInt8('C')) |> reverse .|> Char |> join
end

"""Canonical form of DNA string `dna`"""
canonical(dna::AbstractString) = min(dna, rc(dna))

function gccontent(kmer::AbstractString)
    count(x->x∈"CG", kmer)
end

# Return vector of start=>kmer
function primersfrombed(bed_file)::Dict{Int,String}
    @assert isfile(bed_file)
    split_lines = split.(readlines(bed_file), '\t')
    #primers = [(start=parse(Int, s[2]), stop=parse(Int, s[3]), kmer=String(s[4])) for s in split_lines]
    primers = Dict{Int,String}((parse(Int, s[2]) + 1)=>String(s[4]) for s in split_lines)
end

function findleastupperboundinsortedvec(v::DenseVector{T}, x::T) where T<:Real
    # returns index of least upper bound x in v
    # if x is greater than greatest element of v, return length(v)
    # assumes v is sorted in ascending order

    lower = firstindex(v)
    upper = lastindex(v)
    while lower < upper
        mid = div(upper+lower,2)
        if x <= v[mid]
            upper = mid
        else
            lower = mid+1
        end
    end
    lower
end

function findgreatestlowerboundinsortedvec(v::DenseVector{T}, x::T) where T<:Real
    # returns index of greatest lower bound x in v
    # if x is greater than greatest element of v, return 1
    # assumes v is sorted in ascending order

    lower = firstindex(v)
    upper = lastindex(v)
    while lower < upper
        mid = div(upper+lower,2)+1
        if x < v[mid]
            upper = mid-1
        else
            lower = mid
        end
    end
    lower
end

"""get primer pairs touching any position between `start` and `stop`, max
amplicon length `max_amplicon_length`, max GC content difference
`max_gc_diff`"""
function primerpairs(start::Integer, stop::Integer, max_amplicon_length::Integer, min_amplicon_length::Integer, max_gc_diff::Integer)

    if min_amplicon_length <= 0
        @warn "minimum amplicon length must be positive; using 1 instead"
        min_amplicon_length = 1
    end

    unique_conserved_primers_starts::Vector{Pair{Int,String}} = sort(collect(primersfrombed(unique_conserved_kmers_file)))

    starts = [x[1] for x in unique_conserved_primers_starts]
    f_primers = [x[2] for x in unique_conserved_primers_starts]

    unique_conserved_specific_primers_canonical::Set{String} = Set([canonical(kmer) for kmer in values(primersfrombed(unique_conserved_specific_kmers_file))])

    min_start = max(start - k - max_amplicon_length + 1, 0)
    max_stop = min(stop + max_amplicon_length, ref_genome_length)

    start_index = findleastupperboundinsortedvec(starts, min_start)
    stop_index = findgreatestlowerboundinsortedvec(starts, max_stop)

    # TODO: iterate i_right using fact that must start inside target
    primer_pairs = []
    for i_left in start_index:stop_index
        for i_right in (i_left + 1):stop_index
            f_start = starts[i_left]
            r_start = starts[i_right]
            amplicon_length = r_start - f_start - k
            amplicon_length < min_amplicon_length && continue
            amplicon_length > max_amplicon_length && break # sorted

            # check that overlap target
            ((r_start <= start) || (f_start + k - 1 >= stop)) && continue

            f_primer = f_primers[i_left]
            r_primer = f_primers[i_right]
            abs(gccontent(f_primer) - gccontent(r_primer)) > max_gc_diff && continue
            !((canonical(f_primer) ∈ unique_conserved_specific_primers_canonical) || (canonical(r_primer) ∈ unique_conserved_specific_primers_canonical)) && continue

            push!(primer_pairs, (f_start=f_start, r_start=r_start, f_primer=f_primer, r_primer=r_primer))
        end
    end

    return primer_pairs
end

function writeprimerpairs(file::AbstractString, primer_pairs)

    unique_conserved_specific_primers_canonical::Set{String} = Set([canonical(kmer) for kmer in values(primersfrombed(unique_conserved_specific_kmers_file))])

    @info "found $(length(primer_pairs)) primer pairs"
    @info "writing output to $(file)"

    open(file, "w") do io
        println(io, "forward (5' to 3'),reverse (5' to 3'),forward strand,reverse strand,forward start,forward stop,reverse start,reverse stop,amplicon length,forward GC content (%),reverse GC content (%),forward specific,reverse specific")
        for pp in primer_pairs
            f_primer = pp.f_primer
            r_primer = pp.r_primer
            f_specific = canonical(f_primer) ∈ unique_conserved_specific_primers_canonical
            r_specific = canonical(r_primer) ∈ unique_conserved_specific_primers_canonical
            amplicon_length = pp.r_start - pp.f_start - k
            println(io, join((f_primer, rc(r_primer), '+', '-', pp.f_start, pp.f_start + k - 1, pp.r_start, pp.r_start + k - 1, amplicon_length, gccontent(f_primer), gccontent(r_primer), f_specific, r_specific), ','))
        end
    end
    return
end

function main()
    # check have input files
    for missing_source in missingsources()
        error("missing input file: $(missing_source)")
    end

    parsed_args = ArgParse.parse_args(arg_parse_settings)

    coordinates = parsed_args["coordinates"]
    max_amplicon_length = parsed_args["max-amplicon-length"]
    min_amplicon_length = parsed_args["min-amplicon-length"]
    max_gc_diff = parsed_args["max-gc-diff"]
    #verbose = parsed_args["verbose"]
    output_file = parsed_args["output-file"]

    primer_pairs = primerpairs(coordinates[1], coordinates[2], max_amplicon_length, min_amplicon_length, max_gc_diff)
    writeprimerpairs(output_file, primer_pairs)
end

end # module


text:
using Base: depwarn

const INF_WARN = Ref(false)

function depwarn_timetype_inf()
    if !INF_WARN[]
        depwarn("`infinity` support for $(Dates.TimeType) is deprecated. Use `$InfExtendedTime` instead", :depwarn_timetype_inf)
        INF_WARN[] = true
    end
end


text:
################################################################################
#
#  Generate random examples
#
################################################################################

function _jordan_block(R, n::Int, a)
  z = a * identity_matrix(R, n)
  for i in 1:(n - 1)
    z[i + 1, i] = 1
  end
  return z
end

function _rand_block(R, n::Int, l::Int)
  Rx, x = PolynomialRing(R, "x", cached = false)
  f = rand(Rx, n:n, -100:100)
  if !iszero(f)
    setcoeff!(f, degree(f), one(R))
  end
  while !isirreducible(f) || degree(f) != n
    f = rand(Rx, n:n, -100:100)
    if !iszero(f)
      setcoeff!(f, degree(f), one(R))
    end
  end
  return companion_matrix(f^l)
end

function _random_elementary_operations!(a; type = :rows)
  tr = type == :columns
  @assert issquare(a)
  n = nrows(a)
  d = det(a)
  if n == 1
    return a
  end
  i = rand(1:n)
  j = rand(1:n)
  while j == i
    j = rand(1:n)
  end

  r = rand(base_ring(a), -10:10)
  for k in 1:n
    if tr
      a[k, i] = a[k, i] + r * a[k, j]
    else
      a[i, k] = a[i, k] + r * a[j, k]
    end
  end
  @assert d == det(a)
  return a
end

function _rand_poly_deg(R, n)
  Rx, x = PolynomialRing(R, "x", cached = false)
  f = rand(Rx, n:n, -100:100)
  if !iszero(f)
    setcoeff!(f, degree(f), one(R))
  end
  while !isirreducible(f) || degree(f) != n
    f = rand(Rx, n:n, -100:100)
    if !iszero(f)
      setcoeff!(f, degree(f), one(R))
    end
  end
  return f
end

function _random_matrix_special(R, n, deg)
  z = dense_matrix_type(R)[]
  f = _rand_poly_deg(R, deg)
  for i in 1:4
    l = rand(1:2)
    for j in 1:l
      push!(z, companion_matrix(f^i))
    end
    if sum(nrows.(z)) > n
      break
    end
  end
  return diagonal_matrix(z)
end

function _random_sln(R, n; num_op = 10)
  a = identity_matrix(R, n)
  for i in 1:num_op
    _random_elementary_operations!(a; type = rand() < 0.5 ? :rows : :columns)
  end
  return a
end

function _random_matrix(R, block_shape)
  matrices = dense_matrix_type(R)[]
  for (r, l, ll) in block_shape
    B = _rand_block(R, r, l)
    for j in 1:ll
      push!(matrices, B)
    end
  end
  if length(matrices) == 0
    return zero_matrix(R, 0, 0)
  end
  return diagonal_matrix(matrices)
end

function _similarity_test_setup(R, n; max_block = 4, same = false)
  block_shape = Tuple{Int, Int, Int}[]
  nn = n

  if same
    AA = _random_matrix_special(R, div(n, 3), rand(1:max_block))
    nn = n - nrows(AA)
  end

  while nn > 0
    #@show nn
    r = min(rand(1:nn), max_block)
    #@show r
    l = max(Int(floor(nn/r)), 1)
    if same
      nbl = rand(1:l)
    else
      nbl = l
    end
    #@show nbl
    nnn = nn - r * nbl
    ll = max(Int(floor(nnn/(r * nbl))), 1)
    nbll = rand(1:ll)
    #@show nbll
    push!(block_shape, (r, nbl, nbll))
    nn = nn - r * nbl *  nbll
  end
  b = block_shape
  #@show block_shape
  A = _random_matrix(R, block_shape)
  if same
    A = diagonal_matrix([A, AA])
  end

  z = _random_sln(R, nrows(A))
  @assert isone(det(z))
  zinv = inv(z)
  B = z * A * zinv
  return A, B
end


text:
"""
ДАНО: Робот - в юго-западном углу поля, на котором расставлено некоторое количество маркеров

РЕЗУЛЬТАТ: Функция вернула значение средней температуры всех замаркированных клеток
"""
N=0 #Количество замаркированных клеток
V=0 #Сумма температур
function average(r::Robot)
    side=Ost
    while isborder(r,Nord)==false
        liner(r,side)
        side=inverse(side)
        move!(r,Nord)
    end
    liner(r,side)
    println(V/N)
end

function temper(r::Robot) #Пересчёт суммы температур
    global V
    if ismarker(r)==true
        V=V+temperature(r)
    end
end

function counter(r::Robot) #Пересчёт количества маркеров
    global N
    if ismarker(r)==true
        N=N+1
    end
end

function liner(r::Robot,side::HorizonSide) #Движение по горизонтали
    while isborder(r,side)==false
        move!(r,side)
        temper(r)
        counter(r)
    end
end

function inverse(side::HorizonSide) #Возможно использование библиотечного файла
    HorizonSide(mod(Int(side)+2,4))
end
#= Использованы
глобальные переменные =#

max_stars_repo_name: alexgruzberg/mirea-progs

max_issues_repo_name: alexgruzberg/mirea-progs

max_forks_repo_name: alexgruzberg/mirea-progs


text:
    "reviewerCity": "",
    "reviewerState": ""
  }, {
    "reviewId": "",
    "reviewDate": {
      "$date": "2016-04-09T00:00:00.000Z"
    },
    "reviewTitle": "The location and garden around the hotel",
    "reviewText": "Cleanliness - toilets are not cleaned properly Water leakage from shower, taps are not tight, bath room flooring very old, frequent power cut, internet very bad connection The location and garden around the hotel",
    "reviewerCity": "",
    "reviewerState": ""
  }, {
    "reviewId": "",
    "reviewDate": {
      "$date": "2015-08-09T00:00:00.000Z"
    },
    "reviewTitle": "very gud hotal",
    "reviewText": "nothing Everythings",
    "reviewerCity": "",
    "reviewerState": "",
    "reviewRating": 10
  }, {
    "reviewId": "",
    "reviewDate": {
      "$date": "2016-07-01T00:00:00.000Z"
    },
    "reviewTitle": "Great property, big enough to spend time inside it.",
    "reviewText": "Food was costly here as expected from Hyatt hotel. If you are on leisure trip , you can take cheaper options like boutique hotels like shambaling in boudha, though shambaling's access road is not good but service is good. Great property, big enough to spend time inside it. Its near boudh stupa which is walkable from here. Good food options near boudh stupa",
    "reviewerCity": "",
    "reviewerState": ""
  }, {
    "reviewId": "",
    "reviewDate": {
      "$date": "2016-07-05T00:00:00.000Z"
    },
    "reviewTitle": "Sound of A/",
    "reviewText": "Sound of A/C",
    "reviewerCity": "",
    "reviewerState": ""
  }, {
    "reviewId": "",
    "reviewDate": {
      "$date": "2016-06-09T00:00:00.000Z"
    },
    "reviewTitle": "The best and only place to stay in Kathmandu .. Go for it",
    "reviewText": "Food options Big lawns swimming pool jogging tracks. The property is just awesome",
    "reviewerCity": "",
    "reviewerState": ""
  }, {
    "reviewId": "",
    "reviewDate": {
      "$date": "2016-09-11T00:00:00.000Z"
    },
    "reviewTitle": "Everything",
    "reviewText": "Nothing Everything",
    "reviewerCity": "",
    "reviewerState": "",
    "reviewRating": 10
  }, {
    "reviewId": "",
    "reviewDate": {
      "$date": "2016-03-27T00:00:00.000Z"
    },
    "reviewTitle": "bad",
    "reviewText": "it's not clean bed is not comfortable toilet is totally disgusting stains and not clean after the last guest nothing",
    "reviewerCity": "",
    "reviewerState": ""
  }, {
    "reviewId": "",
    "reviewDate": {
      "$date": "2016-02-04T00:00:00.000Z"
    },
    "reviewTitle": "Great hotel its fun in the shower.!",
    "reviewText": "The water in the showers kept flicking from very hot to cold every 15 to 20 seconds. I have stayed in this hotel 3 times now and in 3 different rooms and it has always been the same. You just can't enjoy your shower because your waiting for the next bust of piping hot water or freezing cold water. The Staff are always helpful and extremely polite.",
    "reviewerCity": "",
    "reviewerState": ""
  }, {
    "reviewId": "",
    "reviewDate": {
      "$date": "2016-07-30T00:00:00.000Z"
    },
    "reviewTitle": "THE place to stay in Kathmandu",
    "reviewText": "excellence in all aspects, room, staff, food, comfort: it is 'home' when in Kathmandu",
    "reviewerCity": "",
    "reviewerState": ""
  }, {
    "reviewId": "",
    "reviewDate": {
      "$date": "2015-10-13T00:00:00.000Z"
    },
    "reviewTitle": "Beautiful hotel only 15 mins to/from the airport and 10 min walk to Boudhanath",
    "reviewText": "The food was fantastic as was the service",
    "reviewerCity": "",
    "reviewerState": ""
  }, {
    "reviewId": "",
    "reviewDate": {
      "$date": "2016-08-09T00:00:00.000Z"
    },
    "reviewTitle": "Excellent hospitality",

max_stars_repo_name: bigdatakid/mongodb-demos

max_issues_repo_name: bigdatakid/mongodb-demos

max_forks_repo_name: bigdatakid/mongodb-demos


text:
@testset "Vec" begin
    @test Vec([1,2,3]) isa Vec{Int}

    v = Vec([3,2,1])

    @test length(v) === 3
    @test indices(v) === Base.OneTo(3)

    @test v[1] === 3
    @test v[2] === 2
    @test v[3] === 1

    @test v[v] isa Vec
    v2 = v[v]
    @test v2 == Vec([1,2,3])
end


text:
@testset "minimal_is_common_interval_test" begin
    @test iscommon(interval(-27.0,-27.0))
    @test iscommon(interval(-27.0, 0.0))
    @test iscommon(interval(0.0,0.0))
    @test iscommon(interval(-0.0,-0.0))
    @test iscommon(interval(-0.0,0.0))
    @test iscommon(interval(0.0,-0.0))
    @test iscommon(interval(5.0, 12.4))
    @test iscommon(interval(-0x1.FFFFFFFFFFFFFp1023, 0x1.FFFFFFFFFFFFFp1023))
    @test !iscommon(entireinterval())
    @test !iscommon(emptyinterval())
    @test !iscommon(interval(-Inf, 0.0))
    @test !iscommon(interval(0.0, Inf))
end

@testset "minimal_is_singleton_test" begin
    @test isthin(interval(-27.0,-27.0))
    @test isthin(interval(-2.0, -2.0))
    @test isthin(interval(12.0,12.0))
    @test isthin(interval(17.1, 17.1))
    @test isthin(interval(-0.0,-0.0))
    @test isthin(interval(0.0,0.0))
    @test isthin(interval(-0.0, 0.0))
    @test isthin(interval(0.0, -0.0))
    @test !isthin(emptyinterval())
    @test !isthin(entireinterval())
    @test !isthin(interval(-1.0, 0.0))
    @test !isthin(interval(-1.0, -0.5))
    @test !isthin(interval(1.0, 2.0))
    @test !isthin(interval(-Inf,-0x1.FFFFFFFFFFFFFp1023))
    @test !isthin(interval(-1.0,Inf))
end

@testset "minimal_is_member_test" begin
    @test in(-27.0, interval(-27.0,-27.0))
    @test in(-27.0, interval(-27.0, 0.0))
    @test in(-7.0, interval(-27.0, 0.0))
    @test in(0.0, interval(-27.0, 0.0))
    @test in(-0.0, interval(0.0,0.0))
    @test in(0.0, interval(0.0,0.0))
    @test in(0.0, interval(-0.0,-0.0))
    @test in(0.0, interval(-0.0,0.0))
    @test in(0.0, interval(0.0,-0.0))
    @test in(5.0, interval(5.0, 12.4))
    @test in(6.3, interval(5.0, 12.4))
    @test in(12.4, interval(5.0, 12.4))
    @test in(0.0, entireinterval())
    @test in(5.0, entireinterval())
    @test in(6.3, entireinterval())
    @test in(12.4, entireinterval())
    @test in(0x1.FFFFFFFFFFFFFp1023, entireinterval())
    @test in(-0x1.FFFFFFFFFFFFFp1023, entireinterval())
    @test in(0x1.0p-1022, entireinterval())
    @test in(-0x1.0p-1022, entireinterval())
    @test !in(-71.0, interval(-27.0, 0.0))
    @test !in(0.1, interval(-27.0, 0.0))
    @test !in(-0.01, interval(0.0,0.0))
    @test !in(0.000001, interval(0.0,0.0))
    @test !in(111110.0, interval(-0.0,-0.0))
    @test !in(4.9, interval(5.0, 12.4))
    @test !in(-6.3, interval(5.0, 12.4))
    @test !in(0.0, emptyinterval())
    @test !in(-4535.3, emptyinterval())
    @test !in(-Inf, emptyinterval())
    @test !in(Inf, emptyinterval())
    @test !in(NaN, emptyinterval())
    @test !in(-Inf, entireinterval())
    @test !in(Inf, entireinterval())
    @test !in(NaN, entireinterval())
end


text:
function gauss_rule{S<:Real,T<:Real}(monic_rc::Tuple{Vector{S}, Vector{T}})
  # W. GAUTSCHI, ORTHOGONAL POLYNOMIALS AND QUADRATURE, Electronic
  # Transactions on Numerical Analysis, Volume 9, 1999, pp. 65-76.

  α = monic_rc[1]
  β = monic_rc[2]
  # set up Jacobi matrix and compute eigenvalues
  J = SymTridiagonal(α, √β)
  #J = float(diagm(α) + diagm(√β,1) + diagm(√β,-1))
  x, V = eig(J)
  w = vec(V[1,:]'.^2)
  w = w / sum(w);

  # symmetrise
  if all(α==0)
      x=0.5*(x-reverse(x));
      w=0.5*(w+reverse(w));
  end
  return (x::Vector{Float64},w::Vector{Float64})
end

export gauss_rule


text:
include("load_cheaplogic.jl")

include("config.jl")
include("misc.jl")
include("utils.jl")
include("types.jl")
include("primitives.jl")

#include("reso.jl")

include("kptype.jl")
include("kpbase.jl")

include("kpcore.jl")

include("coreprint.jl")

include("kpreso.jl")

include("analyzer.jl")
include("repl.jl")
include("merge.jl")

println("kpbase is loaded")

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
using Test, MLJModels

@testset "metadata" begin
    @testset "reading and extracting model metadata" begin
        @test include("metadata.jl")
    end
    @testset "model search" begin
        @test include("model_search.jl")
    end
    @testset "loading model code" begin
        @test include("loading.jl")
    end
end

@testset "built-in models" begin
    @testset "Constant" begin
        @test include("builtins/Constant.jl")
    end
    @testset "Transformers" begin
        @test include("builtins/Transformers.jl")
    end
end

@testset "MultivariateStats  " begin
    @test include("MultivariateStats.jl")
end

@testset "DecisionTree       " begin
    @test include("DecisionTree.jl")
end

@testset "Clustering         " begin
    @test include("Clustering.jl")
end

@testset "GLM                " begin
    @test include("GLM.jl")
end

@testset "LIBSVM             " begin
    @test include("LIBSVM.jl")
end

@testset "NaiveBayes         " begin
    @test include("NaiveBayes.jl")
end

@testset "XGBoost" begin
    @test include("XGBoost.jl")
end

@testset "NearestNeighbors" begin
    @test include("NearestNeighbors.jl")
end


text:
# Imports Sundials and DifferentialEquations to access the CVODE solvers
# Imports DynamicBoundspODEsIneq to use DifferentialInequality integrator
# Imports DynamicBoundsBase for API and problem description functions.
using DynamicBoundspODEsIneq, DynamicBoundsBase, Sundials, DifferentialEquations

# Create problem
y0(u) = [1.2; 1.1]
function f!(dy, y, u, t)
    dy[1] = u[2]*y[1]*(one(typeof(u[1])) - y[2])*u[1]
    dy[2] = u[2]*y[2]*(y[1] - one(typeof(u[1])))
    nothing
end
tspan = (0.0, 1.0)

u_l = Float64[1.1; 2.95]
u_u = Float64[2.2; 3.05]

prob = DynamicBoundsBase.ODERelaxProb(f!, tspan, y0, u_l, u_u)
set!(prob, SupportSet([i for i in 0.0:0.01:1.0]))

# Creates the integrator specifying that relaxations and subgradients
# thereof should be computed.
integrator = DifferentialInequality(prob,
                                    calculate_relax = false,
                                    calculate_subgradient = false)

# defines a point in p at which to compute relaxations
ratio = rand(2)
pstar = u_l.*ratio .+ u_u.*(1.0 .- ratio)
setall!(integrator, ParameterValue(), pstar)

# Computes the relaxation and subgradients
relax!(integrator)

# Computes a trajectory of the pODEs at pstar
integrate!(integrator)


text:
using GNDquery
using Test

@testset "GNDquery.jl" begin
    # Write your tests here.
end


text:
using MechanicalSketch
import MechanicalSketch: °, background, sethue, O, PT, finish, Point
import MechanicalSketch: empty_figure, m, s, arrow, color_with_lumin, PALETTE

let
BACKCOLOR = color_with_lumin(PALETTE[8], 70)
function restart()
    empty_figure(filename = joinpath(@__DIR__, "test_8.png"))
    background(BACKCOLOR)
    sethue(PALETTE[5])
end
restart()
p = O + (0m, 0m)
arrow(p, 30m/s, 3m/s, α = 0°, backgroundcolor = BACKCOLOR)
arrow(p, 30m/s, 3m/s, α = 30°, backgroundcolor = BACKCOLOR)
arrow(p, 30m/s, 0m/s, α = 90°, backgroundcolor = BACKCOLOR)
arrow(p, 30m/s, 0m/s, α = 180°, backgroundcolor = BACKCOLOR)
arrow(p, 30m/s, 0m/s, α = 225°, backgroundcolor = BACKCOLOR)
arrow(p, 30m/s, 0m/s, α = 270°, backgroundcolor = BACKCOLOR)
arrow(p, 30m/s, 0m/s, α = 300°, backgroundcolor = BACKCOLOR)

for i = 1:10
    sethue(PALETTE[i])
    α = 36° * i
    p = O + (-8m, -2m) + Point(1m * cos(α), 1m * sin(α))
    arrow(p, i * 3m/s, 0m/s, α = α, backgroundcolor = BACKCOLOR,  labellength = true)
end

for i = 1:10
    sethue(PALETTE[i])
    α = 36° * i
    p = O + (8m, 4m) + Point(1m * cos(α), 1m * sin(α))
    arrow(p, 0m/s, i * 3m/s, α = α, backgroundcolor = BACKCOLOR, linewidth = 2PT)
end

for i = 1:20
    sethue(get(PALETTE, i / 20))
    α = 18° * i
    p = O + (8m, -4m) + Point(2m * cos(α), 2m * sin(α))
    arrow(p, 0.0m/s, i * 1.5m/s, α = α, backgroundcolor = BACKCOLOR, linewidth = 1.5PT)
end

finish()
nothing
end


text:
module gRPC

using ProtoBuf
using HTTP2
using HTTP2.Session

import Base: TCPServer, close, run
import ProtoBuf: call_method
import HTTP2: Headers
import HTTP2.Session: HTTPConnection, EvtGoaway, EvtRecvHeaders, EvtRecvData

include("common.jl")
include("server.jl")
include("client.jl")

export gRPCController
export gRPCClient, close, stub
export gRPCServer, close, run

end # module


text:
t = (3, "hello")
@assert typeof(t) == Tuple{Int64, String}
@assert isa(t, Tuple{Real, String})
@assert isa(t, Tuple{Integer, AbstractString})
@assert isa(t, Tuple{Number, Any})

@assert Tuple{Symbol, Float16, Union{}} <: Tuple{Any, Real, Any}

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
#### Formatter

function formatter()
    return "using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using JuliaFormatter

formatter_options = (
    indent = 4,
    margin = 80,
    always_for_in = true,
    whitespace_typedefs = true,
    whitespace_ops_in_indices = true,
    remove_extra_newlines = false,
)

filenames = readlines(`git diff --name-only --diff-filter=AM HEAD`)
filter!(f -> endswith(f, \".jl\"), filenames)

format(filenames; formatter_options...)"
end

function formatter_Project_toml()
    return "[deps]
JuliaFormatter = \"98e50ef6-434e-11e9-1051-2b60c6c9e899\"

[compat]
JuliaFormatter = \"0.3\"
"
end

max_stars_repo_name: charleskawczynski/PkgTemplates

max_issues_repo_name: charleskawczynski/PkgTemplates

max_forks_repo_name: charleskawczynski/PkgTemplates


text:
using SigMF
using Documenter

DocMeta.setdocmeta!(SigMF, :DocTestSetup, :(using SigMF); recursive=true)

makedocs(;
    modules=[SigMF],
    authors="Steve Kelly <kd2cca@gmail.com> and contributors",
    repo="https://github.com/juliaIO/SigMF.jl/blob/{commit}{path}#{line}",
    sitename="SigMF.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://juliaIO.github.io/SigMF.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Base Spec" => "sigmf-spec.md",
        "Signal Extension" => "signal.sigmf-ext.md",
    ],
)

deploydocs(;
    repo="github.com/juliaIO/SigMF.jl",
)


text:
# OedipusRex.jl
# - Polaron optical absorption, complex impedence and complex conductivity
# References:
# [1] R. P. Feynman, R. W. Hellwarth, C. K. Iddings, and P. M. Platzman, Mobility of slow electrons in a polar crystal, PhysicalReview127, 1004 (1962)
# [2] Devreese, J. De Sitter, and M. Goovaerts, “Optical absorption of polarons in thefeynman-hellwarth-iddings-platzman approximation,”Phys. Rev. B, vol. 5, pp. 2367–2381, Mar 1972.
# [3] F. Peeters and J. Devreese, “Theory of polaron mobility,” inSolid State Physics,pp. 81–133, Elsevier, 1984.

"""
----------------------------------------------------------------------
Polaron absorption coefficient Γ(Ω).
----------------------------------------------------------------------
"""

"""
    optical_absorption(Ω, β, α, v, w; rtol = 1e-3)

Calculate the absorption coefficient Γ(Ω) for the polaron at at finite temperatures (Eqn. (11a) in DSG 1972) for a given frequency Ω.

# Arguments
- `Ω::Float64`: is the frequency (2π THz) of applied electric field.
- `β::Float64`: is the reduced thermodynamic betas.
- `α::Float64`: is the Frohlich alpha coupling parameter.
- `v::Float64`: is the 'v' variational parameter.
- `w::Float64`: is the 'w' variational parameter.
- `rtol`: relative tolerance for the accuracy of any quadrature integrations.

See DSG 1972: https://doi.org/10.1103/PhysRevB.5.2367.
"""
function optical_absorption(Ω, β, α, v, w; rtol = 1e-3)
    real(complex_conductivity(Ω, β, α, v, w; rtol = rtol))
end

"""
----------------------------------------------------------------------
The Complex Impedence and Conductivity of the Polaron.
----------------------------------------------------------------------
"""

"""
    polaron_complex_impedence(Ω, β, α, v, w; rtol = 1e-3, T = nothing, verbose = false)

Calculate the complex impedence Z(Ω) of the polaron at finite temperatures for a given frequency Ω (Eqn. (41) in FHIP 1962).

# Arguments
- `Ω::Float64`: is the frequency (2π THz) of applied electric field.
- `β::Float64`: is the reduced thermodynamic betas.
- `α::Float64`: is the Frohlich alpha coupling parameter.
- `v::Float64`: is the 'v' variational parameter.
- `w::Float64`: is the 'w' variational parameter.
- `rtol`: relative tolerance for the accuracy of any quadrature integrations.
- `T`: is a token used by `make_polaron()` to keep track of the temperature for printing during a calculation. Do not alter.
- `verbose`: is used by `make_polaron()` to specify whether or not to print. Ignore.

See FHIP 1962: https://doi.org/10.1103/PhysRev.127.1004.

"""
function polaron_complex_impedence(Ω, β, α, v, w; ω = 1.0, rtol = 1e-3, T = nothing, verbose = false)
	impedance = -im * Ω * 2π + im * polaron_memory_function(Ω, β, α, v, w; ω = ω, rtol = rtol)

    if verbose
        println("\e[2K", "Process: $(count) / $processes ($(round.(count / processes * 100, digits = 1)) %) | T = $(round.(T, digits = 3)) | Ω = $(round.(Ω, digits = 3)) | Z = $(round.(impedance, digits = 3))")
        print("\033[F")

        global count += 1
    end

    return impedance
end

"""
    polaron_complex_conductivity(Ω, β, α, v, w; rtol = 1e-3)

Calculate the complex conductivity σ(Ω) of the polaron at finite temperatures for a given frequency Ω (equal to 1 / Z(Ω) with Z the complex impedence).

# Arguments
- `Ω::Float64`: is the frequency (2π THz) of applied electric field.
- `β::Float64`: is the reduced thermodynamic betas.
- `α::Float64`: is the Frohlich alpha coupling parameter.
- `v::Float64`: is the 'v' variational parameter.
- `w::Float64`: is the 'w' variational parameter.
- `rtol`: relative tolerance for the accuracy of any quadrature integrations.

See also [`polaron_complex_impedence`](@ref)
"""
function polaron_complex_conductivity(Ω, β, α, v, w; ω = 1.0, rtol = 1e-3)
	return 1 / polaron_complex_impedence(Ω, β, α, v, w; ω = ω, rtol = rtol)
end

"""
    polaron_mobility(β, α, v, w; rtol = 1e-3, T = nothing, verbose = false)

Calculate the dc mobility μ of the polaron at finite temperatues (Eqn. (11.5) in F. Peeters and J. Devreese 1984) for a given frequency Ω.

# Arguments
- `β::Float64`: is the reduced thermodynamic betas.
- `α::Float64`: is the Frohlich alpha coupling parameter.
- `v::Float64`: is the 'v' variational parameter.
- `w::Float64`: is the 'w' variational parameter.
- `rtol`: relative tolerance for the accuracy of any quadrature integrations.
- `T`: is a token used by `make_polaron()` to keep track of the temperature for printing during a calculation. Do not alter.
- `verbose`: is used by `make_polaron()` to specify whether or not to print. Ignore.

See F. Peeters and J. Devreese 1984: https://doi.org/10.1016/S0081-1947(08)60312-4.

See also [`polaron_complex_conductivity`](@ref)
"""
function polaron_mobility(β, α, v, w; ω = 1.0, rtol = 1e-3, T = nothing, verbose = false)

	mobility = abs(1 / imag(polaron_memory_function_dc(β, α, v, w; ω = ω, rtol = rtol)))

    if verbose
        println("\e[2K", "Process: $(count) / $processes ($(round.(count / processes * 100, digits = 1)) %) | T = $(round.(T, digits = 3)) | μ = $(round.(mobility, digits = 3))")
        print("\033[F")

        global count += 1
    end

    return mobility
end


text:
module Bri17

using LinearAlgebra
using Scapin.Grid
using Scapin.Elasticity

"""
    modal_strain_displacement!(b, n, N, h)

Compute the modal strain-displacement vector `b` in place, for the spatial frequency `n`.

This function returns `b`.

The grid is defined by its size, `N` and spacing, `h`. The spatial frequency is
defined by `n ∈ CartesianIndices(1:N[1], …, 1:N[d])`.

"""
function modal_strain_displacement!(
    b::AbstractVector{Complex{T}},
    n::CartesianIndex{d},
    N::NTuple{d,Int},
    h::NTuple{d,T},
) where {d,T<:Number}
    ∑α = zero(T)
    c = Array{T}(undef, d)
    s = Array{T}(undef, d)
    for i = 1:d
        α = π * (n[i] - 1) / N[i]
        ∑α += α
        c[i] = cos(α)
        s[i] = sin(α) / h[i]
    end

    b₀ = -2sin(∑α) + 2im * cos(∑α)
    if d == 2
        b[1] = b₀ * s[1] * c[2]
        b[2] = b₀ * c[1] * s[2]
    elseif d == 3
        b[1] = b₀ * s[1] * c[2] * c[3]
        b[2] = b₀ * c[1] * s[2] * c[3]
        b[3] = b₀ * c[1] * c[2] * s[3]
    else
        throw(ArgumentError("not implemented"))
    end
    return b
end

"""
    modal_strain_displacement(n, N, h)

Return the modal strain-displacement vector `b`.

See [`modal_strain_displacement!`](@ref) for a description of the parameters.

"""
function modal_strain_displacement(
    n::CartesianIndex{d},
    N::NTuple{d,Int},
    h::NTuple{d,T},
) where {d,T<:Number}
    return modal_strain_displacement!(Array{Complex{T}}(undef, d), n, N, h)
end

"""
    modal_stiffness!(K, n, N, h, C)

Compute the modal stiffness matrix `K` in place, for the spatial frequency `n`.

This function returns `K`.

The grid is defined by its size, `N` and spacing, `h`. The spatial frequency is
defined by `n ∈ CartesianIndices(1:N[1], …, 1:N[d])`. The (homogeneous)
constitutive material is specified by `C`.

!!! danger "Scaling of the modal stiffness matrix"

    The modal stiffness matrix introduced above differs from the matrix initially
    introduced in Ref. [DOI:10.1002/nme.5263](https://doi.org/10.1002/nme.5263)
    by a factor `h[1] * … * h[d]`. More precisely,

    ```
    K̂{Scapin} = h[1] * … * h[d] * K̂{10.1002/nme.5263}
    ```

    Such rescaling makes the relation between modal and nodal stiffness operators
    more natural (the former is the DFT of the latter).

"""
function modal_stiffness!(
    K::AbstractMatrix{Complex{T}},
    n::CartesianIndex{d},
    N::NTuple{d,Int},
    h::NTuple{d,T},
    C::Hooke{d,T},
) where {d,T<:Number}
    # In the notation of [Bri17, see Eq. (B.17)]
    #
    # φ[i] = φ(zᵢ) * hᵢ
    # χ[i] = χ(zᵢ) * hᵢ²
    # ψ[i] = ψ(zᵢ) * hᵢ
    #
    # Which simplifies the expression of Hₖ and allows to express the
    # nodal stiffness matrix as the inverse DFT of the modal stiffness
    # matrix.

    φ = Array{T}(undef, d)
    χ = Array{T}(undef, d)
    ψ = Array{T}(undef, d)
    for i = 1:d
        β = 2π * (n[i] - 1) / N[i]
        φ[i] = 2 * (1 - cos(β)) / h[i]
        χ[i] = (2 + cos(β)) / 3 * h[i]
        ψ[i] = sin(β)
    end

    scaling = C.μ / (1 - 2 * C.ν)
    if d == 2
        H₁₁ = φ[1] * χ[2]
        H₂₂ = χ[1] * φ[2]
        K_diag = C.μ * (H₁₁ + H₂₂)
        K[1, 1] = scaling * H₁₁ + K_diag
        K[1, 2] = scaling * ψ[1] * ψ[2]
        K[2, 1] = K[1, 2]
        K[2, 2] = scaling * H₂₂ + K_diag
    elseif d == 3
        H₁₁ = φ[1] * χ[2] * χ[3]
        H₂₂ = χ[1] * φ[2] * χ[3]
        H₃₃ = χ[1] * χ[2] * φ[3]
        K_diag = C.μ * (H₁₁ + H₂₂ + H₃₃)
        K[1, 1] = scaling * H₁₁ + K_diag
        K[1, 2] = scaling * ψ[1] * ψ[2] * χ[3]
        K[1, 3] = scaling * ψ[1] * χ[2] * ψ[3]
        K[2, 1] = K[1, 2]
        K[2, 2] = scaling * H₂₂ + K_diag
        K[2, 3] = scaling * χ[1] * ψ[2] * ψ[3]
        K[3, 1] = K[1, 3]
        K[3, 2] = K[2, 3]
        K[3, 3] = scaling * H₃₃ + K_diag
    else
        throw(ArgumentError("not implemented"))
    end
    return K
end

"""
    modal_stiffness(n, N, h, C)

Return the modal stiffness matrix `K`.

See [`modal_stiffness!`](@ref) for a description of the parameters.

"""
function modal_stiffness(
    n::CartesianIndex{d},
    N::NTuple{d,Int},
    h::NTuple{d,T},
    C::Hooke{d,T},
) where {T<:Number,d}
    return modal_stiffness!(Array{Complex{T}}(undef, d, d), n, N, h, C)
end

struct DiscreteGreenOperatorBri17{d,T}
    C::Hooke{d,T}
    N::NTuple{d,Int}
    h::NTuple{d,T}
end

Base.size(::DiscreteGreenOperatorBri17{d,T}) where {d,T} =
    ((d * (d + 1)) ÷ 2, (d * (d + 1)) ÷ 2)
Base.eltype(::Type{DiscreteGreenOperatorBri17{d,T}}) where {d,T} = T

function apply!(
    ε̂::AbstractVector{Complex{T}},
    Γ̂::DiscreteGreenOperatorBri17{d,T},
    n::CartesianIndex{d},
    τ̂::AbstractVector{Complex{T}},
) where {T<:Number,d}
    if all(Tuple(n) .== 1)
        ε̂ .= zero(eltype(ε̂))
        return
    end
    s = one(T) / √(2 * one(T))
    b̂ = modal_strain_displacement(n, Γ̂.N, Γ̂.h)
    conj_b̂ = conj.(b̂)
    K̂ = modal_stiffness(n, Γ̂.N, Γ̂.h, Γ̂.C)
    # TODO: use static arrays here!
    τ̂_conj_b̂ = Vector{Complex{T}}(undef, d)
    if d == 2
        #     ┌               ┐
        # τ̂ = │ τ̂₁       τ̂₃/√2│
        #     │ τ̂₃/√2    τ̂₂   │
        #     └               ┘
        τ̂_conj_b̂[1] = τ̂[1] * conj_b̂[1] + s * τ̂[3] * conj_b̂[2]
        τ̂_conj_b̂[2] = τ̂[2] * conj_b̂[2] + s * τ̂[3] * conj_b̂[1]
    elseif d == 3
        #     ┌                        ┐
        #     │ τ̂₁       τ̂₆/√2    τ̂₅/√2│
        # τ̂ = │ τ̂₆/√2    τ̂₂       τ̂₄/√2│
        #     │ τ̂₅/√2    τ̂₄/√2    τ̂₃/√2│
        #     └                        ┘
        τ̂_conj_b̂[1] = τ̂[1] * conj_b̂[1] + s * (τ̂[6] * conj_b̂[2] + τ̂[5] * conj_b̂[3])
        τ̂_conj_b̂[2] = τ̂[2] * conj_b̂[2] + s * (τ̂[6] * conj_b̂[1] + τ̂[4] * conj_b̂[3])
        τ̂_conj_b̂[3] = τ̂[3] * conj_b̂[3] + s * (τ̂[5] * conj_b̂[1] + τ̂[4] * conj_b̂[2])
    end
    # Do not include the `-` sign!
    û = prod(Γ̂.h) .* (K̂ \ τ̂_conj_b̂)
    if d == 2
        ε̂[1] = û[1] * b̂[1]
        ε̂[2] = û[2] * b̂[2]
        ε̂[3] = s * (û[1] * b̂[2] + û[2] * b̂[1])
    elseif d == 3
        ε̂[1] = û[1] * b̂[1]
        ε̂[2] = û[2] * b̂[2]
        ε̂[3] = û[3] * b̂[3]
        ε̂[4] = s * (û[2] * b̂[3] + û[3] * b̂[2])
        ε̂[5] = s * (û[3] * b̂[1] + û[1] * b̂[3])
        ε̂[6] = s * (û[1] * b̂[2] + û[2] * b̂[1])
    end
    return ε̂
end

function apply(
    Γ̂::DiscreteGreenOperatorBri17{d,T},
    n::CartesianIndex{d},
    τ̂::AbstractVector{Complex{T}},
) where {T<:Number,d}
    return apply!(Array{Complex{T}}(undef, size(Γ̂, 1)), n, τ̂)
end

export modal_strain_displacement!,
    modal_strain_displacement,
    modal_stiffness!,
    modal_stiffness,
    DiscreteGreenOperatorBri17,
    apply!,
    apply
end


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

using Actors, Test

@msg Msg A B C
@test supertype(A) == Msg
@test supertype(B) == Msg
@test supertype(C) == Msg

@msg D E F
@test supertype(D) == Any
@test supertype(E) == Any
@test supertype(F) == Any


text:
using Dates

struct SwapIndex{TP <: TenorPeriod, B <: BusinessCalendar, C <: BusinessDayConvention, DC <: DayCount, T <: TermStructure, II <: IborIndex} <: InterestRateIndex
  familyName::String
  tenor::TP
  fixingDays::Int
  currency::Currency
  fixingCalendar::B
  fixedLegTenor::Dates.Period
  fixedLegConvention::C
  fixedLegDayCount::DC
  discount::T
  iborIndex::II
  exogenousDiscount::Bool
  lastSwap::VanillaSwap
  lastFixingDate::Date

  SwapIndex{TP, B, C, DC, T, II}(familyName::String,
                                tenor::TP,
                                fixingDays::Int,
                                currency::Currency,
                                fixingCalendar::B,
                                fixedLegTenor::Dates.Period,
                                fixedLegConvention::C,
                                fixedLegDayCount::DC,
                                discount::T,
                                iborIndex::II,
                                exogenousDiscount::Bool = true) where {TP, B, C, DC, T, II} =
            new{TP, B, C, DC, T, II}(familyName, tenor, fixingDays, currency, fixingCalendar, fixedLegTenor, fixedLegConvention, fixedLegDayCount, discount,
                                    iborIndex, exogenousDiscount)
end

SwapIndex(familyName::String,
          tenor::TP,
          settlementDays::Int,
          currency::Currency,
          fixingCalendar::B,
          fixedLegTenor::Dates.Period,
          fixedLegConvention::C,
          fixedLegDayCount::DC, discount::T,
          iborIndex::II) where {TP <: TenorPeriod, B <: BusinessCalendar, C <: BusinessDayConvention, DC <: DayCount, T <: TermStructure, II <: IborIndex} =
          SwapIndex{TP, B, C, DC, T, II}(familyName, tenor, settlementDays, currency, fixingCalendar, fixedLegTenor, fixedLegConvention, fixedLegDayCount, discount, iborIndex)

function EuriborSwapIsdaFixA(tenor::TenorPeriod, forwardingTS::YieldTermStructure, discTS::YieldTermStructure)
  ibor = tenor.period.value > 1 ? euribor_index(TenorPeriod(Dates.Month(6)), forwardingTS) : euribor_index(TenorPeriod(Dates.Month(3)), forwardingTS)

  return SwapIndex("EuriborSwapIsdaFixA", tenor, 2, EURCurrency(), TargetCalendar(), Dates.Year(1), ModifiedFollowing(), BondThirty360(), discTS, ibor)
end


text:
using Dierckx
using Plots
using RollingFunctions
using LinearAlgebra
using Dierckx
using Plots
using RollingFunctions
using LinearAlgebra
using Dates

function ForwardStep(mtime)

  global Θv,TA,SST,ρA,ρO,ustar,LH,SH,E,LW,TAU,Qnet,SW
  global KAm,KAt,KOm,KOt
  global UA1,qA1,ΘA1,UO1,TO1
  global qA,ΘA,UA,UO,TO,SO

  ########################################################################
  # time
  ########################################################################

  epochdays=Dates.date2epochdays(Date(mtime))
  t=(Dates.datetime2epochms(mtime)/1000-epochdays*24*60*60)/3600;
  #julian day
  jd=Dates.datetime2julian(mtime)-Dates.datetime2julian(DateTime(year(mtime)));

  ########################################################################
  # radiation
  ########################################################################
  AOI=angle_of_incidence(lat,lon,jd,t);
  AOI<0 ? AOI=0 : AOI
  f=orad(ZOF)
  SW=(SunConstant*(1-albedo).*f.*AOI);

  ########################################################################
  # surface
  ########################################################################

  hl,hs,_,_,_,_,_,_,_,_,_,ustar,_,_,_,_=
  bulk(TA[1]+d2k,qA[1],abs(UA[1]-UO[1]),
        SST,ZAC[1],ZAC[1],ZAC[1],ZAC[1],ρA[1]);
  LH=-hl;
  SH=-hs;
  E=-hl/Av;
  ϵa=0.725+0.17*log10(sum(qA/1000*1.2*ΔZA)*100)
  LW=emissivity*sb*((SST+d2k).^4 .- ϵa.*(TAup).^4);
  TAU=ρA[1].*ustar.^2 .*sign(ustar);
  Qnet=SW[1]-LH-SH-LW;
  println(Qnet)

  ########################################################################
  # Atmospheric vertical diffusion
  ########################################################################
  Θv=ΘA.*(1 .+humid_fac*qA);
  KAm,KAt,γcq,γct,γcm=holtslag(Θv,UA,UA.*0,LH,SH,ustar,ρA[1],ZAF);
  KAm,KAt,γcq,γct,γcm=holtslag(Θv,UA,UA.*0,LH,SH,ustar,ρA[1],ZAF);
  #println("----",maximum(γct),"----")
  ########################################################################
  # Oceanic vertical diffusion
  ########################################################################
  ρO=ρO0.*(1 .-alpha.*(TO.-tref).+beta.*(SO.-sref));
  B=gravity_mks.*(ρO0.-ρO)./ρO0;
  KOm,KOt,_= large(TO,TO.*0,UO,VO,ρO,Qnet,0,ustar,ZOF,lat);

  ########################################################################
  # Atmospheric explicit time stepping
  ########################################################################
  TAUS=sign(UA[1]-UO[1]);

  UAst[1]=UA[1]-TAU*TAUS/(ΔZA*ρA[1])*ΔT-(1-Aimp)*(
      KAm[1].*(UA[1]-UA[2])/ΔZA^2*ΔT);
  UAst[end]=UA[end]-(1-Aimp)*(
      KAm[end]*(UA[end]-UAup)/ΔZA^2*ΔT);
  UAst[2:end-1]=UA[2:end-1]+(1-Aimp)*(
      KAm[1:end-2].*(UA[1:end-2]-UA[2:end-1])/ΔZA^2*ΔT-
      KAm[2:end-1].*(UA[2:end-1]-UA[3:end  ])/ΔZA^2*ΔT);

  qAst[1]=qA[1]-E/(ΔZA*ρO[1])*ΔT-
      (1-Aimp)*(KAt[1]*(qA[1]-qA[2])/ΔZA^2*ΔT+
      (-γcq[1].*KAt[1])/ΔZA*ΔT);
  qAst[end]=qA[end]-(1-Aimp)*(KAt[end]*(qA[end]-q0)/ΔZA^2*ΔT+
      (γcq[end-1].*KAt[end-1]-γcq[end].*KAt[end])/ΔZA*ΔT);
  qAst[2:end-1]=qA[2:end-1]+(1-Aimp)*(
      KAt[1:end-2].*(qA[1:end-2]-qA[2:end-1])/ΔZA^2*ΔT-
      KAt[2:end-1].*(qA[2:end-1]-qA[3:end  ])/ΔZA^2*ΔT+
      (γcq[1:end-2].*KAt[1:end-2]-γcq[2:end-1].*KAt[2:end-1])/ΔZA*ΔT);

  ΘAst[1]=ΘA[1]+(SH)/cpa/(ΔZA*qA[1])*ΔT-
        (1-Aimp)*(KAt[1]*(ΘA[1]-ΘA[2])/ΔZA^2*ΔT+
        (-γct[1].*KAt[1])/ΔZA*ΔT);
  ΘAst[end]=ΘA[end]-
        (1-Aimp)*(KAt[end]*(ΘA[end]-ΘA0)/ΔZA^2*ΔT+
        (γct[end-1].*KAt[end-1]-γct[end].*KAt[end])/ΔZA*ΔT);
  ΘAst[2:end-1]=ΘA[2:end-1]+(1-Aimp)*(
        KAt[1:end-2].*(ΘA[1:end-2]-ΘA[2:end-1])/ΔZA^2*ΔT-
        KAt[2:end-1].*(ΘA[2:end-1]-ΘA[3:end  ])/ΔZA^2*ΔT+
        (γct[1:end-2].*KAt[1:end-2]-γct[2:end-1].*KAt[2:end-1])/ΔZA*ΔT);

  ########################################################################
  # Oceanic explicit time stepping
  ########################################################################
  if cpl==1

      TOst[1]=TO[1]+(SW[1]-LH-SH-LW)/cp/(ΔZO*ρO[1])*ΔT-
          (1-Oimp)*(KOt[1]*(TO[1]-TO[2])/ΔZO^2*ΔT);
      TOst[end]=TO[end]+SW[end]/cp/(ΔZO*ρO[end])*ΔT+
          (1-Oimp)*(KOt[end-1].*(TO[end-1]-TO[end])/ΔZO^2*ΔT);
      TOst[2:end-1]=TO[2:end-1]+(1-Oimp)*(
          KOt[1:end-2].*(TO[1:end-2]-TO[2:end-1])/ΔZO^2*ΔT-
          KOt[2:end-1].*(TO[2:end-1]-TO[3:end  ])/ΔZO^2*ΔT)+
          SW[2:end-1]./cp./(ΔZO*ρO[2:end-1])*ΔT;

      UOst[1]=UO[1]+
            (TAU)*TAUS/(ρO[1])/ΔZO*ΔT-
            (1-Oimp)*(KOm[1]*(UO[1]-UO[2])/ΔZO^2*ΔT);
      UOst[end]=UO[end]+
            (1-Oimp)*(KOm[end-1].*(UO[end-1]-UO[end])/ΔZO^2*ΔT);
      UOst[2:end-1]=UO[2:end-1]+(1-Oimp)*(
            KOm[1:end-2].*(UO[1:end-2]-UO[2:end-1])/ΔZO^2*ΔT-
            KOm[2:end-1].*(UO[2:end-1]-UO[3:end  ])/ΔZO^2*ΔT);

  end

  if Aimp==0
      UA1=UAst; qA1=qAst; ΘA1=ΘAst;
  else
      du=-Aimp*KAm[1:end-1].*ΔT./ΔZA^2; dl=-Aimp*KAm[2:end].*ΔT./ΔZA^2;
      d=(1 .+Aimp.*ΔT./ΔZA^2 .*([KAm[1:end-1]; 0]+[0; KAm[2:end]]));
      A=inv(Tridiagonal(dl,d,du));
      UA1=A*UAst;
      du=-Aimp*KAt[1:end-1].*ΔT./ΔZA^2; dl=-Aimp*KAt[2:end].*ΔT./ΔZA^2;
      d=(1 .+Aimp.*ΔT./ΔZA^2 .*([KAt[1:end-1]; 0]+[0; KAt[2:end]]));
      A=inv(Tridiagonal(dl,d,du));
      qA1=A*qAst; ΘA1=A*ΘAst;
  end

  if Oimp==0
      UO1=UOst; TO1=TOst;
  else
      du=-Oimp*KOm[1:end-1].*ΔT./ΔZO^2; dl=-Oimp*KOm[2:end].*ΔT./ΔZO^2;
      d=(1 .+Oimp.*ΔT./ΔZO^2 .*([KOm[1:end-1]; 0]+[0; KOm[2:end]]));
      A=inv(Tridiagonal(dl,d,du));
      UO1=A*UOst;
      du=-Oimp*KOt[1:end-1].*ΔT./ΔZO^2; dl=-Oimp*KOt[2:end].*ΔT./ΔZO^2;
      d=(1 .+Oimp.*ΔT./ΔZO^2 .*([KOt[1:end-1]; 0]+[0; KOt[2:end]]));
      A=inv(Tridiagonal(dl,d,du));
      TO1=A*TOst;
  end
  SST=TO[1]

  TA=(ΘA1).*(PA./PA[1]).^(2/7).-d2k;
  ρA=PA./(Rgas.*(TA.+d2k));

  if moist==1
      qsat=saltsat*cvapor_fac*exp.(-cvapor_exp./(TA.+d2k))./ρA;
      dq=(qA.-qsat).*((qA1.-qsat).>0);
      dq=min.(dq,qA1);
      dTHETA=Av.*dq./cp; #[J/kg][j-1 kg K]
      #println([i maximum(dTHETA) maximum(dq)])
      ind=qA1.>qsat;
      qA1[ind]=qA1[ind]-dq[ind];
      #THETA[i,:]=THETA[i,: ]+dTHETA;
  end

  qA=qA1; ΘA=ΘA1; UA=UA1;
  UO=UO1; TO=TO1;
end


text:
function add_modifier(obj::AbstractMeasurement, arg0::EstimationModifier)
    return jcall(obj, "addModifier", void, (EstimationModifier,), arg0)
end

function compare_to(obj::ComparableMeasurement, arg0::ComparableMeasurement)
    return jcall(obj, "compareTo", jint, (ComparableMeasurement,), arg0)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function estimate(obj::AbstractMeasurement, arg0::jint, arg1::jint, arg2::Vector{SpacecraftState})
    return jcall(obj, "estimate", EstimatedMeasurement, (jint, jint, Vector{SpacecraftState}), arg0, arg1, arg2)
end

function get_base_weight(obj::AbstractMeasurement)
    return jcall(obj, "getBaseWeight", Vector{jdouble}, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_coordinates(::Type{AbstractMeasurement}, arg0::SpacecraftState, arg1::jint, arg2::jint)
    return jcall(AbstractMeasurement, "getCoordinates", TimeStampedFieldPVCoordinates, (SpacecraftState, jint, jint), arg0, arg1, arg2)
end

function get_date(obj::AbstractMeasurement)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_dimension(obj::AbstractMeasurement)
    return jcall(obj, "getDimension", jint, ())
end

function get_modifiers(obj::AbstractMeasurement)
    return jcall(obj, "getModifiers", List, ())
end

function get_observed_value(obj::AbstractMeasurement)
    return jcall(obj, "getObservedValue", Vector{jdouble}, ())
end

function get_parameters_drivers(obj::AbstractMeasurement)
    return jcall(obj, "getParametersDrivers", List, ())
end

function get_satellites(obj::AbstractMeasurement)
    return jcall(obj, "getSatellites", List, ())
end

function get_theoretical_standard_deviation(obj::AbstractMeasurement)
    return jcall(obj, "getTheoreticalStandardDeviation", Vector{jdouble}, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function is_enabled(obj::AbstractMeasurement)
    return jcall(obj, "isEnabled", jboolean, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function set_enabled(obj::AbstractMeasurement, arg0::jboolean)
    return jcall(obj, "setEnabled", void, (jboolean,), arg0)
end

function signal_time_of_flight(::Type{AbstractMeasurement}, arg0::TimeStampedFieldPVCoordinates, arg1::FieldVector3D, arg2::FieldAbsoluteDate)
    return jcall(AbstractMeasurement, "signalTimeOfFlight", RealFieldElement, (TimeStampedFieldPVCoordinates, FieldVector3D, FieldAbsoluteDate), arg0, arg1, arg2)
end

function signal_time_of_flight(::Type{AbstractMeasurement}, arg0::TimeStampedPVCoordinates, arg1::Vector3D, arg2::AbsoluteDate)
    return jcall(AbstractMeasurement, "signalTimeOfFlight", jdouble, (TimeStampedPVCoordinates, Vector3D, AbsoluteDate), arg0, arg1, arg2)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
function solve(
    instance::UniformMatroidInstance{Float64, Maximise},
    ::GreedyAlgorithm,
)
    # Algorithm: sort the weights, take the m largest ones, this is the
    # optimum solution.
    # Implementation: no need for sorting, partialsortperm returns the
    # largest items.
    items = collect(partialsortperm(instance.rewards, 1:(instance.m), rev=true))
    return UniformMatroidSolution(instance, items)
end


text:
using Test, BridgeSampling, Turing

## Test on bivariate standard Normal distribution
samples = rand(MvNormal([1.0 0.0; 0.0 1.0]), 10_000)

log_posterior(x) = -0.5 * x' * x

Bridge = bridgesampling(samples, log_posterior, [-Inf, -Inf], [Inf, Inf])
Analytical = log(2*pi)

@test isapprox(value(Bridge), Analytical, atol = 1e-3)


text:
#####
##### Building blocks for traversing a Hamiltonian deterministically, using the leapfrog
##### integrator.
#####

####
#### kinetic energy
####

"""
$(SIGNATURES)
Return kinetic energy `κ`, at momentum `p`.
"""
function kinetic_energy(κ::GaussianKineticEnergy{D,T,L}, p::PtrVector{D,T,L}, q = nothing) where {D,T,L}
    M⁻¹ = κ.M⁻¹.diag
    # @show M⁻¹
    # @show p
    ke = zero(T)
    @avx for d ∈ eachindex(p)
        pᵈ = p[d]
        ke += pᵈ * M⁻¹[d] * pᵈ
    end
    T(0.5) * ke
end
# kinetic_energy(κ::GaussianKineticEnergy, p, q = nothing) = dot(p, κ.M⁻¹ * p) / 2

"""
$(SIGNATURES)
Return ``p♯ = M⁻¹⋅p``, used for turn diagnostics.
"""
function calculate_p♯(sptr::StackPointer, κ::GaussianKineticEnergy, p::PtrVector{P,T,L}, q = nothing) where {P,T,L}
    M⁻¹ = κ.M⁻¹
    sptr, M⁻¹p = PtrVector{P,T}(sptr)
    @avx for l ∈ eachindex(M⁻¹p)
        M⁻¹p[l] = M⁻¹[l] * p[l]
    end
    sptr, M⁻¹p
end
function calculate_p♯(tree::Tree{P,T,L}, κ::GaussianKineticEnergy, p::PtrVector{P,T,L}, q = nothing) where {P,T,L}
    M⁻¹ = κ.M⁻¹.diag
    M⁻¹p = undefined_ρ♯( tree )
    @avx for l ∈ eachindex(M⁻¹p)
        M⁻¹p[l] = M⁻¹[l] * p[l]
    end
    M⁻¹p
end

# """
# $(SIGNATURES)
# Calculate the gradient of the logarithm of kinetic energy in momentum `p`.
# """
# ∇kinetic_energy(sptr::StackPointer, κ::GaussianKineticEnergy, p, q = nothing) = calculate_p♯(sptr, κ, p)

"""
$(SIGNATURES)
Generate a random momentum from a kinetic energy at position `q`.
"""
function rand_p(rng::AbstractRNG, sptr::StackPointer, κ::GaussianKineticEnergy{P,T,L}, q = nothing) where {P,T,L}
    W = κ.W
    sptr, r = PtrVector{P,T,L}(sptr)
    sptr, randn!(rng, r, W)
end
rand_p!(rng::VectorizedRNG.AbstractPCG, r::PaddedMatrices.AbstractMutableFixedSizeVector{P,T,L}, κ::GaussianKineticEnergy{P,T,L}, q = nothing) where {P,T,L} = randn!(rng, r, κ.W.diag)
rand_p(rng::AbstractRNG, κ::GaussianKineticEnergy{P,T,L}, q = nothing) where {P,T,L} = randn!(rng, FixedSizeVector{P,T,L}(undef), κ)

"""
$(SIGNATURES)
Evaluate log density and gradient and save with the position. Preferred interface for
creating `EvaluatedLogDensity` instances.
"""
function evaluate_ℓ!(sptr::StackPointer, ∇ℓq::PtrVector{P,T,L,false}, ℓ::AbstractProbabilityModel{P}, q::PtrVector{P,T,L,false}) where {P,T,L}
    ℓq = logdensity_and_gradient!(∇ℓq, ℓ, q, sptr)
    sp = reinterpret(Int, pointer(sptr, Float64))
    dlq = reinterpret(Int, pointer(∇ℓq))
    dq = reinterpret(Int, pointer(q))
    # @show sp
    # @show (sp - dlq) >> 3
    # @show (sp - dq) >> 3
    if isfinite(ℓq)
        EvaluatedLogDensity(q, ℓq, ∇ℓq)
    else
        EvaluatedLogDensity(q, oftype(ℓq, -Inf), q) # second q used just as a placeholder
    end
end

function evaluate_ℓ(sptr::StackPointer, ℓ::AbstractProbabilityModel{P}, q::PtrVector{P,T,L}) where {P,T,L}
    sptr2, ∇ℓq = PtrVector{P,T,L}(sptr)
    ℓq = logdensity_and_gradient!(∇ℓq, ℓ, q, sptr2)
    if isfinite(ℓq)
        sptr2, EvaluatedLogDensity(q, ℓq, ∇ℓq)
    else
        sptr, EvaluatedLogDensity(q, oftype(ℓq, -Inf), q) # second q used just as a placeholder
    end
end

"""
$(SIGNATURES)
Log density for Hamiltonian `H` at point `z`.
If `ℓ(q) == -Inf` (rejected), skips the kinetic energy calculation.
Non-finite values (incl `NaN`, `Inf`) are automatically converted to `-Inf`. This can happen
if
1. the log density is not a finite value,
2. the kinetic energy is not a finite value (which usually happens when `NaN` or `Inf` got
mixed in the leapfrog step, leading to an invalid position).
"""
function QuasiNewtonMethods.logdensity(H::Hamiltonian{D,T,L}, z::PhasePoint{D,T,L}) where {D,T,L}
    @unpack ℓq = z.Q
    isfinite(ℓq) || return oftype(ℓq, -Inf)
    K = kinetic_energy(H.κ, z.p)
    ℓq - (isfinite(K) ? K : oftype(K, Inf))
end

function calculate_p♯(stack, H::Hamiltonian{D,T,L}, z::PhasePoint{D,T,L}) where {D,T,L}
    calculate_p♯(stack, H.κ, z.p)
end

"""
    leapfrog(H, z, ϵ)
Take a leapfrog step of length `ϵ` from `z` along the Hamiltonian `H`.
Return the new phase point.
The leapfrog algorithm uses the gradient of the next position to evolve the momentum. If
this is not finite, the momentum won't be either, `logdensity` above will catch this and
return an `-Inf`, making the point divergent.
"""
function leapfrog(tree::Tree{P,T,L},
        H::Hamiltonian{P,T,L},
        z::PhasePoint{P,T,L}, ϵ::T
) where {P,L,T}
    @unpack ℓ, κ = H
    @unpack p, Q = z
    @unpack q, ∇ℓq = Q
#    @argcheck isfinite(Q.ℓq) "Internal error: leapfrog called from non-finite log density"
    # @show bitstring(z.flag)
    treeptr, flag = undefined_z(tree)
    # @show bitstring(flag)
    LT = L*sizeof(T) # counting on this being aligned.
    pₘ = PtrVector{P,T,L}(treeptr)
    q′ = PtrVector{P,T,L}(treeptr + LT)

    # @show ϵ
    # @show z
    # @show bitstring(z.flag)
    M⁻¹ = κ.M⁻¹.diag
    ϵₕ = T(0.5) * ϵ
    @avx for l ∈ eachindex(q)
        pₘₗ = p[l] + ϵₕ * ∇ℓq[l]
        pₘ[l] = pₘₗ
        q′[l] = q[l] + ϵ * M⁻¹[l] * pₘₗ
    end
    # @show q′
    # Variables that escape:
    # p′, Q′ (q′, ∇ℓq)
    Q′ = evaluate_ℓ!(tree.sptr, PtrVector{P,T,L}(treeptr + 2LT), H.ℓ, q′) # ∇ℓq is sorted second
    # @show Q′.q
    # isfinite(Q′.ℓq) || return PhasePoint(Q′, pₘ, flag)
    # p′ = pₘ # PtrVector{P,T,L}(sptr + 3LT)
    ∇ℓq′ = Q′.∇ℓq
    @avx for l ∈ eachindex(pₘ)
        pₘ[l] += ϵₕ * ∇ℓq′[l]
    end
    PhasePoint(Q′, pₘ, flag)
end
function leapfrog(
    sp::StackPointer,
    H::Hamiltonian{P,T,L},
    z::PhasePoint{P,T,L}, ϵ::T
) where {P,L,T}
    @unpack ℓ, κ = H
    @unpack p, Q = z
    @unpack q, ∇ℓq = Q

    sptr = pointer(sp, T)
    LT = L*sizeof(T) # counting on this being aligned.
    pₘ = PtrVector{P,T,L}(sptr)
    q′ = PtrVector{P,T,L}(sptr + LT)

    M⁻¹ = κ.M⁻¹.diag
    ϵₕ = T(0.5) * ϵ
    @avx for l ∈ eachindex(p)
        pₘₗ = p[l] + ϵₕ * ∇ℓq[l]
        pₘ[l] = pₘₗ
        q′[l] = q[l] + ϵ * M⁻¹[l] * pₘₗ
    end
    # Variables that escape:
    # p′, Q′ (q′, ∇ℓq)
    sp, Q′ = evaluate_ℓ(sp + 2LT, H.ℓ, q′) # ∇ℓq is sorted second
    # isfinite(Q′.ℓq) || return PhasePoint(Q′, pₘ)
    ∇ℓq′ = Q′.∇ℓq
    # p′ = pₘ # PtrVector{P,T,L}(sptr + 3LT)
    @avx for l ∈ eachindex(pₘ)
        pₘ[l] += ϵₕ * ∇ℓq′[l]
    end
    sp, PhasePoint(Q′, pₘ)
end


text:
module TrashData

include("util/util.jl")
include("matrix-data/triangle.jl")
include("dataframe/dataframe.jl")
include("interpolation/functions.jl")
include("interpolation/spline.jl")
include("interpolation/smoothing.jl")
include("table-data/merge.jl")
include("table-data/resample.jl")
include("io/file.jl")
include("io/utils.jl")

export ResampleStrategy

export data_size_symmetric_matrix, data_column_upper_matrix, data_column_lower_matrix
export data_symmetric_lowerTriangle_to_upperTriangle
export file_data_transform
export load_data, save_data, load_multiple_data
export round_data!, round_data_file

export createDataFrameWithColumns

export mirror_antisymmetric_data

export sigmoid_of, sigmoid

# ---
export convert_functionDataFrame_to_Functions
export merge_data_piece, merge_data_frames_sorted, join_data_files_sequencially
export merge_data_file_columns_to_main

export resample_data_frame, resample_data_file

export create_table_with_constants, create_table_with_constants!

export mpos, mvec

end # module


text:
module ReadStatTables

using DataValues: DataValueVector
using Dates
using PrettyTables: pretty_table
using ReadStat: read_data_file
using Tables

import DataAPI: refarray, unwrap
import PrettyTables: compact_type_str
import Tables: columnnames

export refarray, unwrap
export Date, DateTime
export columnnames

export LabeledValue,
       LabeledArray,
       LabeledVector,
       labels,

       ReadStatMeta,
       varlabels,
       varformats,
       val_label_keys,
       val_label_dict,
       filelabel,
       filetimestamp,
       fileext,
       ReadStatTable,
       getmeta,

       readstat

include("LabeledArrays.jl")
include("datetime.jl")
include("table.jl")
include("readstat.jl")

end


text:
input = joinpath(@__DIR__, "input")
answers = read(input, String)

groups = map(group -> split(group, "\n"), split(answers, "\n\n"))

p1 = sum(map(qs -> length(qs), map(group -> reduce(union, group), groups)))
p2 = sum(map(qs -> length(qs), map(group -> reduce(intersect, group), groups)))

@assert(p1 == 7120)
@assert(p2 == 3570)

println("-----------------------------------------------------------------------")
println("custom customs -- part one :: $p1")
println("custom customs -- part two :: $p2")
println("-----------------------------------------------------------------------")


text:
__precompile__(true)

module MathParse

import Base.string
import Base.show

export string
export show

function replacewithfunc(s::String, m::Regex, f)
	while (i = match(m, s)) != nothing
		s = replace(s, m => f(i.captures), count=1)
	end
	return s
end

function mparse(s::String)
	s = replace(s, r"\s+" => "")
	s = replace(s, r"(\d)\(" => s"\g<1>*(")
	s = replace(s, r"\)(\d)" => s")*\g<1>")
	s = replace(s, r"\)\(" => s")*(")
	s = replacewithfunc(s, r"\(([^()]*)\)", (t) -> mparse(String(t[1])))
	@show s

	r = [(r"(\-?\d+\.?\d*)/(\-?\d+\.?\d*)", /),
		(r"(\-?\d+\.?\d*)\*(\-?\d+\.?\d*)", *),
		(r"(\-?\d+\.?\d*)\+(\-?\d+\.?\d*)", +),
		(r"(\-?\d+\.?\d*)\-(\-?\d+\.?\d*)", -)]

	apply(t, f) = string(f(parse(Float64, t[1]), parse(Float64, t[2])))

	for i in r
		s = replacewithfunc(s, i[1], (t) -> apply(t, i[2]))
	end

	return s
end

ops = [(+, '+'), (-, '-'), (*, '*'), (/, '/'), (^, '^')]

mutable struct Mathobj
	op::Int
	left::Union{Mathobj, Float64}
	right::Union{Mathobj, Float64}
end

function treeifyRec(s::String)
	if length(s) == 0
		return Float64(0)
	end
	found = false
	op = 0
	split_ind = -1
	while op < length(ops) && !found
		op += 1
		par_count::Int = 0
		for i in (op == 5 ? (1 : length(s)) : (length(s) : -1 : 1))
			if s[i] == ')'
				par_count += 1
			elseif s[i] == '('
				par_count -= 1
			elseif par_count == 0 && s[i] == ops[op][2]
				found = true
				split_ind = i
				break
			end
		end
	end

	if op > length(ops) || split_ind == -1
		if s[1] == '('
			return treeifyRec(s[2 : end - 1])
		else
			return parse(Float64, s)
		end
	end

	return Mathobj(op, treeifyRec(s[1 : split_ind - 1]), treeifyRec(s[split_ind + 1 : end]))
end

function treeify(s::String)
	s = replace(s, r"\s+" => "")
	s = replace(s, r"(\d)\(" => s"\g<1>*(")
	s = replace(s, r"\)(\d)" => s")*\g<1>")
	s = replace(s, r"\)\(" => s")*(")
	treeifyRec(s)
end

eval(obj::Mathobj)::Float64 = ops[obj.op][1](eval(obj.left), eval(obj.right))

string(mo::Mathobj) = (ops[mo.op][2], string(mo.left), string(mo.right))

function show(io::IO, mo::Mathobj)
	print(io, string(mo))
end

mtparse(s::String)::Float64 = eval(treeify(s))

end


text:
using GeometricIntegrators
using LinearAlgebra
using Aqua

# Aqua.test_all(GeometricIntegrators)

# ambiguities are checked for Base, LinearAlgebra, etc., which fails
# Aqua.test_all(GeometricIntegrators; ambiguities=false)

# undefined_exports doesn't seem to work well with reexport
# Aqua.test_all(GeometricIntegrators; ambiguities=false, undefined_exports=false)

# stale_deps catches test-only dependencies
Aqua.test_all(GeometricIntegrators; ambiguities=false, undefined_exports=false, stale_deps=false)


text:
function Kprec( pw_grid::PWGrid, psi::Array{ComplexF64,2} )
    Ncol = size(psi)[2]
    G2 = pw_grid.G2
    Npoints = pw_grid.Npoints
    Kpsi = zeros( ComplexF64, size(psi) )
    for ic = 1:Ncol
        for ip = 1:Npoints
          Kpsi[ip,ic] = psi[ip,ic] / ( 1.0 + G2[ip] )
        end
    end
    return Kpsi
end

function Kprec(psi)
  return psi
end


text:
"""
    PERK

Module implementing parameter estimation via regression with kernels (PERK).

# Exports
- `EuclideanKernel`: Euclidean inner product kernel (for ridge regression
  instead of kernel ridge regression)
- `GaussianKernel`: Gaussian kernel used in kernel ridge regression
- `GaussianRFF`: Approximation of Gaussian kernel using random Fourier features
- `generatenoisydata`: Function for generating noisy data
- `perk`: Function for running PERK
"""
module PERK

import ForwardDiff # ForwardDiff.Dual, ForwardDiff.partials, ForwardDiff.value
using LinearAlgebra: I, Diagonal, norm
using Random: AbstractRNG
import Random # Random.GLOBAL_RNG
using Statistics: mean

include("utils.jl")
include("kernels.jl")
include("training.jl")
include("krr.jl")
include("estimation.jl")
include("holdout.jl")

export EuclideanKernel
export GaussianKernel
export GaussianRFF
export generatenoisydata
export perk

end


text:
using Base.Dates
using BusinessDays
using Base.Test

bd = BusinessDays

println("Perftests")

d0 = Date(2015, 06, 29) ; d1 = Date(2100, 12, 20)

cal_type = bd.Brazil()
cal_sym = :Brazil
cal_str = "Brazil"

@time bd.initcache(cal_type)
bdays(cal_type, d0, d1) # force JIT compilation
@time bdays(cal_type, d0, d1)
@time for i in 1:1000000 bdays(cal_type, d0, d1) end

bdays(cal_sym, d0, d1) # force JIT compilation
bdays(cal_str, d0, d1) # force JIT compilation
@time for i in 1:1000000 bdays(cal_sym, d0, d1) end
@time for i in 1:1000000 bdays(cal_str, d0, d1) end

#=
Perftests
  0.010197 seconds (54.70 k allocations: 1.100 MB)
  0.000003 seconds (9 allocations: 240 bytes)
  0.464481 seconds (5.00 M allocations: 76.294 MB, 0.79% gc time)
=#

N = 1000
d0vec = fill(d0, N)
d1vec = fill(d1, N)
cal_type_vec = fill(cal_type, N)
cal_sym_vec = fill(cal_sym, N)
cal_str_vec = fill(cal_str, N)

# Warmup
bd.initcache(bd.Brazil())
bd.initcache("USSettlement")
bd.initcache(:UKSettlement)
bd.bdays(cal_type, d0, d1)
bd.bdays(cal_sym, d0, d1)
bd.bdays(cal_str, d0, d1)
bd.bdays(cal_type, d0vec, d1vec)
bd.bdays(cal_sym, d0vec, d1vec)
bd.bdays(cal_str, d0vec, d1vec)
bd.bdays(cal_type_vec, d0vec, d1vec)
bd.bdays(cal_sym_vec, d0vec, d1vec)
bd.bdays(cal_str_vec, d0vec, d1vec)

println("type")
@time for i in 1:1000 bd.bdays(cal_type, d0, d1) end
@time bd.bdays(cal_type, d0vec, d1vec)
@time bd.bdays(cal_type_vec, d0vec, d1vec)
println("sym")
@time for i in 1:1000 bd.bdays(cal_sym, d0, d1) end
@time bd.bdays(cal_sym, d0vec, d1vec)
@time bd.bdays(cal_sym_vec, d0vec, d1vec)
println("str")
@time for i in 1:1000 bd.bdays(cal_str, d0, d1) end
@time bd.bdays(cal_str, d0vec, d1vec)
@time bd.bdays(cal_str_vec, d0vec, d1vec)

bd.cleancache()

println("initcache type")
@time bd.initcache(cal_type)
bd.cleancache()

println("initcache sym")
@time bd.initcache(cal_sym)
bd.cleancache()

println("initcache str")
@time bd.initcache(cal_str)
bd.cleancache()


text:
mutable struct Handle
    shape::Shape
    label::Label{String}
end

Handle( ; shape, label) = Handle(shape, label)