text:
function initialize!(integrator, cache::PDIRK44ConstantCache) end

@muladd function perform_step!(integrator, cache::PDIRK44ConstantCache, repeat_step=false)
  @unpack dt,uprev,u = integrator
  alg = unwrap_alg(integrator, true)
  @unpack nlsolver, tab = cache
  @unpack γs,cs,α1,α2,b1,b2,b3,b4 = tab

  if isthreaded(alg.threading)
    k2 = Array{typeof(u)}(undef,2)
    k1 = Array{typeof(u)}(undef,2)
    let nlsolver=nlsolver, u=u, uprev=uprev, integrator=integrator, cache=cache, dt=dt, repeat_step=repeat_step,
      k1=k1
      @threaded alg.threading for i in 1:2
        nlsolver[i].z = zero(u)
        nlsolver[i].tmp = uprev
        nlsolver[i].γ = γs[i]
        nlsolver[i].c = cs[i]
        markfirststage!(nlsolver[i])
        k1[i] = nlsolve!(nlsolver[i], integrator, cache, repeat_step)
      end
    end
    nlsolvefail(nlsolver[1]) && return
    nlsolvefail(nlsolver[2]) && return
    let nlsolver=nlsolver, u=u, uprev=uprev, integrator=integrator, cache=cache, dt=dt, repeat_step=repeat_step,
      k1=k1, k2=k2
      @threaded alg.threading for i in 1:2
        nlsolver[i].c = cs[2+i]
        nlsolver[i].z = zero(u)
        nlsolver[i].tmp = uprev + α1[i] * k1[1] + α2[i] * k1[2]
        k2[i] = nlsolve!(nlsolver[i], integrator, cache, repeat_step)
      end
    end
    nlsolvefail(nlsolver[1]) && return
    nlsolvefail(nlsolver[2]) && return
    integrator.u = uprev + b1 * k1[1] + b2 * k2[1] + b3 * k1[2] + b4 * k2[2]
  else
    _nlsolver = nlsolver[1]
    _nlsolver.z = zero(u)
    _nlsolver.tmp = uprev
    _nlsolver.γ = γs[1]
    _nlsolver.c = cs[1]
    markfirststage!(_nlsolver)
    k11 = nlsolve!(_nlsolver, integrator, γs[1]*dt, repeat_step)
    nlsolvefail(_nlsolver) && return
    _nlsolver.z = zero(u)
    _nlsolver.tmp = uprev
    _nlsolver.γ = γs[2]
    _nlsolver.c = cs[2]
    markfirststage!(_nlsolver)
    k12 = nlsolve!(_nlsolver, integrator, γs[2]*dt, repeat_step)
    nlsolvefail(_nlsolver) && return
    _nlsolver.z = zero(u)
    _nlsolver.tmp = uprev + α1[1] * k11 + α2[1] * k12
    _nlsolver.γ = γs[1]
    _nlsolver.c = cs[3]
    markfirststage!(_nlsolver)
    k21 = nlsolve!(_nlsolver, integrator, γs[1]*dt, repeat_step)
    nlsolvefail(_nlsolver) && return
    _nlsolver.z = zero(u)
    _nlsolver.tmp = uprev + α1[2] * k11 + α2[2] * k12
    _nlsolver.γ = γs[2]
    _nlsolver.c = cs[4]
    markfirststage!(_nlsolver)
    k22 = nlsolve!(_nlsolver, integrator, γs[2]*dt, repeat_step)
    nlsolvefail(_nlsolver) && return
    integrator.u = uprev + b1 * k11 + b2 * k21 + b3 * k12 + b4 * k22
  end
end

function initialize!(integrator, cache::PDIRK44Cache) end

@muladd function perform_step!(integrator, cache::PDIRK44Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  alg = unwrap_alg(integrator, true)
  @unpack nlsolver,k1,k2,tab = cache
  @unpack γs,cs,α1,α2,b1,b2,b3,b4 = tab
  if isthreaded(alg.threading)
    let nlsolver=nlsolver, u=u, uprev=uprev, integrator=integrator, cache=cache, dt=dt, repeat_step=repeat_step,
      k1=k1
      @threaded alg.threading for i in 1:2
        nlsolver[i].z .= zero(eltype(u))
        nlsolver[i].tmp .= uprev
        nlsolver[i].γ = γs[i]
        nlsolver[i].c = cs[i]
        markfirststage!(nlsolver[i])
        k1[i] .= nlsolve!(nlsolver[i], integrator, cache, repeat_step)
      end
    end
    nlsolvefail(nlsolver[1]) && return
    nlsolvefail(nlsolver[2]) && return
    let nlsolver=nlsolver, u=u, uprev=uprev, integrator=integrator, cache=cache, dt=dt, repeat_step=repeat_step,
      k1=k1, k2=k2
      @threaded alg.threading for i in 1:2
        nlsolver[i].c = cs[2+i]
        nlsolver[i].z .= zero(eltype(u))
        @.. broadcast=false nlsolver[i].tmp = uprev + α1[i] * k1[1] + α2[i] * k1[2]
        k2[i] .= nlsolve!(nlsolver[i], integrator, cache, repeat_step)
      end
    end
    nlsolvefail(nlsolver[1]) && return
    nlsolvefail(nlsolver[2]) && return
  else
    _nlsolver = nlsolver[1]
    _nlsolver.z .= zero(eltype(u))
    _nlsolver.tmp .= uprev
    _nlsolver.γ = γs[1]
    _nlsolver.c = cs[1]
    markfirststage!(_nlsolver)
    k1[1] .= nlsolve!(_nlsolver, integrator, cache, repeat_step)
    nlsolvefail(_nlsolver) && return
    _nlsolver.z .= zero(eltype(u))
    _nlsolver.tmp .= uprev
    _nlsolver.γ = γs[2]
    _nlsolver.c = cs[2]
    markfirststage!(_nlsolver)
    k1[2] .= nlsolve!(_nlsolver, integrator, cache, repeat_step)
    nlsolvefail(_nlsolver) && return
    _nlsolver.z .= zero(eltype(u))
    @.. broadcast=false _nlsolver.tmp = uprev + α1[1] * k1[1] + α2[1] * k1[2]
    _nlsolver.γ = γs[1]
    _nlsolver.c = cs[3]
    markfirststage!(_nlsolver)
    k2[1] .= nlsolve!(_nlsolver, integrator, cache, repeat_step)
    nlsolvefail(_nlsolver) && return
    _nlsolver.z .= zero(eltype(u))
    @.. broadcast=false _nlsolver.tmp = uprev + α1[2] * k1[1] + α2[2] * k1[2]
    _nlsolver.γ = γs[2]
    _nlsolver.c = cs[4]
    markfirststage!(_nlsolver)
    k2[2] .= nlsolve!(_nlsolver, integrator, cache, repeat_step)
    nlsolvefail(_nlsolver) && return
  end
  @.. broadcast=false u = uprev + b1 * k1[1] + b2 * k2[1] + b3 * k1[2] + b4 * k2[2]
end


text:
module TSPLIB
  using DataStructures
  using Match

  export TSP, readTSP, readTSPLIB, TSPLIB95_path, findTSP, Optimals

  struct TSP
    name::AbstractString
    dimension::Integer
    weight_type::AbstractString
    weights::Matrix
    nodes::Matrix
    Dnodes::Bool
    ffx::Function
    pfx::Function
    optimal::Float64
  end

  const TSPLIB95_path = joinpath(pkgdir(TSPLIB), "data", "TSPLIB95", "tsp")

  const tsp_keys = ["NAME",
                    "TYPE",
                    "COMMENT",
                    "DIMENSION",
                    "EDGE_WEIGHT_TYPE",
                    "EDGE_WEIGHT_FORMAT",
                    "EDGE_DATA_FORMAT",
                    "NODE_COORD_TYPE",
                    "DISPLAY_DATA_TYPE",
                    "NODE_COORD_SECTION",
                    "DEPOT_SECTION",
                    "DEMAND_SECTION",
                    "EDGE_DATA_SECTION",
                    "FIXED_EDGES_SECTION",
                    "DISPLAY_DATA_SECTION",
                    "TOUR_SECTION",
                    "EDGE_WEIGHT_SECTION",
                    "EOF"]

  include("reader.jl")
  include("distances.jl")
  include("utils.jl")
  include("fitness.jl")
  include("optimals.jl")

end # module

ext: jl


text:
###############################################################################
"""

# General

This module contains Julia code for performing extrinsic camera calibration, as
described in Schweitzer & Cowen, 2021 (references below). It follows the
methodology develocped by Holland et al. (1997).

# References
K. T. Holland, R. A. Holman, T. C. Lippmann, J. Stanley and N. Plant, "Practical
use of video imagery in nearshore oceanographic field studies," in IEEE Journal
of Oceanic Engineering, vol. 22, no. 1, pp. 81-92, Jan. 1997, doi:10.1109/48.557542.

Schweitzer, S. A., & Cowen, E. A. (2021). "Instantaneous river-wide water surface
velocity field measurements at centimeter scales using infrared quantitative
image velocimetry". Water Resources Research, 57,
e2020WR029279. https://doi.org/10.1029/2020WR029279

# License
This software is released under the MIT license, however, the authors request to
be credited if you use the software, preferably by citing the above-mentioned
Schweitzer & Cowen (2021) paper.

MIT License

Copyright (c) 2022 Seth A. Schweitzer

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

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

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

"""
module ExtrinsicCameraCalibration
###############################################################################

###############################################################################
using Statistics
using Printf # This modules is only required to enable the optional text output of HollandEtAlExtrinsicCalibration()
###############################################################################

###############################################################################
export IntrinsicCalibrationParameters, ExtrinsicCalibrationParameters, # structs
    HollandEtAlExtrinsicCalibration, CalculateDLTcoefficients # functions
###############################################################################

###############################################################################
"""
    IntrinsicCalibrationParameters

Data type for intrinsic camera calibration parameters.
"""
struct IntrinsicCalibrationParameters <: AbstractFloat
    u₀::Float64                 # u-coordinate (horizontal) of the image center (pixels)
    v₀::Float64                 # v-coordinate (vertical) of the image center (pixels)
    λᵤ::Float64                 # horizontal magnification factor
    λᵥ::Float64                 # vertical magnification factor
    PixelPitch::Float64
    radial_distortion_poly_degrees::Array{Float64,1}
    radial_distortion_poly_coeffs::Array{Float64,1}
end
###############################################################################

###############################################################################
"""
    ExtrinsicCalibrationParameters

Data type for extrinsic calibration parameters.
"""
struct ExtrinsicCalibrationParameters <: AbstractFloat
    φ::Float64                  # azimuth. clockwise from the north (y-axis), radians.
    τ::Float64                  # tilt. angle away from the vertical of the camera centeral ray, radians.
    σ::Float64                  # roll. rotation of the camera around the centeral ray, radians.
    f::Float64                  # distance b/w FPA and optical center of camera, pixels. Note that this will often be a negative number (since the FPA is behind the focal point)
    x_c::Float64                # x coordinate of camera's optical center, meters
    y_c::Float64                # y coordinate of camera's optical center, meters
    z_c::Float64                # z coordinate of camera's optical center, meters
end
### Convenience constructor:
# Input as NamedTuple:
ExtrinsicCalibrationParameters(x::NamedTuple) = ExtrinsicCalibrationParameters(x.φ, x.τ, x.σ, x.f, x.x_c, x.y_c, x.z_c)
###############################################################################

##############################################################################
struct RotationMatrix <: AbstractFloat
    m11::Float64
    m12::Float64
    m13::Float64
    m21::Float64
    m22::Float64
    m23::Float64
    m31::Float64
    m32::Float64
    m33::Float64
end
""" Initialize a RotationMatrix from a set of ExtrinsicCalibrationParameters """
function RotationMatrix(χ::ExtrinsicCalibrationParameters)
    return RotationMatrix(f_m11(χ), f_m12(χ), f_m13(χ),
                          f_m21(χ), f_m22(χ), f_m23(χ),
                          f_m31(χ), f_m32(χ), f_m33(χ))
end
##############################################################################

###############################################################################
"""
    CalculateDLTcoefficients(χ, ι)

Calculates the DLT coefficients from the camera calibration coefficients

Uses the equations in the appendix of Holland et al. 1997 to calculate the
Direct Linear Transformation coefficients directly from the physical parameters
of the camera's installation - location, orientation, magnification.
"""
function CalculateDLTcoefficients(χ::ExtrinsicCalibrationParameters, ι::IntrinsicCalibrationParameters)::NamedTuple
    Rot = RotationMatrix(χ)
    L = fill(NaN, 11)
    LL = -( χ.x_c * Rot.m31 + χ.y_c * Rot.m32 + χ.z_c * Rot.m33 )
    abs(LL) < 1.0 ? @warn("LL has a small value! - results might be unstable.") : nothing
    L₁ = ( ι.u₀ * Rot.m31 -χ.f * Rot.m11 ) / (ι.λᵤ * LL)
    L₂ = ( ι.u₀ * Rot.m32 -χ.f * Rot.m12 ) / (ι.λᵤ * LL)
    L₃ = ( ι.u₀ * Rot.m33 -χ.f * Rot.m13 ) / (ι.λᵤ * LL)
    L₄ = -(L₁ * χ.x_c + L₂ * χ.y_c + L₃ * χ.z_c)
    L₅ = ( ι.v₀ * Rot.m31 -χ.f * Rot.m21 ) / (ι.λᵥ * LL)
    L₆ = ( ι.v₀ * Rot.m32 -χ.f * Rot.m22 ) / (ι.λᵥ * LL)
    L₇ = ( ι.v₀ * Rot.m33 -χ.f * Rot.m23 ) / (ι.λᵥ * LL)
    L₈ = -(L₅ * χ.x_c + L₆ * χ.y_c + L₇ * χ.z_c)
    L₉ = Rot.m31 / LL
    L₁₀ = Rot.m32 / LL
    L₁₁ = Rot.m33 / LL
    return (; L₁, L₂, L₃, L₄, L₅, L₆, L₇, L₈, L₉, L₁₀, L₁₁)
end
###############################################################################

###############################################################################
"""
    HollandEtAlExtrinsicCalibration(GCPuv, GCPxyz, χ, ι)

Returns the extrinsic camera calibration coefficients, calculated from a set of
GCPs in real-world and pixel-space coordinates, using the methodology of Holland
et al. 1997.
"""
function HollandEtAlExtrinsicCalibration(GCPuv::Array{Float64,2}, GCPxyz::Array{Float64,2},
                                         χinitial::ExtrinsicCalibrationParameters,
                                         ι::IntrinsicCalibrationParameters;
                                         MaxIterations = 50,
                                         PrintoutStats=true)
    # in practical experience if the method doesn't converge after about 30
    # iterations it probably will not converge within 101 either.
    χ = deepcopy(χinitial)      # this prevents modification of the original input
    nGCPs = minimum([size(GCPuv,1), size(GCPxyz,1)])
    A = fill(NaN, nGCPs*2,7)
    b = fill(NaN, nGCPs*2)
    CorrectionTerm = fill(Inf, 7)
    PrevCorrectionTerm = copy(CorrectionTerm)
    iteration,ContinueFlag = (1,true)
    while (ContinueFlag==true)
        for gcp in 1:nGCPs
            x,y,z = GCPxyz[gcp,:]
            u,v = GCPuv[gcp,:]
            F₀, G₀, Matrix_row1, Matrix_row2 = UpdateConstants(x, y, z, u, v, χ, ι)
            A[gcp,:]       = Matrix_row1
            A[gcp+nGCPs,:] = Matrix_row2
            b[gcp]         = -F₀
            b[gcp+nGCPs]   = -G₀
        end
        CorrectionTerm = A\b
        χ = ExtrinsicCalibrationParameters(
            χ.φ + CorrectionTerm[1],
            χ.τ + CorrectionTerm[2],
            χ.σ + CorrectionTerm[3],
            χ.f + CorrectionTerm[4],
            χ.x_c + CorrectionTerm[5],
            χ.y_c + CorrectionTerm[6],
            χ.z_c + CorrectionTerm[7],
        )
        if (iteration > MaxIterations)
            @warn join(["exceeded maximum  iterations ($iteration)",
                        "max(|ε|)=$(map(y -> @sprintf("%.0e",y),maximum(abs.(CorrectionTerm))))",
                        "rms(ε₁-ε₀)=$(map(y -> @sprintf("%.0e",y),sqrt(mean(CorrectionTerm-PrevCorrectionTerm).^2)))"],
                       ", ")
            ContinueFlag = false
        end
        if ( (maximum(abs.(CorrectionTerm)) < 1e-10) || # convergenece conditions
             (sqrt(mean(CorrectionTerm-PrevCorrectionTerm).^2) <= 1e-10) )
            ContinueFlag = false
        end
        iteration <= 3 ? (ContinueFlag = true) : nothing # but makes sure there have been at least 3 iterations
        ContinueFlag == true ? PrevCorrectionTerm = copy(CorrectionTerm) : nothing # store this iteration's results unless we're done - in that case keep the last iteration's results for reporting.
        iteration += 1
    end
    if PrintoutStats == true
        println(join(["Finished extrinsic calibration:",
                      "$iteration iterations",
                      "max(|ε|)=$(map(y -> @sprintf("%.0e",y),maximum(abs.(CorrectionTerm))))",
                      "rms(ε₁-ε₀)=$(map(y -> @sprintf("%.0e",y),sqrt(mean(CorrectionTerm-PrevCorrectionTerm).^2)))"],
                     ", "))
    end
    return χ
end
###############################################################################

##############################################################################
""" helper function for HollandEtAlExtrinsicCalibration() """
function UpdateConstants(x, y, z, u, v, χ, ι)
    # ---------- distances from the GCP to the camera:
    Δx  = f_Δx(x,χ)
    Δy  = f_Δy(y,χ)
    Δz  = f_Δz(z,χ)
    # ---------- undistorted, scale corrected image coordinates of the GCP:
    uˢ  = f_uˢ(u,ι)
    vˢ  = f_vˢ(v,ι)
    # ---------- elements of the 3x3 orthogonal rotation matrix (aka direction cosines):
    m11  = f_m11(χ)
    m12  = f_m12(χ)
    m13  = f_m13(χ)
    # --
    m21  = f_m21(χ)
    m22  = f_m22(χ)
    m23  = f_m23(χ)
    # --
    m31  = f_m31(χ)
    m32  = f_m32(χ)
    m33  = f_m33(χ)
    # ---------- results using rotation matrix:
    o = m11 * Δx + m12 * Δy + m13 * Δz
    p = m21 * Δx + m22 * Δy + m23 * Δz
    q = m31 * Δx + m32 * Δy + m33 * Δz
    # ---------- bᵢⱼ are a set of constants used in LHS of the iterative calibration calculation:
    # Note that as calculated here they are multiplied by q.
    # b₁ⱼ are in the u coordinate direction
    qb11 = uˢ * (Δx*m32 -Δy*m31 ) + χ.f  * (Δx * m12 -Δy * m11 )
    qb12 = uˢ * (Δx * cos(χ.τ) * sin(χ.φ) + Δy * cos(χ.φ) * cos(χ.τ) + Δz * sin(χ.τ)) + χ.f  * (-Δx * sin(χ.φ) * sin(χ.τ) * sin(χ.σ) - Δy * cos(χ.φ) * sin(χ.τ) * sin(χ.σ) + Δz * cos(χ.τ) * sin(χ.σ))
    qb13 = χ.f * (Δx * m21 + Δy * m22 + Δz * m23 )
    qb14 = o
    qb15 = uˢ * m31 + χ.f * m11
    qb16 = uˢ * m32 + χ.f * m12
    qb17 = uˢ * m33 + χ.f * m13
    # ----------
    # b₂ⱼ are in the v coordinate direction
    qb21 = vˢ * (Δx * m32 -Δy * m31 ) + χ.f  * (Δx * m22 -Δy * m21 )
    qb22 = vˢ * (Δx * cos(χ.τ) * sin(χ.φ) + Δy * cos(χ.φ) * cos(χ.τ) + Δz * sin(χ.τ)) + χ.f  * (-Δx * sin(χ.φ) * sin(χ.τ) * cos(χ.σ)-Δy * cos(χ.φ) * sin(χ.τ) * cos(χ.σ) + Δz * cos(χ.τ) * cos(χ.σ))
    qb23 = -χ.f * (Δx * m11 + Δy * m12 + Δz * m13 ) #
    qb24 = p                   #
    qb25 = vˢ * m31 + χ.f * m21
    qb26 = vˢ * m32 + χ.f * m22
    qb27 = vˢ * m33 + χ.f * m23
    ###############################################################################
    # F₀, G₀ form the RHS of the system of equations:
    F₀ = uˢ + o/q * χ.f
    G₀ = vˢ + p/q * χ.f
    ###############################################################################
    Matrix_row1 = [qb11, qb12, qb13, qb14, -qb15, -qb16, -qb17] ./ q
    Matrix_row2 = [qb21, qb22, qb23, qb24, -qb25, -qb26, -qb27] ./ q
    return F₀, G₀, Matrix_row1, Matrix_row2
end
###############################################################################

###############################################################################
#
# Functions to calculate parameters defined in Holland et al. 1997:
#
# ---------- distances from the GCP to the camera:
f_Δx(x,χ) = x - χ.x_c                  # [meters]
f_Δy(y,χ) = y - χ.y_c                  # [meters]
f_Δz(z,χ) = z - χ.z_c                  # [meters]
# ---------- undistorted, scale corrected image coordinates of the GCP:
f_uˢ(u,ι) = (u - ι.u₀)ι.λᵤ               # [pixels]
f_vˢ(v,ι) = (v - ι.v₀)ι.λᵥ               # [pixels]
# ---------- elements of the 3x3 orthogonal rotation matrix (aka direction cosines):
f_m11(χ) =  cos(χ.φ)cos(χ.σ) + sin(χ.φ)cos(χ.τ)sin(χ.σ)  # [dimensionless]
f_m12(χ) = -sin(χ.φ)cos(χ.σ) + cos(χ.φ)cos(χ.τ)sin(χ.σ) # [dimensionless]
f_m13(χ) =  sin(χ.τ)sin(χ.σ)                       # [dimensionless]
# --
f_m21(χ) = -cos(χ.φ)sin(χ.σ) + sin(χ.φ)cos(χ.τ)cos(χ.σ) # [dimensionless]
f_m22(χ) =  sin(χ.φ)sin(χ.σ) + cos(χ.φ)cos(χ.τ)cos(χ.σ)  # [dimensionless]
f_m23(χ) =  sin(χ.τ)cos(χ.σ)                       # [dimensionless]
# --
f_m31(χ) =  sin(χ.φ)sin(χ.τ)                       # [dimensionless]
f_m32(χ) =  cos(χ.φ)sin(χ.τ)                       # [dimensionless]
f_m33(χ) = -cos(χ.τ)                            # [dimensionless]
# ---------- scaling coefficients from direction cosines and distances:
f_o(x,y,z,χ) = f_m11(χ)*f_Δx(x,χ) + f_m12(χ)*f_Δy(y,χ) + f_m13(χ)*f_Δz(z,χ) # [meters]
f_p(x,y,z,χ) = f_m21(χ)*f_Δx(x,χ) + f_m22(χ)*f_Δy(y,χ) + f_m23(χ)*f_Δz(z,χ) # [meters]
f_q(x,y,z,χ) = f_m31(χ)*f_Δx(x,χ) + f_m32(χ)*f_Δy(y,χ) + f_m33(χ)*f_Δz(z,χ) # [meters]
###############################################################################

###############################################################################
###############################################################################
###############################################################################
end                             # end of module
###############################################################################
###############################################################################
###############################################################################

max_stars_repo_name: saschweitzer/camera-calibration

max_issues_repo_name: saschweitzer/camera-calibration

max_forks_repo_name: saschweitzer/camera-calibration


text:
using Descartes
using GeometryBasics: Mesh

c2 = Cuboid([4,4.0,4])
cyl1 = rotate(pi/6, [1,0,0])translate([2,2,0])Cylinder(1.0,10.0)

u = Shell(0.5)RadiusedCSGUnion(1,c2, cyl1)
u2 = diff(u, Cuboid([2,2,2]))

m = Mesh(u2)

#Descartes.visualize(u2)

save("radiused_union.ply",m)


text:
# Create a model using Agents.jl

function initialize_model(param_file::AbstractString)
	# Load params
	d = YAML.load_file(param_file)

	r = Researcher(1, d["risk_taking_average"], d["generality_average"])
	d["time"] = 1

	model = ABM(typeof(r), properties=d, scheduler=Schedulers.randomly)

	researchers = start_population(d["population_size"], d["risk_taking_average"], d["generality_average"])
	for res in researchers
		add_agent!(res, model)
	end

	start_experienced_population(model, true)
	model.properties["max_grants"] = d["max_grants"] # restore this param
	model.properties["max_grants_org"] = deepcopy(d["max_grants"]) # restore this param

	return model
end

function initialize_model(;population_size, max_grants, years, risk_taking_average, generality_average, min_prof_exp, average_students, max_exp, mutation_rate, K, max_years, risk_range, A, B, C, save_gens, replicates)

	d = Dict(
		"risk_taking_average" => risk_taking_average,
		"generality_average" => generality_average,
		"population_size" => population_size,
		"max_grants" => max_grants,
		"years" => years,
		"min_prof_exp" => min_prof_exp,
		"average_students" => average_students,
		"max_exp" => max_exp,
		"mutation_rate" => mutation_rate,
		"K" => K,
		"max_years" => max_years,
		"risk_range" => risk_range,
		"A" => A,
		"B" => B,
		"C" => C,
		"save_gens" => save_gens,
		"replicates" => replicates
	)

	r = Researcher(1, d["risk_taking_average"], d["generality_average"])
	d["time"] = 1

	model = ABM(typeof(r), properties=d, scheduler=Schedulers.randomly)

	researchers = start_population(d["population_size"], d["risk_taking_average"], d["generality_average"])
	for res in researchers
		add_agent!(res, model)
	end

	start_experienced_population(model, true)
	model.properties["max_grants"] = d["max_grants"] # restore this param
	model.properties["max_grants_org"] = deepcopy(d["max_grants"]) # restore this param

	return model
end

function model_step!(model::ABM)
	retire!(model)
	exclude_unproductive!(model)
	get_grant!(model)
	# gain experience, do one year of research, publish finished work, receive citations
	update_researcher!(model)
	train_students!(model)
	# update time
	model.properties["time"] += 1
	# update the available grants per year
	model.properties["max_grants"] = deepcopy(model.properties["max_grants_org"])
end

"""
Create a homogeneous population of researchers.

parameters:
population_size: number of researchers.
number_of_grants: how many years of grants exists.
max_years: maximum number of years a researcher can receive a grant.
"""
function start_population(population_size, risk_taking::Float64=0.5, generality::Float64=0.5)
	researchers = [Researcher(i, risk_taking, generality) for i in 1:population_size]
	return researchers
end

"""
Create a heterogenous population of experienced researchers.
"""
function start_experienced_population(model, randomage=true)
	experience = model.properties["min_prof_exp"]
	population_size = model.properties["population_size"]

	if !randomage
		for year in 1:experience
			get_grant!(model)
			for resid in allids(model)
				update_researcher!(model[resid], K)  # gain experience, do one year of research, publish finished work, receive citations
			end
		end
	else
		age_groups = Int(round(population_size/experience))
		allind = collect(allids(model))
		ress = sample(allind, age_groups, replace=false, ordered=true)
		for ind in reverse(ress)
			splice!(allind, findfirst(a-> a==ind, allind)[1])
		end

		for year in 1:experience
			get_grant!(model)
			for resid in ress
				update_researcher!(model[resid], model)  # gain experience, do one year of research, publish finished work, receive citations
			end
			# remove the ress indices from allind and create new ress indices
			if age_groups >= length(allind)
				ress2 = allind
				ress = vcat(ress, allind)
				continue
			else
				ress2 = sample(allind, age_groups, replace=false, ordered=true)
				ress = vcat(ress, ress2)
			end
			for ind in reverse(ress2)
				splice!(allind, findfirst(a-> a==ind, allind))
			end
		end
	end
end


text:
using Distributions
using FileIO
using GeometryBasics
using ImageMagick
using Images
using OffsetArrays
using Plots
using Random
using RecipesPipeline
using RDatasets
using SparseArrays
using StaticArrays
using Statistics
using StatsPlots
using Test
using TestImages

@testset "unzip" begin
    x, y, z = unzip([(1., 2., 3.), (1., 2., 3.)])
    @test all(x .== 1.) && all(y .== 2.) && all(z .== 3.)
    x, y, z = unzip(Tuple{Float64, Float64, Float64}[])
    @test isempty(x) && isempty(y) && isempty(z)
end

@testset "group" begin
    include("test_group.jl")
end

@testset "plots" begin
    for i in eachindex(Plots._examples)
        if i ∉ Plots._backend_skips[:gr]
            @test Plots.test_examples(:gr, i, disp=false) isa Plots.Plot
        end
    end
end

#=
makie_test_dir = joinpath(@__DIR__, "test_makie")
mkpath(makie_test_dir)

@testset "Makie integration" begin
    cd(makie_test_dir) do
        include("makie.jl")
    end
end
=#


text:
using StatsBase

abstract type Species end
abstract type Agent{S<:Species} end

abstract type PlantSpecies <: Species end
abstract type Grass <: PlantSpecies end

abstract type AnimalSpecies <: Species end
abstract type Sheep <: AnimalSpecies end
abstract type Wolf <: AnimalSpecies end

abstract type Sex end
abstract type Male <: Sex end
abstract type Female <: Sex end

id(a::Agent) = a.id

##########  World Definition  ##################################################

mutable struct World{A<:Agent}
    agents::Dict{Int,A}
    max_id::Int
end

function World(agents::Vector{<:Agent})
    ids = id.(agents)
    length(unique(ids)) == length(agents) || error("Not all agents have unique IDs!")
    World(Dict(id(a)=>a for a in agents), maximum(ids))
end

function world_step!(world::World)
    for id in deepcopy(keys(world.agents))
        !haskey(world.agents,id) && continue
        a = world.agents[id]
        agent_step!(a,world)
    end
end

function Base.show(io::IO, w::World)
    println(io, typeof(w))
    for (_,a) in w.agents
        println(io,"  $a")
    end
end

##########  Plant Definition  ##################################################

mutable struct Plant{P<:PlantSpecies} <: Agent{P}
    id::Int
    size::Int
    max_size::Int
end

Base.size(a::Plant) = a.size
max_size(a::Plant) = a.max_size
grow!(a::Plant) = a.size += 1

# constructor for all Plant{<:PlantSpecies} callable as PlantSpecies(...)
(A::Type{<:PlantSpecies})(id, s, m) = Plant{A}(id,s,m)
(A::Type{<:PlantSpecies})(id, m) = (A::Type{<:PlantSpecies})(id,rand(1:m),m)

function agent_step!(a::Plant, w::World)
    if size(a) != max_size(a)
        grow!(a)
    end
end

function Base.show(io::IO, p::Plant{P}) where P
    x = size(p)/max_size(p) * 100
    print(io,"$P  #$(id(p)) $(round(Int,x))% grown")
end

Base.show(io::IO, ::Type{Grass}) = print(io,"🌿")

##########  Animal Definition  #################################################

mutable struct Animal{A<:AnimalSpecies,S<:Sex} <: Agent{A}
    id::Int
    energy::Float64
    Δenergy::Float64
    reprprob::Float64
    foodprob::Float64
end

energy(a::Animal) = a.energy
Δenergy(a::Animal) = a.Δenergy
reprprob(a::Animal) = a.reprprob
foodprob(a::Animal) = a.foodprob
energy!(a::Animal, e) = a.energy = e
incr_energy!(a::Animal, Δe) = energy!(a, energy(a)+Δe)

function (A::Type{<:AnimalSpecies})(id::Int, E, ΔE, pr, pf, S=rand(Bool) ? Female : Male)
    Animal{A,S}(id,E,ΔE,pr,pf)
end

function agent_step!(a::Animal, w::World)
    incr_energy!(a,-1)
    if rand() <= foodprob(a)
        dinner = find_food(a,w)
        eat!(a, dinner, w)
    end
    if energy(a) <= 0
        kill_agent!(a,w)
        return
    end
    if rand() <= reprprob(a)
        reproduce!(a,w)
    end
    return a
end

function find_rand(f, w::World)
    xs = filter(f, w.agents |> values |> collect)
    isempty(xs) ? nothing : sample(xs)
end

find_food(a::Animal, w::World) = find_rand(x->eats(a,x),w)

eats(::Animal{Sheep},p::Plant{Grass}) = size(p)>0
eats(::Animal{Wolf},::Animal{Sheep}) = true
eats(::Agent,::Agent) = false

function eat!(a::Animal{Wolf}, b::Animal{Sheep}, w::World)
    incr_energy!(a, energy(b)*Δenergy(a))
    kill_agent!(b,w)
end
function eat!(a::Animal{Sheep}, b::Plant{Grass}, w::World)
    incr_energy!(a, size(b)*Δenergy(a))
    b.size = 0
end
eat!(::Animal,::Nothing,::World) = nothing

function reproduce!(a::A, w::World) where A<:Animal
    b = find_mate(a,w)
    if !isnothing(b)
        energy!(a, energy(a)/2)
        a_vals = [getproperty(a,n) for n in fieldnames(A) if n!=:id]
        new_id = w.max_id + 1
        â = A(new_id, a_vals...)
        w.agents[id(â)] = â
        w.max_id = new_id
    end
end

find_mate(a::Animal, w::World) = find_rand(x->mates(a,x),w)

function mates(a,b)
    error("""You have to specify the mating behaviour of your agents by overloading `mates` e.g. like this:

        mates(a::Animal{S,Female}, b::Animal{S,Male}) where S<:Species = true
        mates(a::Animal{S,Male}, b::Animal{S,Female}) where S<:Species = true
        mates(a::Agent, b::Agent) = false
    """)
end

kill_agent!(a::Animal, w::World) = delete!(w.agents, id(a))

Base.show(io::IO, ::Type{Sheep}) = print(io,"🐑")
Base.show(io::IO, ::Type{Wolf}) = print(io,"🐺")
Base.show(io::IO, ::Type{Male}) = print(io,"♂")
Base.show(io::IO, ::Type{Female}) = print(io,"♀")
function Base.show(io::IO, a::Animal{A,S}) where {A,S}
    e = energy(a)
    d = Δenergy(a)
    pr = reprprob(a)
    pf = foodprob(a)
    print(io,"$A$S #$(id(a)) E=$e ΔE=$d pr=$pr pf=$pf")
end


text:
@testset "$TEST $G" begin

f, fio = mktemp()

p = readgraph(joinpath(testdir,"testdata","tutte.gml"), G)
g = graph(:tutte, G)
@test p == g

gml1 = readgraph(joinpath(testdir,"testdata", "twographs-10-28.gml"), G)
@test nv(gml1) == 10
@test ne(gml1) == 28
writegraph(f, gml1, :gml)
gml1a = readgraph(f, :gml, G)
@test gml1a == gml1

gml1 = readgraph(joinpath(testdir,"testdata", "twounnamedgraphs.gml"), G)
@test typeof(gml1) == G
@test nv(gml1) == 4
@test ne(gml1) == 6
writegraph(f, gml1, :gml)
gml1a = readgraph(f, :gml, G)
@test gml1a == gml1

g = G(10,0)
writegraph(f, g, :gml)
ga = readgraph(f, :gml, G)
@test g == ga

g = G(10,20)
writegraph(f, g, :gml)
ga = readgraph(f, :gml, G)
@test g == ga

g = DG(10,20)
writegraph(f, g, :gml)
ga = readgraph(f, :gml, G)
@test g == ga

if G <: ANetwork ##################

g = readnetwork(joinpath(testdir, "testdata", "testnet.gml"), G)
writenetwork(f, g, :gml)
h = readnetwork(f, :gml, G)
test_networks_eq(g, h)

for gname in [:lesmis, :karate]
    g = readnetwork(gname, G)
    rem_gprop!(g, "readme")
    rem_vprop!(g, "pos") # TODO vector attributes not supported
    writenetwork(f, g, :gml)
    h = readnetwork(f, :gml, G)
    test_networks_eq(g, h)
end

end #if network ####################

isfile(f) && rm(f)
 end #testset


text:
struct SSHKey
    id::Integer
    fingerprint::String
    public_key::String
    name::String

    function SSHKey(data::Dict{String})
        new(
            data["id"],
            data["fingerprint"],
            data["public_key"],
            data["name"]
        )
    end
end

function show(io::IO, s::SSHKey)
    print(io, "SSHKey ($(s.name))")
end

function getallsshkeys!(client::AbstractClient)
    response = getdata!(client, joinpath(ENDPOINT, "account", "keys"))

    if floor(response.status/100) == 2 # OK
        body = parse(String(response.body))
        meta = body["meta"]
        links = body["links"]
        data = body["ssh_keys"]

        ssh_keys = Array{SSHKey, 1}(UndefInitializer(), meta["total"])

        for (i, ssh_key) in enumerate(data)
            ssh_keys[i] = SSHKey(ssh_key)
        end
    else
        error("Received error $(response.status)")
    end

    ssh_keys
end

function getsshkey!(client::AbstractClient, key_id::Union{Integer, String})
    response = getdata!(client, joinpath(ENDPOINT, "account", "keys", "$key_id"))

    if floor(response.status/100) == 2 # OK
        body = parse(String(response.body))
        data = body["ssh_key"]

        ssh_key = SSHKey(data)
    else
        error("Received error $(response.status)")
    end
end

function getsshkey!(client::AbstractClient, key::SSHKey)
    getsshkey!(client, key.id)
end

function createsshkey!(client::AbstractClient; kwargs...)
    body = Dict([String(k[1]) => k[2] for k in kwargs])

    if !haskey(body, "name")
        error("'name' is a required argument")
    end

    if !haskey(body, "public_key")
        error("'public_key' is a required argument")
    end

    response = postdata!(client, joinpath(ENDPOINT, "account", "keys"), body)

    if floor(response.status/100) == 2 # OK
        body = parse(String(response.body))
        data = body["ssh_key"]

        ssh_key = SSHKey(data)
    else
        error("Received error $(response.status)")
    end
end

function update_ssh_key(client::AbstractClient, key_id::Union{Integer, String}; kwargs...)
    body = Dict([String(k[1]) => k[2] for k in kwargs])

    if !haskey(body, "name")
        error("'name' is a required argument")
    end

    uri = joinpath(ENDPOINT, "account", "keys", "$(key_id)")
    response = putdata!(client, uri, body)

    if floor(response.status/100) == 2 # OK
        body = parse(String(response.body))
        data = body["ssh_key"]

        ssh_key = SSHKey(data)
    else
        error("Received error $(response.status)")
    end
end

function update_ssh_key(client::AbstractClient, key::SSHKey; kwargs...)
    update_ssh_key(client, key.id; kwargs...)
end

function deletesshkey!(client::AbstractClient, key_id::Union{Integer, String})
    deletedata!(client, joinpath(ENDPOINT, "account", "keys", "$(key_id)"))
end

function deletesshkey!(client::AbstractClient, key::SSHKey)
    deletesshkey!(client, key.id)
end


text:
module Multigrid
	using LinearAlgebra
	using SparseArrays
	using Statistics
	using ParSpMatVec;
	using jInv.LinearSolvers
	using jInv.Utils

	const global hasParSpMatVec = ParSpMatVec.isBuilt();
	if hasParSpMatVec == false
		println("hasParSpMatVec==false. ParSpMatVec has failed to build!!!!!!!");
	end

	import jInv.LinearSolvers.JuliaSolver, jInv.LinearSolvers.getJuliaSolver
	import jInv.LinearSolvers.AbstractSolver
	import jInv.LinearSolvers.solveLinearSystem!,jInv.LinearSolvers.solveLinearSystem, jInv.LinearSolvers.setupSolver
	import jInv.Utils.clear!

	SparseCSCTypes = Union{SparseMatrixCSC{ComplexF64,Int64},SparseMatrixCSC{Float64,Int64},SparseMatrixCSC{ComplexF32,Int64},SparseMatrixCSC{Float32,Int64}}
	ArrayTypes = Union{Array{ComplexF64},Array{ComplexF32},Array{Float64},Array{Float32}}

	export SparseCSCTypes,ArrayTypes
	include("ParallelJuliaSolver/parallelJuliaSolver.jl");
	include("DomainDecomposition/DomainDecomposition.jl");
	include("Multigrid/MGdef.jl");
end

# check if MUMPS can be used
# const minMUMPSversion = VersionNumber(0,0,1)
# hasMUMPS=false
# vMUMPS = VersionNumber(0,0,0)
# try
	# vMUMPS = Pkg.installed("MUMPS")
	# hasMUMPS = vMUMPS >= minMUMPSversion
	# if hasMUMPS
		# using MUMPS;
	# end
# catch
# end

# end # module


text:
using KernelLibrary, KernelRecursiveLeastSquares, Winston

x = randn(5000,20);
h = randn(20,1);
h2 = randn(20,1);
y = (x*h + 1).^3 + 4*(x*h2 + 1).^2+ 500*randn(5000,1);

x2 = randn(200,20);
y2 = (x2*h + 1).^3 + 4*(x2*h2 + 1).^2;

x=x'; y=y';
x2=x2'; y2=y2';

k(x,y) = polynomial_kernel(x,y,2);

(alpha1, dict1, Kinv1, dict1_idx) = krls(x, y, kernelfunc=k, nu=1., maxdict=30);

@time (alpha1, dict1, Kinv1, dict1_idx) = krls(x, y, kernelfunc=k, nu=1., maxdict=300);

ypred1 = zeros(size(x2,2),1);
for ii in [1:size(x2,2)]
	tmp = k(dict1, x2[:,ii])*alpha1;
	ypred1[ii]=tmp[1];
end

cor(ypred1, y2')

k(x,y) = multiquad_kernel(x,y,1);

@time (alpha2, dict2, Kinv2, dict2_idx) = krls(x, y, kernelfunc=k, nu=3., maxdict=300, index="rand");

ypred2 = zeros(size(x2,2),1);
for ii in [1:size(x2,2)]
	tmp = k(dict2, x2[:,ii])*alpha2;
	ypred2[ii]=tmp[1];
end

cor(ypred2, y2')

k(x,y) = power_kernel(x,y,3);

@time (alpha3, dict3, Kinv3, dict3_idx) = krls(x, y, kernelfunc=k, nu=.5, maxdict=300, index="rand");

ypred3 = zeros(size(x2,2),1);
for ii in [1:size(x2,2)]
	tmp = k(dict3, x2[:,ii])*alpha3;
	ypred3[ii]=tmp[1];
end

cor(ypred3, y2')

hold(false)
plot(y2')
hold(true)
plot(ypred1, "r");
plot(ypred2, "g")
plot(ypred3, "b")

mp = mean([ypred1 ypred2 ypred3],2);
plot(mp,"m")
cor(mp, y2')

max_stars_repo_name: the-moliver/KernelRecursiveLeastSquares

max_issues_repo_name: the-moliver/KernelRecursiveLeastSquares

max_forks_repo_name: the-moliver/KernelRecursiveLeastSquares


text:
# change this so it works on windows
function download_with_curl(  url, dest)
    println("Running curl_cmd")
    if !isfile(dest)
        run(`curl -o $dest -L $url`)
        println("Downloaded $url to $dest")
    else
        println("$dest exists no need to download")
    end
end

function list_annotations()
 # UCSC
 # GENCODE
 # ENSEMBL
end


text:
#=
    Ninjas, Bears, hunters

    Port of the PSI model
    https://github.com/eth-sri/psi/blob/master/test/ninjasBearsHunters.psi

    Similar to Rock, Paper, Scissors but there are some differences:
    - At start there are certain amount of ninjas, bears, and hunters. It
      might be different for each type.
    - In each turn, two of separare types are selected and
      one of them is eliminated according to the elimination rules
      (see below).
    - We continue until there is just one left.
    - The elimination rules are
        * bear kills ninja
        * hunter kills bear
        * ninja kills hunter

   What are the probabilities that the winner is a specific type?

   The PSI model starts with
      10 ninjas, 2 bears, and 5 hunters
   and the the exact probabilities of survival are:
    (probNinjaSurvive,probBearsSurvive,probHuntersSurvive)=
    (0.221480658306459,0.767229372971558,0.0112899687219824)

   It is - at first sight - a little surprising that the probability
   of the ninja as the winner is so low (0.22). But since there are so many
   ninjas, the hunters are killed quite soon which mean that there is no hunter
   left that can kill the bears (which kills the ninjas). Thus the bears have the greatest
   probability of survival (0.77).

   This Turing model:
   Summary Statistics
            parameters      mean       std   naive_se      mcse         ess       rhat   ess_per_sec
               Symbol   Float64   Float64    Float64   Float64    Float64?   Float64?       Float64

      probNinjaSurvive    0.2240    0.4171     0.0132    0.0167    951.6074     0.9992     5228.6121
      probBearsSurvive    0.7590    0.4279     0.0135    0.0151    963.7017     0.9990     5295.0643
    probHuntersSurvive    0.0170    0.1293     0.0041    0.0047   1037.0165     1.0005     5697.8929

   Quantiles
            parameters      2.5%     25.0%     50.0%     75.0%     97.5%
               Symbol   Float64   Float64   Float64   Float64   Float64

      probNinjaSurvive    0.0000    0.0000    0.0000    0.0000    1.0000
      probBearsSurvive    0.0000    1.0000    1.0000    1.0000    1.0000
    probHuntersSurvive    0.0000    0.0000    0.0000    0.0000    0.0000

    Some other examples:
    * start with 3 ninjas, 4 bears, and 3 hunters.
      The probability of survival is
      probNinjaSurvive  : 0.2520
      probBearsSurvive  : 0.2880
      probHuntersSurvive: 0.4600
      (Exact probability from PSI:
       3927004/15690675,2846092/9984975,3402579/7322315 = 0.250276294678209,0.285037468796867,0.464686236524924
       )

    * start with 1 ninjas, 2 bears, and 3 hunters.
      The probability of survival is
      probNinjaSurvive  : 0.4420
      probBearsSurvive  : 0.0810
      probHuntersSurvive: 0.4770

      (Exact probability from PSI:
       339/770,111/1540,751/1540 = 0.44025974025974,0.0720779220779221,0.487662337662338
      )

    Remember;
    * bear kills ninja
    * hunter kills bear
    * ninja kills hunter

    Also see
    * https://en.wikipedia.org/wiki/Rock_paper_scissors
    * It seems to have been popularized by the FedEx commercial: Bear, hunter ninja.
      https://www.youtube.com/watch?v=AhFbbq0zpHY
    * (The PSI model refers to a link that's not accessible.)

=#
using Turing, StatsPlots, DataFrames, Distributions, Random
include("jl_utils.jl")

@model function ninjas_bears_hunters(ninjas=10, bears=2, hunters=5)
    n = ninjas + bears + hunters
    NINJA = 1
    BEAR = 2
    HUNTER = 3
    p1 = tzeros(1)
    p2 = tzeros(1)
    i = 0
    while true
      i += 1

      s = ninjas + bears + hunters
      v = [ninjas,bears,hunters]
      # println("i:$i v:$v")
      if s <= 1 || count(t -> t > 0,v) == 1
         break
      end

      a = [ninjas/s,bears/s,hunters/s]

      push!(p1,0)
      push!(p2,0)
      p1[i] ~ Categorical(simplex(a))
      p2[i] ~ Categorical(simplex(a))
      # true ~ Dirac(p1[i] != p2[i])

      # there must be at least two remaining civilizations.
      if p1[i] != p2[i] && s > 1 && ninjas + bears > 0 && ninjas + hunters > 0 && bears + hunters > 0

         # true ~ Dirac(!(ninjas  == 0 && (p1[i] == NINJA  || p2[i] == NINJA)))
         # true ~ Dirac(!(bears	  == 0 && (p1[i] == BEAR   || p2[i] == BEAR)))
         # true ~ Dirac(!(hunters == 0 && (p1[i] == HUNTER || p2[i] == HUNTER)))

         if		((p1[i] == NINJA && p2[i] == BEAR)   || (p1[i] == BEAR   && p2[i] == NINJA) && ninjas  > 0) ninjas  -=1 # Bear kills ninja
         elseif ((p1[i] == NINJA && p2[i] == HUNTER) || (p1[i] == HUNTER && p2[i] == NINJA) && hunters > 0) hunters -=1 # Ninja kills hunter
         elseif ((p1[i] == BEAR  && p2[i] == HUNTER) || (p1[i] == HUNTER && p2[i] == BEAR)  && bears   > 0) bears   -=1 # Hunter kills bear
         end
      end
   end

   # true ~ Dirac(!(ninjas == 0 && bears == 0 && hunters == 0))
   # true ~ Dirac(ninjas + bears == 0 || ninjas + hunters == 0 || bears + hunters == 0)
   probNinjaSurvive   ~ Dirac(ninjas > 0)
   probBearsSurvive   ~ Dirac(bears > 0)
   probHuntersSurvive ~ Dirac(hunters > 0)
   return
end

num_ninjas=10
num_bears=2
num_hunters=5
# num_ninjas=3
# num_bears=4
# num_hunters=3

# num_ninjas=1
# num_bears=2
# num_hunters=3
model = ninjas_bears_hunters(num_ninjas,num_bears,num_hunters)

# Since we don't do any observations, we can use Prior()
chns = sample(model, Prior(), 1_000)
# chns = sample(model, MH(), 1_000)
# chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 1_000)
# chns = sample(model, IS(), 1_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns[[:probNinjaSurvive,:probBearsSurvive,:probHuntersSurvive]])

# show_var_dist_pct(chns,:probNinjaSurvive)
# show_var_dist_pct(chns,:probBearsSurvive)
# show_var_dist_pct(chns,:probHuntersSurvive)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
mutable struct CoefplotOption
    xticklabel_angle::Float64
    xticklabel_size::Float64
    yticklabel_angle::Float64
    yticklabel_size::Float64
    other_coefplot_options::PGFPlotsX.Options
end

mutable struct SinglecoefplotOption
    dotsize::Float64
    dotcolor::Color
    linewidth::Float64
    linecolor::Color
    other_dot_options::PGFPlotsX.Options
    other_line_options::PGFPlotsX.Options
end

default_singlecoefplot_options() = SinglecoefplotOption(4, first(COLOR_PALATTE), 0.5, first(COLOR_PALATTE), PGFPlotsX.Options(), PGFPlotsX.Options())

default_coefplot_options() = CoefplotOption(0, 5, 0, 5, PGFPlotsX.Options())

get_line_options(option::SinglecoefplotOption) = merge(option.other_line_options,
                                                       color_as_option(option.linecolor),
                                                       PGFPlotsX.Options(Symbol("line width") => option.linewidth))

get_dot_options(option::SinglecoefplotOption) = merge(option.other_dot_options,
                                                      color_as_fill_option(option.dotcolor),
                                                      PGFPlotsX.Options(:circle => nothing, # TODO consider make this flexible
                                                                        Symbol("inner sep") => 0,
                                                                        Symbol("minimum size") => option.dotsize))

get_coefplot_options(option::CoefplotOption) = merge(option.other_coefplot_options,
                                                     PGFPlotsX.Options(Symbol("xticklabel style")=> "font=\\fontsize{$(option.xticklabel_size)}{$(option.xticklabel_size)}\\selectfont, rotate=$(option.xticklabel_angle)"),
                                                     PGFPlotsX.Options(Symbol("yticklabel style")=> "font=\\fontsize{$(option.yticklabel_size)}{$(option.yticklabel_size)}\\selectfont, rotate=$(option.yticklabel_angle)"))

Base.:(==)(o1::PGFPlotsX.Options,o2::PGFPlotsX.Options) = (o1.dict == o2.dict) && (o1.print_empty == o2.print_empty)
function Base.:(==)(so1::SinglecoefplotOption,so2::SinglecoefplotOption)
    for fieldname in fieldnames(typeof(so1))
        if getfield(so1,fieldname) == getfield(so2,fieldname)
            continue
        else
            return false
        end
    end
    return true
end
function all_equal(v::Vector{SinglecoefplotOption})
    for i in 1:(length(v)-1)
        if v[i] == v[i+1]
            continue
        else
            return false
        end
    end
    return true
end

mutable struct LegendsOption
    fontsize::Real
    at::Tuple{Real, Real}
    anchor::Symbol

    function LegendsOption(fontsize::Real=5,at::Tuple{Real, Real}=(0.98,0.98),anchor::Symbol=Symbol("north east"))
        new(fontsize,at,anchor)
    end
end

function gen_legend_options(lo::LegendsOption)
    legend_stype = PGFPlotsX.Options(:font => "\\fontsize{$(lo.fontsize)}{$(lo.fontsize)}\\selectfont",
                                     :at => "($(first(lo.at)),$(last(lo.at)))",
                                     :anchor => string(lo.anchor))
    PGFPlotsX.Options(Symbol("legend style") => legend_stype)
end


text:
"""
The barrel vault (Scordelis-Lo) roof is one of the benchmarks for linear elastic
analysis of shells.

The candidate element's usefulness in irregular geometries (and most practical
cases involve a high degree of geometric irregularity) is tested. As would be
expected,the irregular mesh results are not as good as those provided by a
regular meshwith the same number of variables.

Problem description

The physical basis of the problem is a deeply arched roof supported only
bydiaphragms at its curved edges (an aircraft hanger), deforming under its own
weight. It is interesting to observe that the geometry is such that the
centerpoint of the roof moves upward under the self-weight(downwardly directed)
load. Perhaps this is one reason why the problem is not straightforward
numerically.
"""
module barrel_vault_examples

using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using FinEtools.MeshExportModule.VTKWrite: vtkwrite

function _execute(input = "barrelvault_s3r_fineirreg.inp", visualize = true)
    E = 3.0e6;
    nu = 0.0;
    thickness = 3.0;
    tolerance = thickness/20
    # analytical solution for the vertical deflection under the load
    analyt_sol = -0.3024*12;
    R = 25.0*12;
    L = 50.0*12;
    formul = FEMMShellT3FFModule

    output = import_ABAQUS(joinpath(dirname(@__FILE__()), input))
    fens = output["fens"]
    fes = output["fesets"][1]

     # boundingbox(fens.xyz)
    connected = findunconnnodes(fens, fes);
    fens, new_numbering = compactnodes(fens, connected);
    fes = renumberconn!(fes, new_numbering);

    fens, fes = mergenodes(fens, fes, thickness/10)

    # plots = cat(plot_space_box([[0 0 -R/2]; [R/2 R/2 R/2]]),
    #     plot_nodes(fens),
    #     plot_midsurface(fens, fes);
    # dims = 1)
    # pl = render(plots)

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

    @info "Mesh: $input"

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    femm.drilling_stiffness_scale = 0.1
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # vtkwrite("barrel_vault_examples-mesh.vtu", fens, fes)
    # bfes = meshboundary(fes)
    # vtkwrite("barrel_vault_examples-boundary-mesh.vtu", fens, bfes)

    # Apply EBC's
    # Rigid diaphragm end
    l1 = selectnode(fens; box = Float64[-Inf Inf -Inf Inf 0 0], inflate = tolerance)
    for i in [1,2]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to Z
    l1 = selectnode(fens; box = Float64[-Inf Inf -Inf Inf L/2 L/2], inflate = tolerance)
    for i in [3,4,5]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to Y (apex)
    l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,4,6]
        setebc!(dchi, l1, true, i)
    end

    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Load
    # Midpoint of the free edge
    nl = selectnode(fens; box = Float64[-Inf Inf sin(40/360*2*pi)*R sin(40/360*2*pi)*R L/2 L/2], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))
    fi = ForceIntensity(FFlt[-0.625, 0, 0, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);

    # @infiltrate
    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])

    targetu =  dchi.values[nl, 1][1]
    @info "Solution: $(round(targetu, digits=8)),  $(round(targetu/analyt_sol, digits = 4)*100)%"

    # Generate a graphical display of resultants
    cylindrical!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) = begin
        r = -vec(XYZ); r[3] = 0.0
        csmatout[:, 3] .= vec(r)/norm(vec(r))
        csmatout[:, 2] .= (0.0, 0.0, 1.0)
        cross3!(view(csmatout, :, 1), view(csmatout, :, 2), view(csmatout, :, 3))
        return csmatout
    end
    ocsys = CSys(3, 3, cylindrical!)
    scalars = []
    for nc in 1:3
        fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("m$nc", fld.values))
    end
    vtkwrite("$(input)-m.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
    scalars = []
    for nc in 1:3
        fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("n$nc", fld.values))
    end
    vtkwrite("$(input)-n.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
    scalars = []
    for nc in 1:2
        fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("q$nc", fld.values))
    end
    vtkwrite("$(input)-q.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])

    # Visualization
    if !visualize
        return true
    end
    scattersysvec!(dchi, (R/2)/maximum(abs.(U)).*U)
    update_rotation_field!(Rfield0, dchi)
    plots = cat(plot_space_box([[0 0 -R]; [R R R]]),
        plot_nodes(fens),
        plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
    dims = 1)
    pl = render(plots)
    return true
end

function test_convergence()

    @info "Scordelis-Lo Abaqus model"

    _execute("barrelvault_stri3_irreg.inp", false)

    _execute("barrelvault_s3r_fineirreg.inp", false)

    return true
end

function allrun()
    println("#####################################################")
    println("# test_convergence ")
    test_convergence()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
import Base: *

export sortVarNested, sortDFG
export isPrior, lsfPriors
export getVariableType, getSofttype
export getFactorType, getfnctype
export lsTypes, lsfTypes
export lsWho, lsfWho
export *

*(a::Symbol, b::AbstractString)::Symbol = Symbol(string(a,b))

## Utility functions for getting type names and modules (from IncrementalInference)
function _getmodule(t::T) where T
  T.name.module
end
function _getname(t::T) where T
  T.name.name
end

"""
    $(SIGNATURES)
Test if all elements of the string is a number:  Ex, "123" is true, "1_2" is false.
"""
allnums(str::S) where {S <: AbstractString} = occursin(Regex(string(["[0-9]" for j in 1:length(str)]...)), str)
# occursin(r"_+|,+|-+", node_idx)

isnestednum(str::S; delim='_') where {S <: AbstractString} = occursin(Regex("[0-9]+$(delim)[0-9]+"), str)

function sortnestedperm(strs::Vector{<:AbstractString}; delim='_')
  str12 = split.(strs, delim)
  sp1 = sortperm(parse.(Int,getindex.(str12,2)))
  sp2 = sortperm(parse.(Int,getindex.(str12,1)[sp1]))
  return sp1[sp2]
end

function getFirstNumericalOffset(st::AS) where AS <: AbstractString
  i = 1
  while !allnums(st[i:i])  i+=1; end
  return i
end

"""
    $SIGNATURES

Sort a variable list which may have nested structure such as `:x1_2` -- does not sort for alphabetic characters.
"""
function sortVarNested(vars::Vector{Symbol})::Vector{Symbol}
	# whos nested and first numeric character offset
	sv = string.(vars)
	offsets = getFirstNumericalOffset.(sv)
	masknested = isnestednum.(sv)
	masknotnested = true .⊻ masknested

	# strip alphabetic characters from front
	msv = sv[masknotnested]
	msvO = offsets[masknotnested]
	nsv = sv[masknested]
	nsvO = offsets[masknested]

	# do nonnested list separately
	nnreducelist = map((s,o) -> s[o:end], msv, msvO)
	nnintlist = parse.(Int, nnreducelist)
	nnp = sortperm(nnintlist)
	nnNums = nnintlist[nnp] # used in mixing later
	nonnested = msv[nnp]
	smsv = vars[masknotnested][nnp]

	# do nested list separately
	nestedreducelist = map((s,o) -> s[o:end], nsv, nsvO)
	nestedp = sortnestedperm(nestedreducelist)
	nesNums = parse.(Int, getindex.(split.(nestedreducelist[nestedp], '_'),1)) # used in mixing later
	nested = nsv[nestedp]
	snsv = vars[masknested][nestedp]

	# mix back together, pick next sorted item from either pile
	retvars = Vector{Symbol}(undef, length(vars))
	nni = 1
	nesi = 1
	lsmsv = length(smsv)
	lsnsv = length(snsv)
	MAXMAX = 999999999999
	for i in 1:length(vars)
	  # inner ifs to ensure bounds and correct sorting at end of each list
	  if (nni<=lsmsv ? nnNums[nni] : MAXMAX) <= (nesi<=lsnsv ? nesNums[nesi] : MAXMAX)
	    retvars[i] = smsv[nni]
		nni += 1
	  else
		retvars[i] = snsv[nesi]
		nesi += 1
	  end
	end
	return retvars
end

"""
    $SIGNATURES

Sort variable (factor) lists in a meaningful way, for example `[:april;:x1_3;:x1_6;]`

Notes
- Not fool proof, but does better than native sort.

Example

`sortDFG(ls(dfg))`

Related

ls, lsf
"""
sortDFG(vars::Vector{Symbol})::Vector{Symbol} = sortVarNested(vars)

"""
    $SIGNATURES

Return the factor type used in a `::DFGFactor`.

Notes:
- OBSOLETE, use newer getFactorType instead.

Related

getFactorType
"""
function getfnctype(data::GenericFunctionNodeData)
  # TODO what is this?
  if typeof(data).name.name == :VariableNodeData
    return VariableNodeData
  end

  # this looks right
  return data.fnc.usrfnc!
end
function getfnctype(fact::DFGFactor; solveKey::Symbol=:default)
  data = solverData(fact) # TODO , solveKey=solveKey)
  return getfnctype(data)
end
function getfnctype(dfg::T, lbl::Symbol; solveKey::Symbol=:default) where T <: AbstractDFG
  getfnctype(getFactor(dfg, exvertid))
end

"""
    $SIGNATURES

Return user factor type from factor graph identified by label `::Symbol`.

Notes
- Replaces older `getfnctype`.
"""
getFactorType(data::GenericFunctionNodeData) = data.fnc.usrfnc!
getFactorType(fct::DFGFactor) = getFactorType(solverData(fct))
function getFactorType(dfg::G, lbl::Symbol) where G <: AbstractDFG
  getFactorType(getFactor(dfg, lbl))
end

"""
    $SIGNATURES

Return `::Bool` on whether given factor `fc::Symbol` is a prior in factor graph `dfg`.
"""
function isPrior(dfg::G, fc::Symbol)::Bool where G <: AbstractDFG
  fco = getFactor(dfg, fc)
  getfnctype(fco) isa FunctorSingleton
end

"""
    $SIGNATURES

Return vector of prior factor symbol labels in factor graph `dfg`.
"""
function lsfPriors(dfg::G)::Vector{Symbol} where G <: AbstractDFG
  priors = Symbol[]
  fcts = lsf(dfg)
  for fc in fcts
    if isPrior(dfg, fc)
      push!(priors, fc)
    end
  end
  return priors
end

"""
   $(SIGNATURES)

Variable nodes softtype information holding a variety of meta data associated with the type of variable stored in that node of the factor graph.

Related

getVariableType
"""
function getSofttype(vnd::VariableNodeData)
  return vnd.softtype
end
function getSofttype(v::DFGVariable; solveKey::Symbol=:default)
  return getSofttype(solverData(v, solveKey))
end

"""
    $SIGNATURES

Return the DFGVariable softtype in factor graph `dfg<:AbstractDFG` and label `::Symbol`.
"""
getVariableType(var::DFGVariable; solveKey::Symbol=:default) = getSofttype(var, solveKey=solveKey)
function getVariableType(dfg::G, lbl::Symbol; solveKey::Symbol=:default) where G <: AbstractDFG
  getVariableType(getVariable(dfg, lbl), solveKey=solveKey)
end

"""
    $SIGNATURES

Return `::Dict{Symbol, Vector{String}}` of all unique factor types in factor graph.
"""
function lsfTypes(dfg::G)::Dict{Symbol, Vector{String}} where G <: AbstractDFG
  alltypes = Dict{Symbol,Vector{String}}()
  for fc in lsf(dfg)
    Tt = typeof(getFactorType(dfg, fc))
    sTt = string(Tt)
    name = Symbol(Tt.name)
    if !haskey(alltypes, name)
      alltypes[name] = String[string(Tt)]
    else
      if sum(alltypes[name] .== sTt) == 0
        push!(alltypes[name], sTt)
      end
    end
  end
  return alltypes
end

"""
    $SIGNATURES

Return `::Dict{Symbol, Vector{String}}` of all unique variable types in factor graph.
"""
function lsTypes(dfg::G)::Dict{Symbol, Vector{String}} where G <: AbstractDFG
  alltypes = Dict{Symbol,Vector{String}}()
  for fc in ls(dfg)
    Tt = typeof(getVariableType(dfg, fc))
    sTt = string(Tt)
    name = Symbol(Tt.name)
    if !haskey(alltypes, name)
      alltypes[name] = String[string(Tt)]
    else
      if sum(alltypes[name] .== sTt) == 0
        push!(alltypes[name], sTt)
      end
    end
  end
  return alltypes
end

function ls(dfg::G, ::Type{T}; solveKey::Symbol=:default) where {G <: AbstractDFG, T <: InferenceVariable}
  xx = getVariables(dfg)
  mask = getVariableType.(xx, solveKey=solveKey) .|> typeof .== T
  vxx = view(xx, mask)
  map(x->x.label, vxx)
end

function ls(dfg::G, ::Type{T}) where {G <: AbstractDFG, T <: FunctorInferenceType}
  xx = getFactors(dfg)
  names = getfield.(typeof.(getFactorType.(xx)), :name) .|> Symbol
  vxx = view(xx, names .== Symbol(T))
  map(x->x.label, vxx)
end

function lsf(dfg::G, ::Type{T}) where {G <: AbstractDFG, T <: FunctorInferenceType}
  ls(dfg, T)
end

"""
    $(SIGNATURES)
Gives back all factor labels that fit the bill:
	lsWho(dfg, :Pose3)

Dev Notes
- Cloud versions will benefit from less data transfer
 - `ls(dfg::C, ::T) where {C <: CloudDFG, T <: ..}`

Related

ls, lsf, lsfPriors
"""
function lsWho(dfg::AbstractDFG, type::Symbol; solveKey::Symbol=:default)::Vector{Symbol}
    vars = getVariables(dfg)
	labels = Symbol[]
    for v in vars
		varType = typeof(getVariableType(v, solveKey=solveKey)).name |> Symbol
		varType == type && push!(labels, v.label)
	end
	return labels
end

"""
    $(SIGNATURES)
Gives back all factor labels that fit the bill:
	lsfWho(dfg, :Point2Point2)

Dev Notes
- Cloud versions will benefit from less data transfer
 - `ls(dfg::C, ::T) where {C <: CloudDFG, T <: ..}`

Related

ls, lsf, lsfPriors
"""
function lsfWho(dfg::AbstractDFG, type::Symbol)::Vector{Symbol}
	facs = getFactors(dfg)
	labels = Symbol[]
    for f in facs
		facType = typeof(getFactorType(f)).name |> Symbol
		facType == type && push!(labels, f.label)
	end
	return labels
end


text:
"""
    FactorRotationAlgorithm

An abstract type for factor rotation algorithms.
"""
abstract type FactorRotationAlgorithm end

"""
    Orthomax <: FactorRotationAlgorithm

A type representing the orthomax factor rotation algorithm.

The positive parameter `γ::Real` determines which type of rotation is
performed. For a `n x p` matrix, the default `γ = 1.0` leads to varimax
rotation, `γ = 0.0` is quartimax rotation, `γ = p / 2` is equamax
rotation, and `γ = n (p - 1) / (n + p - 2)` is parsimax rotation.

The parameter `maxiter::Integer` controls the maximum number of iterations
to perform (default `1000`), `miniter::Integer` controls the minimum number
of iterations taken before convergence is checked, and `ϵ::Real` is a small
positive constant determining convergence (default `1e-12`).
"""
struct Orthomax <: FactorRotationAlgorithm
    γ::Real
    miniter::Integer
    maxiter::Integer
    ϵ::Real

    Orthomax(;γ::Union{Real, Integer} = 1.0,
              miniter::Integer = 20,
              maxiter::Integer = 1000,
              ϵ::Real = 1e-12) = begin
        γ ≥ zero(eltype(γ)) ||
            throw(DomainError("Orthomax: γ needs to be non-negative"))
        miniter > zero(eltype(miniter)) ||
            throw(DomainError("Orthomax: miniter needs to be positive"))
        maxiter > zero(eltype(maxiter)) ||
            throw(DomainError("Orthomax: maxiter needs to be positive"))
        ϵ > zero(eltype(ϵ)) ||
            throw(DomainError("Orthomax: ϵ needs to be positive"))

        new(float(γ), miniter, maxiter, ϵ)
    end
end

"""
    Varimax() -> Orthomax

Creates an orthomax factor rotation algorithm object for a matrix of
size `n x p` with `γ = 1.0`.
Remaining keyword parameters as for [`Orthomax`](@ref).
"""
function Varimax(;miniter::Integer = 20,
                  maxiter::Integer = 1000,
                  ϵ::Real = 1e-12)
    Orthomax(γ = 1.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ)
end

"""
    Quartimax() -> Orthomax

Creates an orthomax factor rotation algorithm object for a matrix of
size `n x p` with `γ = 0.0`.
Remaining keyword parameters as for [`Orthomax`](@ref).
"""
function Quartimax(;miniter::Integer = 20,
                    maxiter::Integer = 1000,
                    ϵ::Real = 1e-12)
    Orthomax(γ = 0.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ)
end

"""
    Equamax(p) -> Orthomax

Creates an orthomax factor rotation algorithm object for a matrix of
size `n x p` with `γ = p / 2`.
Remaining keyword parameters as for [`Orthomax`](@ref).
"""
function Equamax(p::Integer; miniter::Integer = 20,
                             maxiter::Integer = 1000,
                             ϵ::Real = 1e-12)
    Orthomax(γ = float(p) / 2.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ)
end

"""
    Parsimax(n, p) -> Orthomax

Creates an orthomax factor rotation algorithm object for a matrix of
size `n x p` with `γ = n (p - 1) / (n + p - 2)`.
Remaining keyword parameters as for [`Orthomax`](@ref).
"""
function Parsimax(n::Integer, p::Integer; miniter::Integer = 20,
                                          maxiter::Integer = 1000,
                                          ϵ::Real = 1e-12)
    Orthomax(γ = float(n) * (float(p) - 1) / (float(n) + float(p) - 2),
             miniter = miniter, maxiter = maxiter, ϵ = ϵ)
end

function show(io::IO, mime::MIME{Symbol("text/plain")}, alg::Orthomax)
    summary(io, alg); println(io)
    println(io, "γ = $(alg.γ)")
    println(io, "miniter = $(alg.miniter)")
    println(io, "maxiter = $(alg.maxiter)")
    println(io, "ϵ = $(alg.ϵ)")
end

"""
    FactorRotation{T <: Real}

The return type for factor rotations.

`F` contains the rotated factors and `R` is the rotation matrix that
was applied to the original factors.
"""
struct FactorRotation{T <: Real, Ta <: FactorRotationAlgorithm}
    F::Matrix{T}
    R::Matrix{T}

    alg::Ta

    FactorRotation{T, Ta}(F, R, alg) where {T <: Real, Ta <: FactorRotationAlgorithm} = new{T, Ta}(F, R, alg)
end

eltype(::Type{FactorRotation{T, Ta}}) where {T,Ta} = T

FactorRotation(F::Matrix{T}, R::Matrix{T}, alg::Ta) where {T <: Real, Ta <: FactorRotationAlgorithm} = FactorRotation{T, Ta}(F, R, alg)
FactorRotation(F::Matrix{T}, alg::Ta) where {T <: Real, Ta <: FactorRotationAlgorithm} = FactorRotation(F, Matrix{eltype(F)}(I, size(F, 2), size(F, 2)), alg)

function FactorRotation(F::Matrix, R::Matrix, alg::Ta) where {Ta <: FactorRotationAlgorithm}
    return FactorRotation(promote(F, R)..., alg)
end

function show(io::IO, mime::MIME{Symbol("text/plain")}, FR::FactorRotation{<:Any,<:Any})
    summary(io, FR); println(io)
end

"""
    loadings(FR::FactorRotation)

Returns the loading matrix.
"""
loadings(FR::FactorRotation{T, Ta}) where {T, Ta} = FR.F

"""
    rotation(FR::FactorRotation)

Returns the rotation matrix.
"""
rotation(FR::FactorRotation{T, Ta}) where {T, Ta} = FR.R

## Comparison to other implementations
# The implementation of varimax in R row-normlises the input matrix before
# application of the algorithm and rescales the rows afterwards.
## Reference
# Mikkel B. Stegmann, Karl Sjöstrand, Rasmus Larsen, "Sparse modeling of
# landmark and texture variability using the orthomax criterion,"
# Proc. SPIE 6144, Medical Imaging 2006: Image Processing, 61441G
# (10 March 2006); doi: 10.1117/12.651293
function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ)
    n, p = size(F)
    if n < 2
        return (F, Matrix{eltype(F)}(I, p, p))
    end

    # Warm up step
    # Do one step. If that first step did not lead away from the identity
    # matrix enough use a random orthogonal matrix as a starting point.
    M = svd(F' * (F .^ 3 - γ / n * F * Diagonal(vec(sum(F .^ 2, dims=1)))))
    R = M.U * M.Vt
    if norm(R - Matrix{eltype(R)}(I, p, p)) < ϵ
        R = qr(randn(p, p)).Q
    end

    # Main iterations
    d = 0
    lastchange = NaN
    converged = false
    for i in 1:maxiter
        dold  = d
        B = F * R
        M = svd(F' * (B .^ 3 - γ / n * B * Diagonal(vec(sum(B .^ 2, dims=1)))))
        R = M.U * M.Vt
        d = sum(M.S)
        lastchange = abs(d - dold) / d
        if lastchange < ϵ && i >= miniter
            converged = true
            break
        end
    end

    converged || throw(ConvergenceException(maxiter, lastchange, ϵ))

    (F * R, R)
end

"""
    fit(::Type{FactorRotation}, F::AbstractMatrix; ...) -> FactorRotation

Fit a factor rotation to the matrix `F` and apply it. The algorithm used to
perform the factor rotation is by default [`Orthomax`](@ref) and can be changed
with the keyword argument `alg` which is of type `<:FactorRotationAlgorithm`.
"""
function fit(::Type{FactorRotation}, F::AbstractMatrix;
             alg::T = Orthomax()) where {T <: FactorRotationAlgorithm}
    if isa(alg, Orthomax)
        F, R = orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ)
        return FactorRotation(F, R, alg)
    end
end

"""
    fit(::Type{FactorRotation}, F::FactorAnalysis; ...) -> FactorRotation

Fit a factor rotation to the loading matrix of the  [`FactorAnalysis`](@ref)
object and apply it.

See [`fit(::Type{FactorRotation}, F::AbstractMatrix)`](@ref) for keyword arguments.
"""
function fit(::Type{FactorRotation}, F::FactorAnalysis;
             alg::T = Orthomax()) where {T <: FactorRotationAlgorithm}
    return fit(FactorRotation, F.W; alg = alg)
end

## Alternative interface

"""
    rotatefactors(F::AbstractMatrix, [alg::FactorRotationAlgorithm]) -> FactorRotation

Rotate the factors in the matrix `F`. The algorithm to be used can be passed in
via the second argument `alg`. By default [`Orthomax`](@ref) is used.
"""
function rotatefactors(F::AbstractMatrix, alg::FactorRotationAlgorithm = Orthomax())
    F, R = _rotatefactors(F, alg)
    return FactorRotation(F, R, alg)
end

# Use multiple dispatch to decide on the algorithm implementation

function _rotatefactors(F::AbstractMatrix, alg::Orthomax)
    return orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ)
end

"""
    rotatefactors(F::FactorAnalysis, [alg::FactorRotationAlgorithm]) -> FactorAnalysis

Rotate the factors in the loading matrix of `F` which is of type [`FactorAnalysis`](@ref).
The algorithm to be used can be passed in via the second argument `alg`.
By default [`Orthomax`](@ref) is used.
"""
function rotatefactors(F::FactorAnalysis, alg::FactorRotationAlgorithm = Orthomax())
    FR = rotatefactors(F.W, alg)
    return FactorAnalysis(F.mean, FR.F, F.Ψ)
end


text:
include("NfRel/NfRel.jl")
include("NfRel/NfRelNS.jl")
include("NfRel/NfRelNSMor.jl")
include("NfRel/NfRelOrd.jl")
include("NfRel/OrdElem.jl")
include("NfRel/Ideal.jl")
include("NfRel/FracIdeal.jl")
include("../Map/NfRelOrd.jl")
include("NfRel/ResidueRing.jl")
include("NfRel/NEQ.jl")
include("NfRel/Conjugates.jl")
include("NfRel/Simplify.jl")
include("NfRel/absolute_field.jl")
include("NfRel/relative_extension.jl")


text:
############################
# Lorenz System  (MIMO)    #
# by VSSM (all directions) #
############################

###########################################
# Controll parameters for VSSM controller #
###########################################

K=500
w=1

#############################################
# Parameters to design a Nominal Trajectory #
#############################################

A₁=2
ω₁=0.5
A₂=3.8
ω₂=0.7
A₃=1.5
ω₃=0.9

#################
# Time variable #
#################

# cycle time in seconds
δt=1e-3
# The lenght of the simulation in seconds
LONG=Int(2e4)

Λ=12

##########################
# Exact model parameters #
##########################
σₑ=5
βₑ=8/3
ρₑ=40

############################
# Approx. model parameters #
############################
σₐ=5
βₐ=8/3
ρₐ=40

#Arrays

xN=zeros(LONG)
xN_p=zeros(LONG)
yN=zeros(LONG)
yN_p=zeros(LONG)
zN=zeros(LONG)
zN_p=zeros(LONG)

x_p_des=zeros(LONG)
y_p_des=zeros(LONG)
z_p_des=zeros(LONG)

u_x=zeros(LONG)
u_y=zeros(LONG)
u_z=zeros(LONG)

t=zeros(LONG)

S_x=zeros(LONG)
S_p_x=zeros(LONG)
S_y=zeros(LONG)
S_p_y=zeros(LONG)
S_z=zeros(LONG)
S_p_z=zeros(LONG)

x=zeros(LONG)
y=zeros(LONG)
z=zeros(LONG)
x_p=zeros(LONG)
y_p=zeros(LONG)
z_p=zeros(LONG)

#initial conditions
hint_x=0
hint_y=0
hint_z=0

x[1]=A₁*sin(ω₁*δt)
y[1]=A₂*sin(ω₂*δt)
z[1]=A₃*sin(ω₃*δt)

##############
# SIMULATION #
##############

for i=1:LONG-1 #i counts the cycles of the simulation
  t[i]=δt*i
  #create Nominal Trajectory
  xN[i]=A₁*sin(ω₁*t[i])
  xN_p[i]=A₁*ω₁*cos(ω₁*t[i])
  #xN_p[i]=-A₁*ω₁^2*sin(ω₁*t[i])

  yN[i]=A₂*sin(ω₂*t[i])
  yN_p[i]=A₂*ω₂*cos(ω₂*t[i])
  #yN_p[i]=-A₂*ω₂^2*sin(ω₂*t[i])

  zN[i]=A₃*sin(ω₃*t[i])
  zN_p[i]=A₃*ω₃*cos(ω₃*t[i])
  #zN_p[i]=-A₃*ω₃^2*sin(ω₃*t[i])
  # Compute  error (h) and error integral (hint) and the first derivative of the error (h_p)
  h_x=xN[i]-x[i]
  h_y=yN[i]-y[i]
  h_z=zN[i]-z[i]

  hint_x=hint_x+δt*h_x
  hint_y=hint_y+δt*h_y
  hint_z=hint_z+δt*h_z

  h_p_x=xN_p[i]-x_p[i]
  h_p_y=yN_p[i]-y_p[i]
  h_p_z=zN_p[i]-z_p[i]

  #Error metric (S), and the first derivative of the Error metric (S_p)
  S_x[i]=Λ*hint_x+h_x
  S_p_x[i]=Λ*h_x+h_p_x
  S_y[i]=Λ*hint_y+h_y
  S_p_y[i]=Λ*h_y+h_p_y
  S_z[i]=Λ*hint_z+h_z
  S_p_z[i]=Λ*h_z+h_p_z
  #Kinematic Block (PID type feedback)
  x_p_des[i]=xN_p[i]+Λ*h_x+K*tanh(S_x[i]/w)
  y_p_des[i]=yN_p[i]+Λ*h_y+K*tanh(S_y[i]/w)
  z_p_des[i]=zN_p[i]+Λ*h_z+K*tanh(S_z[i]/w)
  # Create the control signal
  u_x[i]=x_p_des[i]-σₐ*(y[i]-x[i])
  u_y[i]=y_p_des[i]-x[i]*(ρₐ-z[i])
  u_z[i]=z_p_des[i]-y[i]+βₐ*z[i]
  # Get the system "reaction" (the exact model of the system)
  x_p[i]=σₑ*(y[i]-x[i])+u_x[i]
  y_p[i]=x[i]*(ρₑ-z[i])-y[i]+u_y[i]
  z_p[i]=x[i]*y[i]-βₑ*z[i]+u_z[i]
  # Integrate back to x, y, z
  x[i+1]=x[i]+δt*x_p[i]
  y[i+1]=y[i]+δt*y_p[i]
  z[i+1]=z[i]+δt*z_p[i]
end #for

using PyPlot

figure(1)
grid("on")
title("Nominal and Realized trajectorys")
xlabel(L"Time $[s]$")
ylabel(L"q $[m]$")
plot(t[3:LONG-1],xN[3:LONG-1],color="#7684FF")
plot(t[3:LONG-1],yN[3:LONG-1],color="#2A40FF")
plot(t[3:LONG-1],zN[3:LONG-1],color="#3B427F")

plot(t[3:LONG-1],x[3:LONG-1],color="#FFAA41","r--")
plot(t[3:LONG-1],y[3:LONG-1],color="#FF9F28","r--")
plot(t[3:LONG-1],z[3:LONG-1],color="#B2690E","r--")

figure(2)
grid("on")
title("Tracking Error")
xlabel(L"Time $[s]$")
ylabel(L"error $[m]$")
plot(t[3:LONG-1],xN[3:LONG-1]-x[3:LONG-1],color="red")
plot(t[3:LONG-1],yN[3:LONG-1]-y[3:LONG-1],color="green")
plot(t[3:LONG-1],zN[3:LONG-1]-z[3:LONG-1],color="blue")

figure(3)
grid("on")
title("Phase Space of Error Metric (S_p vs S)")
plot(S_x[3:LONG-1],S_p_x[3:LONG-1],color="red")
plot(S_y[3:LONG-1],S_p_y[3:LONG-1],color="green","r--")
plot(S_z[3:LONG-1],S_p_z[3:LONG-1],color="blue","r-.")

figure(4)
grid("on")
title("Phase Space of motion (x_p vs x)")
plot(xN[1:LONG-1],xN_p[1:LONG-1],color="red")
plot(x[1:LONG-1],x_p[1:LONG-1],color="blue","r--")
plot(yN[1:LONG-1],yN_p[1:LONG-1],color="#7684FF")
plot(y[1:LONG-1],y_p[1:LONG-1],color="#FFAA41","r--")
plot(zN[1:LONG-1],zN_p[1:LONG-1],color="#2A40FF")
plot(z[1:LONG-1],z_p[1:LONG-1],color="#FF9F28","r-.")

max_stars_repo_name: krisztiankosi/Control

max_issues_repo_name: krisztiankosi/Control

max_forks_repo_name: krisztiankosi/Control


text:
"""
    gams = parsegams(filename)

Do low-level lexing and parsing, but do not try to compile into a module.
"""
function parsegams(lexed::Vector{AbstractLex})
    gams = Dict{String,Any}()
    i = 1
    while i <= length(lexed)
        stmt = lexed[i]
        if stmt isa StatementEnd
            # do nothing
        elseif stmt isa Keyword
            kwname = stmt.text
            if kwname == "sets"
                sets = getdefault!(gams, kwname, Dict{String,Any})
                newsets = Dict{String,Vector{String}}()
                i = parse_slashed!(newsets, lexed, i)
                # Scan for i*j range declarations
                for (k, s) in newsets
                    if length(s) == 1
                        m = match(r"(\d+)\s*\*\s*(\d+)", s[1])
                        if m != nothing
                            lo, hi = parse(Int, m.captures[1]), parse(Int, m.captures[2])
                            sets[k] = lo == 1 ? Base.OneTo(hi) : (lo:hi)
                        else
                            sets[k] = s
                        end
                    else
                        sets[k] = s
                    end
                end
            elseif kwname == "alias"
                alias = lexed[i+=1]
                @assert(alias isa Parens && length(alias.args) == 3 && alias.args[2] == GText(","))
                sets = gams["sets"]
                symold, symnew = string(alias.args[1]), string(alias.args[3])
                if !haskey(sets, symold)  # GAMS allows arbitrary ordering
                    symnew, symold = symold, symnew
                end
                sets[symnew] = sets[symold]
            elseif kwname == "scalars" || kwname == "parameters"
                params = getdefault!(gams, "parameters", Dict{Any,Any})
                i = parse_slashed!(params, lexed, i; parseitems=false)
            elseif kwname == "table"
                table = getdefault!(gams, "tables", Dict{Any,Any})
                item = lexed[i+=1]
                table[item.name] = item.body
                # Some GAMS files declare it a parameter and then a table; remove it
                # from parameters
                if haskey(gams, "parameters")
                    params = gams["parameters"]
                    key = GText(getname(item.name))
                    if haskey(params, key)
                        delete!(params, key)
                    end
                end
            elseif kwname == "variables" || kwname ∈ vartypes
                attr = "free"
                if kwname ∈ vartypes
                    attr = kwname
                    kw = lexed[i+=1]
                    @assert(kw == Keyword("variables"))
                end
                sets = get(gams, "sets", nothing)
                variables = getdefault!(gams, "variables", OrderedDict{Any,VarInfo})
                while i < length(lexed) && !isa(lexed[i+1], StatementEnd) && !isa(lexed[i+1], Keyword)
                    var = lexed[i+=1]
                    var == GText(",") && continue
                    @assert(var isa GText || var isa GArray)
                    push!(variables, varinfo(var, attr, sets))
                end
            elseif kwname == "models"
                models = getdefault!(gams, kwname, Dict{String,ModelInfo})
                newmodels = Dict{String,Any}()
                i = parse_slashed!(newmodels, lexed, i)
                for (k,v) in newmodels
                    models[k] = ModelInfo(v)
                end
            elseif kwname == "options"
                options = getdefault!(gams, kwname)
                str = lexed[i+=1]
                while !isa(lexed[i+1], StatementEnd)
                    if lexed[i+1] == GText("=")
                        options[str.text] = string(lexed[i+2])
                        i+=1
                    else
                        options[str.text] = ""
                    end
                    i += 1
                end
            elseif kwname == "solve"
                newkw = lexed[i+=1]
                while !isa(newkw, Keyword) || newkw.text == "using"
                    newkw = lexed[i+=1]
                end
                @assert newkw.text == "minimizing" || newkw.text == "maximizing"
                modelname = lexed[i+=1]
                gams[newkw.text] = modelname.text
            elseif kwname == "equations"
                eqs = getdefault!(gams, kwname, OrderedDict{Any,Any})
                while i < length(lexed) && !isa(lexed[i+1], StatementEnd) && !isa(lexed[i+1], Keyword)
                    eq = lexed[i+=1]
                    @assert(eq isa GText || eq isa GArray)
                    # push!(eqs, eq=>"")  # add later so the ordering respects the definition order, not decl order
                end
            elseif kwname == "files"
                filename = lexed[i+=1]
                newkw = lexed[i+=1]
                while newkw != GText("putclose") || lexed[i+1] != filename
                    newkw = lexed[i+=1]
                end
                newkw = lexed[i+=1]
                while !isa(newkw, StatementEnd)
                    newkw = lexed[i+=1]
                end
            elseif kwname ∈ Set(["display"])
                newkw = lexed[i+=1]
                while !isa(newkw, StatementEnd)
                    newkw = lexed[i+=1]
                end
            else
                error("unhandled keyword ", kwname)
            end
        else
            @assert(stmt isa GText || stmt isa GArray)
            if length(lexed) < i+1
                # some files seem to have junk after the last semicolon
                println("skipping trailing ", stmt)
                break
            end
            if lexed[i+1] == Dots("..")
                # an equation definition
                eqs = gams["equations"]
                iend = seek_to_end(lexed, i+=2)
                eqs[stmt] = parseexprs1(lexed[i:iend])
                i = iend
            elseif lexed[i+1] == Dots(".")
                # a variable or model attribute
                @assert(lexed[i+3] == GText("="))
                attr = lexed[i+=2]
                iend = seek_to_end(lexed, i+=2)
                rhs = parseexprs1(lexed[i:iend])
                if attr isa GText && string(attr) ∈ modelattributes
                    m = gams["models"][string(stmt)]
                    m.assignments[string(attr)] = rhs
                else
                    attrstr = isa(attr, GText) ? string(attr) : attr.name
                    varkey = isa(attr, GText) ? stmt : GArray(string(stmt), attr.indices)
                    if attrstr ∈ varattributes
                        v = gams["variables"][varkey]  # gams["variables"][stmt] ?
                        push!(v.assignments, attr=>rhs)
                    else
                        error(attr, " not a recognized attribute")
                    end
                end
                i = iend
            else
                # general assignment (to, e.g., a parameter)
                @assert lexed[i+1] == GText("=")
                iend = seek_to_end(lexed, i+=2)
                assign = getdefault!(gams, "assignments", Vector{Pair{Any,Any}})
                push!(assign, stmt=>parseexprs1(lexed[i:iend]))
                i = iend
            end
        end
        i += 1
    end
    return gams
end
parsegams(arg) = parsegams(lex(arg))

function parse_slashed!(dest, lexed, i; parseitems::Bool=true)
    function parse_slashed1!(dest, lexed, i)
        name = lexed[i+=1]
        @assert(name isa GText || name isa GArray)
        if lexed[i+1] isa StatementEnd
            dest[name] = ""
            return i
        elseif lexed[i+1] isa Slashed
            val = lexed[i+=1]
        elseif lexed[i+1] isa GText && lexed[i+2] isa Slashed
            val = lexed[i+=2]
        else
            dest[name] = ""
            return i
        end
        @assert(val isa Slashed)
        if parseitems
            items = strip.(split(val.text, (',','\n')))
            # drop text for each item
            for j = 1:length(items)
                itemname, itemtext = splitws(items[j])
                items[j] = itemname
            end
            dest[convert(keytype(dest), name)] = items
        else
            dest[name] = val.text
        end
        return i
    end
    i = parse_slashed1!(dest, lexed, i)
    while i < length(lexed) && !isa(lexed[i+1], StatementEnd) && !isa(lexed[i+1], Keyword)
        i = parse_slashed1!(dest, lexed, i)
    end
    return i
end

function parseexprs(lexed, debug=false)
    # Step 1: split at commas
    args = []
    thisarg = []
    for i = 1:length(lexed)
        tok = lexed[i]
        if tok == GText(",")
            @assert(!isempty(thisarg))
            push!(args, copy(thisarg))
            empty!(thisarg)
        else
            push!(thisarg, tok)
        end
    end
    if !isempty(thisarg)
        push!(args, copy(thisarg))
        empty!(thisarg)
    end
    # Step 2: parse each arg
    n = length(args)
    parsed = Vector{AbstractLex}(undef, n)
    todelete = Int[]
    for j = 1:n
        thisarg = args[j]
        # Step 2a: recurse into calls and parse their args
        debug && sexprlf.(thisarg)
        for i = 1:length(thisarg)
            tok = thisarg[i]
            if tok isa GCall
                callargs = parseexprs(tok.args, false)
                thisarg[i] = GCall(tok.name, callargs)
            elseif tok isa Parens
                pargs = parseexprs(tok.args, false)
                thisarg[i] = Parens(pargs)
            end
        end
        # Step 2b: process ^ (highest-precedence operator)
        debug && sexprlf.(thisarg)
        process_ops!(thisarg, ("^",), debug)
        # Step 2c: process * and / (next-highest)
        debug && sexprlf.(thisarg)
        process_ops!(thisarg, ("*", "/"), debug)
        # Step 2d: process + and -
        debug && println("prior to +/-")
        debug && sexprlf.(thisarg)
        process_ops!(thisarg, ("+", "-"), debug)
        # Step 2e: process "and" and "or"
        debug && println("post to +/-")
        debug && sexprlf.(thisarg)
        process_ops!(thisarg, ("and", "or"), debug)
        # Step 2f: process other boolean-returning operations
        process_ops!(thisarg, ("<", "<=", ">", ">="), debug)
        # Step 2g: process relations
        debug && sexprlf.(thisarg)
        process_ops!(thisarg, ("=e=", "=g=", "=l="), debug)
        if length(thisarg) != 1
            println(stderr, "thisarg:")
            foreach(x->(print(stderr, "  "); sexpr(stderr, x); println(stderr)), thisarg)
        end
        @assert(length(thisarg) == 1)
        parsed[j] = stripnewcall!(thisarg[1])
    end
    return parsed
end

function parseexprs1(lexed)
    parsed = parseexprs(lexed)
    @assert(length(parsed) == 1)
    return parsed[1]
end

# A "temporary type" for marking calls that arose from operator-precedence resolution
# Future lower-priority operators will need to recurse into any NewCalls
struct NewCall <: AbstractLex
    call
end

function process_ops!(tokens::Vector, ops, debug)
    debug = false
    # i = length(tokens)
    i = 1
    # while i > 0
    while i < length(tokens)
        debug && @show i
        tok = tokens[i]
        debug && @show tok
        if tok isa NewCall
            newtok = tok.call
            for arg in newtok.args
                process_ops!(arg, ops)
            end
        end
        if !isa(tok, GText) || tok.text ∉ ops
            i += 1
            continue
        end
        if i > 1 && !isoperator(tokens[i-1])
            debug && @show tokens[i-1] isoperator(tokens[i-1])
            tokens[i-1] = NewCall(GCall(tok.text, [tokens[i-1], tokens[i+1]]))
            deleteat!(tokens, (i,i+1))
            i -= 1
        else
            # unary + or -
            @assert(tok.text ∈ ("+", "-"))
            if tok.text == "+"
                tokens[i] = tokens[i+1]
            else
                nexttok = tokens[i+1]
                if nexttok isa GNumber
                    tokens[i] = GNumber(-nexttok.val)
                else
                    tokens[i] = NewCall(GCall("-", [nexttok]))
                end
            end
            deleteat!(tokens, (i+1,))
        end
        i += 1
    end
    return tokens
end
process_ops!(tokens, ops) = nothing

isoperator(op) = op isa GText && op.text ∈ operatortokens

function sexpr(io::IO, l::NewCall)
    print(io, "NewCall(")
    sexpr(io, l.call)
    print(io, ")")
end

function stripnewcall!(expr::NewCall)
    call = expr.call
    for i = 1:length(call.args)
        call.args[i] = stripnewcall!(call.args[i])
    end
    return call
end
stripnewcall!(expr) = expr

function deparen!(expr)
    if expr isa GCall
        for i = 1:length(expr.args)
            expr.args[i] = deparen!(expr.args[i])
        end
    elseif expr isa Parens && length(expr.args) == 1
        expr = expr.args[1]
        expr = deparen!(expr)
    end
    return expr
end

function varinfo(name::GText, attr, sets)
    return name=>VarInfo(attr, ())
end
function varinfo(name::GArray, attr, sets)
    return name=>VarInfo(attr, getaxes(name.indices, sets))
end

function seek_to_end(lexed, i)
    iend = i
    while iend < length(lexed) && !isa(lexed[iend+1], Union{StatementEnd, Keyword})
        iend += 1
    end
    return iend
end

function getwithkey(dict, key::AbstractString)
    for (k, v) in dict
        if getname(k) == key
            return v
        end
    end
    error(key, " not found")
end

function getdefault!(gams, tag, ::Type{T}) where T
    if !haskey(gams, tag)
        gams[tag] = T<:AbstractDict ? T() : (T<:Vector ? T(undef, 0) : error("type $T not handled"))
    end
    return gams[tag]
end
getdefault!(gams, tag) = getdefault!(gams, tag, Dict{String,String})