text:
abstract type StochasticDiffEqAlgorithm <: AbstractSDEAlgorithm end
abstract type StochasticDiffEqAdaptiveAlgorithm <: StochasticDiffEqAlgorithm end
abstract type StochasticDiffEqCompositeAlgorithm <: StochasticDiffEqAlgorithm end

abstract type StochasticDiffEqRODEAlgorithm <: AbstractRODEAlgorithm end
abstract type StochasticDiffEqRODEAdaptiveAlgorithm <: StochasticDiffEqRODEAlgorithm end
abstract type StochasticDiffEqRODECompositeAlgorithm <: StochasticDiffEqRODEAlgorithm end

abstract type StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,FDT,ST,Controller} <: StochasticDiffEqAdaptiveAlgorithm end
abstract type StochasticDiffEqNewtonAlgorithm{CS,AD,FDT,ST,Controller} <: StochasticDiffEqAlgorithm end

abstract type StochasticDiffEqJumpAlgorithm <: StochasticDiffEqAlgorithm end
abstract type StochasticDiffEqJumpAdaptiveAlgorithm <: StochasticDiffEqAlgorithm end
abstract type StochasticDiffEqJumpNewtonAdaptiveAlgorithm{CS,AD,FDT,ST,Controller} <: StochasticDiffEqJumpAdaptiveAlgorithm end

abstract type StochasticDiffEqJumpDiffusionAlgorithm <: StochasticDiffEqAlgorithm end
abstract type StochasticDiffEqJumpDiffusionAdaptiveAlgorithm <: StochasticDiffEqAlgorithm end
abstract type StochasticDiffEqJumpNewtonDiffusionAdaptiveAlgorithm{CS,AD,FDT,ST,Controller} <: StochasticDiffEqJumpDiffusionAdaptiveAlgorithm end

abstract type IteratedIntegralApprox end
struct IICommutative <:  IteratedIntegralApprox end
struct IIWiktorsson <:  IteratedIntegralApprox end

################################################################################

# Basics
"""
EM: Nonstiff Method
The Euler-Maruyama method. Strong Order 0.5 in the Ito sense.
Has an optional argument split=true for controlling step splitting.
When splitting is enabled, the stability with large diffusion eigenvalues is improved.
Can handle all forms of noise, including non-diagonal, scalar, and colored noise.
Fixed time step only.
"""
struct EM{split} <: StochasticDiffEqAlgorithm end
EM(split=true) = EM{split}()

struct SplitEM <: StochasticDiffEqAlgorithm end
"""
EulerHeun: Nonstiff Method
The Euler-Heun method.
Strong Order 0.5 in the Stratonovich sense.
Can handle all forms of noise, including non-diagonal, scalar, and colored noise.
Fixed time step only.
"""
struct EulerHeun <: StochasticDiffEqAlgorithm end
"""
LambaEM: Nonstiff Method
A modified Euler-Maruyama method with adaptive time stepping with an error estimator based on Lamba and Rackauckas.
Has an optional argument split=true for controlling step splitting.
When splitting is enabled, the stability with large diffusion eigenvalues is improved.
Strong Order 0.5 in the Ito sense. Can handle all forms of noise, including non-diagonal, scalar, and colored noise
"""
struct LambaEM{split} <: StochasticDiffEqAdaptiveAlgorithm end
LambaEM(split=true) = LambaEM{split}()
"""
LambaEulerHeun: Nonstiff Method
A modified Euler-Heun method with adaptive time stepping with an error estimator based on Lamba due to Rackauckas.
Strong order 0.5 in the Stratonovich sense.
Can handle all forms of noise, including non-diagonal, scalar, and colored noise.
"""
struct LambaEulerHeun <: StochasticDiffEqAdaptiveAlgorithm end

"""
Kloeden, P.E., Platen, E., Numerical Solution of Stochastic Differential Equations.
Springer. Berlin Heidelberg (2011)

SimplifiedEM: High Weak Order Method
A simplified Euler-Maruyama method with weak order 1.0 and fixed step size.
Can handle all forms of noise, including non-diagonal, scalar, and colored noise.
"""
struct SimplifiedEM <: StochasticDiffEqAlgorithm end

"""
Kloeden, P.E., Platen, E., Numerical Solution of Stochastic Differential Equations.
Springer. Berlin Heidelberg (2011)

RKMil: Nonstiff Method
An explicit Runge-Kutta discretization of the strong order 1.0 Milstein method.
Defaults to solving the Ito problem, but RKMil(interpretation=:Stratonovich) makes it solve the Stratonovich problem.
Only handles scalar and diagonal noise.
"""
struct RKMil{interpretation} <: StochasticDiffEqAdaptiveAlgorithm end
RKMil(;interpretation=:Ito) = RKMil{interpretation}()

"""
Kloeden, P.E., Platen, E., Numerical Solution of Stochastic Differential Equations.
Springer. Berlin Heidelberg (2011)

RKMilCommute: Nonstiff Method
An explicit Runge-Kutta discretization of the strong order 1.0 Milstein method for commutative noise problems.
Defaults to solving the Ito problem, but RKMilCommute(interpretation=:Stratonovich) makes it solve the Stratonovich problem.
Uses a 1.5/2.0 error estimate for adaptive time stepping.
"""
struct RKMilCommute{interpretation} <: StochasticDiffEqAdaptiveAlgorithm end
RKMilCommute(;interpretation=:Ito) = RKMilCommute{interpretation}()

"""
Kloeden, P.E., Platen, E., Numerical Solution of Stochastic Differential Equations.
Springer. Berlin Heidelberg (2011)

RKMilGeneral: Nonstiff Method
RKMilGeneral(;interpretation=:Ito, ii_approx=IIWiktorsson()
An explicit Runge-Kutta discretization of the strong order 1.0 Milstein method for general non-commutative noise problems.
Allows for a choice of interpretation between :Ito and :Stratonovich.
Allows for a choice of iterated integral approximation.
"""
struct RKMilGeneral{T<:IteratedIntegralApprox, TruncationType} <: StochasticDiffEqAdaptiveAlgorithm
  interpretation::Symbol
  ii_approx::T
  c::Int
  p::TruncationType
end
function RKMilGeneral(;interpretation=:Ito, ii_approx=IIWiktorsson(), c=1, p=nothing, dt=nothing)
  γ = 1//1
  p==true && (p = Int(floor(c*dt^(1//1-2//1*γ)) + 1))
  RKMilGeneral(interpretation, ii_approx, c, p)
end
"""
WangLi3SMil_A: Nonstiff Method
Fixed step-size explicit 3-stage Milstein methods for Ito problem with strong and weak order 1.0
"""
struct WangLi3SMil_A <: StochasticDiffEqAlgorithm end
"""
WangLi3SMil_B: Nonstiff Method
Fixed step-size explicit 3-stage Milstein methods for Ito problem with strong and weak order 1.0
"""
struct WangLi3SMil_B <: StochasticDiffEqAlgorithm end
"""
WangLi3SMil_C: Nonstiff Method
Fixed step-size explicit 3-stage Milstein methods for Ito problem with strong and weak order 1.0
"""
struct WangLi3SMil_C <: StochasticDiffEqAlgorithm end
"""
WangLi3SMil_D: Nonstiff Method
Fixed step-size explicit 3-stage Milstein methods for Ito problem with strong and weak order 1.0
"""
struct WangLi3SMil_D <: StochasticDiffEqAlgorithm end
"""
WangLi3SMil_E: Nonstiff Method
Fixed step-size explicit 3-stage Milstein methods for Ito problem with strong and weak order 1.0
"""
struct WangLi3SMil_E <: StochasticDiffEqAlgorithm end
"""
WangLi3SMil_F: Nonstiff Method
Fixed step-size explicit 3-stage Milstein methods for Ito problem with strong and weak order 1.0
"""
struct WangLi3SMil_F <: StochasticDiffEqAlgorithm end

#SROCK methods
"""
SROCK1: S-ROCK Method
Is a fixed step size stabilized explicit method for stiff problems.
Defaults to solving th Ito problem but SROCK1(interpretation=:Stratonovich) can make it solve the Stratonovich problem.
Strong order of convergence is 0.5 and weak order 1, but is optimised to get order 1 in case os scalar/diagonal noise.
"""
struct SROCK1{interpretation,E} <: StochasticDiffEqAlgorithm
  eigen_est::E
end
SROCK1(;interpretation=:Ito,eigen_est=nothing) = SROCK1{interpretation,typeof(eigen_est)}(eigen_est)

# Weak Order 2
for Alg in [:SROCK2, :KomBurSROCK2, :SROCKC2]
  @eval begin
    struct $Alg{E} <: StochasticDiffEqAlgorithm
      eigen_est::E
    end
    $Alg(;eigen_est=nothing) = $Alg(eigen_est)
  end
end

# ROCK stabilization for EM
"""
SROCKEM: S-ROCK Method
Is fixed step Euler-Mayurama with first order ROCK stabilization thus can handle stiff problems.
Only for Ito problems. Defaults to strong and weak order 1.0, but can solve with weak order 0.5 as SROCKEM(strong_order_1=false).
This method can handle 1-dimensional, diagonal and multi-dimensional noise.
"""
struct SROCKEM{E} <: StochasticDiffEqAlgorithm
  strong_order_1::Bool
  eigen_est::E
end
SROCKEM(;strong_order_1=true,eigen_est=nothing) = SROCKEM(strong_order_1,eigen_est)
"""
SKSROCK: S-ROCK Method
Is fixed step stabilized explicit method for stiff Ito problems.
Strong order 0.5 and weak order 1.
This method has a better stability domain then SROCK1.
Also it allows special post-processing techniques in case of ergodic dynamical systems, in the context of ergodic Brownian dynamics, to achieve order 2 accuracy.
SKSROCK(;post_processing=true) will make use of post processing.
By default it doesn't use post processing.
Post processing is optional and under development.
The rest of the method is completely functional and can handle 1-dimensional, diagonal and multi-dimensional noise.
"""
struct SKSROCK{E} <: StochasticDiffEqAlgorithm
  post_processing::Bool
  eigen_est::E
end
SKSROCK(;post_processing=false,eigen_est=nothing) = SKSROCK(post_processing,eigen_est)
"""
TangXiaoSROCK2: S-ROCK Method
Is a fixed step size stabilized expicit method for stiff problems.
Only for Ito problems. Weak order of 2 and strog order of 1.
Has 5 versions with different stability domains which can be used as TangXiaoSROCK2(version_num=i) where i is 1-5. Under Development.
"""
struct TangXiaoSROCK2{E} <: StochasticDiffEqAlgorithm
  version_num::Int
  eigen_est::E
end
TangXiaoSROCK2(;version_num=5,eigen_est=nothing) = TangXiaoSROCK2(version_num,eigen_est)
###############################################################################

# Predictor Corrector
struct PCEuler{T<:Real, F} <: StochasticDiffEqAlgorithm
  theta::T
  eta::T
  ggprime::F
end

"""
    PCEuler(ggprime; theta=1/2, eta=1/2)

Predictor Corrector Euler

# Arguments
- `ggprime::Function`:
  For scalar problems, `ggprime` ``= b\\partial_x(b)``
  For multi-dimensional problems
  `bbprime_k` ``= \\sum_{j=1...M, i=1...D} b^(j)_i \\partial_i b^(j)_k``
  where ``b^(j)`` correspond to the noise vector due to the j'th noise channel.
  If problem is in place - a in place ggprime should be supplied - and
  vice versa for not in place speicification of problem.
- `theta::Real`:
  Degree of implicitness in the drift term. Set to 0.5 by default.
- `eta::Real`:
  Degree of implicitness in the diffusion term. Set to 0.5 by default.

Reference: Stochastics and Dynamics, Vol. 8, No. 3 (2008) 561–581
Note that the original paper has a typo in the definition of ggprime...
"""
PCEuler(ggprime; theta=1/2, eta=1/2) = PCEuler(theta,eta,ggprime)

################################################################################

# Rossler

"""
Rößler A., Runge–Kutta Methods for the Strong Approximation of Solutions of
Stochastic Differential Equations, SIAM J. Numer. Anal., 48 (3), pp. 922–952.
DOI:10.1137/09076636X

SRA: Nonstiff Method
Adaptive strong order 1.5 methods for additive Ito and Stratonovich SDEs.
Default tableau is for SRA1. Can handle diagonal, non-diagonal and scalar additive noise.
"""
struct SRA{TabType} <: StochasticDiffEqAdaptiveAlgorithm
  tableau::TabType
end
SRA(;tableau=constructSRA1()) = SRA(tableau)

"""
Rößler A., Runge–Kutta Methods for the Strong Approximation of Solutions of
Stochastic Differential Equations, SIAM J. Numer. Anal., 48 (3), pp. 922–952.
DOI:10.1137/09076636X

SRI: Nonstiff Method
Adaptive strong order 1.5 methods for diagonal/scalar Ito SDEs.
Default tableau is for SRIW1.
"""
struct SRI{TabType} <: StochasticDiffEqAdaptiveAlgorithm
  tableau::TabType
  error_terms::Int
end
SRI(;tableau=constructSRIW1(),error_terms=4) = SRI(tableau,error_terms)

"""
Rößler A., Runge–Kutta Methods for the Strong Approximation of Solutions of
Stochastic Differential Equations, SIAM J. Numer. Anal., 48 (3), pp. 922–952.
DOI:10.1137/09076636X

SRIW1: Nonstiff Method
Adaptive strong order 1.5 and weak order 2.0 for diagonal/scalar Ito SDEs.
"""
struct SRIW1 <: StochasticDiffEqAdaptiveAlgorithm end

"""
Rößler A., Runge–Kutta Methods for the Strong Approximation of Solutions of
Stochastic Differential Equations, SIAM J. Numer. Anal., 48 (3), pp. 922–952.
DOI:10.1137/09076636X

SRIW2: Nonstiff Method
Adaptive strong order 1.5 and weak order 3.0 for diagonal/scalar Ito SDEs.
"""
struct SRIW2 <: StochasticDiffEqAdaptiveAlgorithm end
"""
SOSRI: Nonstiff Method
Stability-optimized adaptive strong order 1.5 and weak order 2.0 for diagonal/scalar Ito SDEs.
Stable at high tolerances and robust to stiffness.
"""
struct SOSRI <: StochasticDiffEqAdaptiveAlgorithm end
"""
SOSRI2: Nonstiff Method
Stability-optimized adaptive strong order 1.5 and weak order 2.0 for diagonal/scalar Ito SDEs.
Stable at high tolerances and robust to stiffness.
"""
struct SOSRI2 <: StochasticDiffEqAdaptiveAlgorithm end

"""
Rößler A., Runge–Kutta Methods for the Strong Approximation of Solutions of
Stochastic Differential Equations, SIAM J. Numer. Anal., 48 (3), pp. 922–952.
DOI:10.1137/09076636X

SRA1: Nonstiff Method
Adaptive strong order 1.5 for additive Ito and Stratonovich SDEs with weak order 2.
Can handle diagonal, non-diagonal, and scalar additive noise.
"""
struct SRA1 <: StochasticDiffEqAdaptiveAlgorithm end

"""
Rößler A., Runge–Kutta Methods for the Strong Approximation of Solutions of
Stochastic Differential Equations, SIAM J. Numer. Anal., 48 (3), pp. 922–952.
DOI:10.1137/09076636X

SRA2: Nonstiff Method
Adaptive strong order 1.5 for additive Ito and Stratonovich SDEs with weak order 2.
Can handle diagonal, non-diagonal, and scalar additive noise.
"""
struct SRA2 <: StochasticDiffEqAdaptiveAlgorithm end

"""
Rößler A., Runge–Kutta Methods for the Strong Approximation of Solutions of
Stochastic Differential Equations, SIAM J. Numer. Anal., 48 (3), pp. 922–952.
DOI:10.1137/09076636X

SRA3: Nonstiff Method
Adaptive strong order 1.5 for additive Ito and Stratonovich SDEs with weak order 3.
Can handle non-diagonal and scalar additive noise.
"""
struct SRA3 <: StochasticDiffEqAdaptiveAlgorithm end
"""
SOSRA: Nonstiff Method
A stability-optimized adaptive SRA. Strong order 1.5 for additive Ito and Stratonovich SDEs with weak order 2.
Can handle diagonal, non-diagonal, and scalar additive noise.
Stable at high tolerances and robust to stiffness.
"""
struct SOSRA <: StochasticDiffEqAdaptiveAlgorithm end
"""
SOSRA2: Nonstiff Method
A stability-optimized adaptive SRA. Strong order 1.5 for additive Ito and Stratonovich SDEs with weak order 2.
Can handle diagonal, non-diagonal, and scalar additive noise.
Stable at high tolerances and robust to stiffness.
"""
struct SOSRA2 <: StochasticDiffEqAdaptiveAlgorithm end

################################################################################

# Rossler second order for weak approx.

"""
Debrabant, K. and Rößler A., Families of efficient second order Runge–Kutta methods
for the weak approximation of Itô stochastic differential equations,
Applied Numerical Mathematics 59, pp. 582–594 (2009)
DOI:10.1016/j.apnum.2008.03.012

DRI1: High Weak Order Method
Adaptive step weak order 2.0 for Ito SDEs with minimized error constants (deterministic order 3).
Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise.
"""
struct DRI1 <: StochasticDiffEqAdaptiveAlgorithm end

"""
Debrabant, K. and Rößler A., Families of efficient second order Runge–Kutta methods
for the weak approximation of Itô stochastic differential equations,
Applied Numerical Mathematics 59, pp. 582–594 (2009)
DOI:10.1016/j.apnum.2008.03.012

DRI1NM: High Weak Order Method
Adaptive step weak order 2.0 for Ito SDEs with minimized error constants (deterministic order 3).
Can handle non-mixing diagonal (i.e., du[k] = f(u[k])) and scalar additive noise.
"""
struct DRI1NM <: StochasticDiffEqAdaptiveAlgorithm end

"""
Rößler A., Second Order Runge–Kutta Methods for Itô Stochastic Differential Equations,
SIAM J. Numer. Anal., 47, pp. 1713-1738 (2009)
DOI:10.1137/060673308

RI1: High Weak Order Method
Adaptive step weak order 2.0 for Ito SDEs (deterministic order 3).
Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise.
"""
struct RI1 <: StochasticDiffEqAdaptiveAlgorithm end

"""
Rößler A., Second Order Runge–Kutta Methods for Itô Stochastic Differential Equations,
SIAM J. Numer. Anal., 47, pp. 1713-1738 (2009)
DOI:10.1137/060673308

RI3: High Weak Order Method
Adaptive step weak order 2.0 for Ito SDEs (deterministic order 3).
Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise.
"""
struct RI3 <: StochasticDiffEqAdaptiveAlgorithm end

"""
Rößler A., Second Order Runge–Kutta Methods for Itô Stochastic Differential Equations,
SIAM J. Numer. Anal., 47, pp. 1713-1738 (2009)
DOI:10.1137/060673308

RI5: High Weak Order Method
Adaptive step weak order 2.0 for Ito SDEs (deterministic order 3).
Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise.
"""
struct RI5 <: StochasticDiffEqAdaptiveAlgorithm end

"""
Rößler A., Second Order Runge–Kutta Methods for Itô Stochastic Differential Equations,
SIAM J. Numer. Anal., 47, pp. 1713-1738 (2009)
DOI:10.1137/060673308

RI6: High Weak Order Method
Adaptive step weak order 2.0 for Ito SDEs (deterministic order 2).
Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise.
"""
struct RI6 <: StochasticDiffEqAdaptiveAlgorithm end

"""
Debrabant, K. and Rößler A., Classification of Stochastic Runge–Kutta Methods for
the Weak Approximation of Stochastic Differential Equations,
Mathematics and Computers in Simulation 77, pp. 408-420 (2008)
DOI:10.1016/j.matcom.2007.04.016

RDI1WM: High Weak Order Method

Fixed step weak order 1.0 for Ito SDEs (deterministic order 2).
Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise.
"""
struct RDI1WM <: StochasticDiffEqAlgorithm end

"""
Debrabant, K. and Rößler A., Classification of Stochastic Runge–Kutta Methods for
the Weak Approximation of Stochastic Differential Equations,
Mathematics and Computers in Simulation 77, pp. 408-420 (2008)
DOI:10.1016/j.matcom.2007.04.016

RDI2WM: High Weak Order Method
Adaptive step weak order 2.0 for Ito SDEs (deterministic order 2).
Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise.
"""
struct RDI2WM <: StochasticDiffEqAdaptiveAlgorithm end

"""
Debrabant, K. and Rößler A., Classification of Stochastic Runge–Kutta Methods for
the Weak Approximation of Stochastic Differential Equations,
Mathematics and Computers in Simulation 77, pp. 408-420 (2008)
DOI:10.1016/j.matcom.2007.04.016

RDI3WM: High Weak Order Method
Adaptive step weak order 2.0 for Ito SDEs (deterministic order 3).
Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise.†
"""
struct RDI3WM <: StochasticDiffEqAdaptiveAlgorithm end

"""
Debrabant, K. and Rößler A., Classification of Stochastic Runge–Kutta Methods for
the Weak Approximation of Stochastic Differential Equations,
Mathematics and Computers in Simulation 77, pp. 408-420 (2008)
DOI:10.1016/j.matcom.2007.04.016

RDI4WM: High Weak Order Method
Adaptive step weak order 2.0 for Ito SDEs (deterministic order 3).
Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise.
"""
struct RDI4WM <: StochasticDiffEqAdaptiveAlgorithm end

# Stratonovich sense

"""
Rößler A., Second order Runge–Kutta methods for Stratonovich stochastic differential
equations, BIT Numerical Mathematics 47, pp. 657-680 (2007)
DOI:10.1007/s10543-007-0130-3

RS1: High Weak Order Method
Fixed step weak order 2.0 for Stratonovich SDEs (deterministic order 2).
Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise.
"""
struct RS1 <: StochasticDiffEqAlgorithm end

"""
Rößler A., Second order Runge–Kutta methods for Stratonovich stochastic differential
equations, BIT Numerical Mathematics 47, pp. 657-680 (2007)
DOI:10.1007/s10543-007-0130-3

RS2: High Weak Order Method
Fixed step weak order 2.0 for Stratonovich SDEs (deterministic order 3).
Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise.
"""
struct RS2 <: StochasticDiffEqAlgorithm end

"""
Kloeden, P.E., Platen, E., Numerical Solution of Stochastic Differential Equations.
Springer. Berlin Heidelberg (2011)

PL1WM: High Weak Order Method
Fixed step weak order 2.0 for Ito SDEs (deterministic order 2).
Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise.
"""
struct PL1WM <: StochasticDiffEqAlgorithm end

"""
Kloeden, P.E., Platen, E., Numerical Solution of Stochastic Differential Equations.
Springer. Berlin Heidelberg (2011)

PL1WMA: High Weak Order Method
Fixed step weak order 2.0 for Ito SDEs (deterministic order 2).
Can handle additive noise.
"""
struct PL1WMA <: StochasticDiffEqAlgorithm end

"""
Komori, Y., Weak second-order stochastic Runge–Kutta methods for non-commutative
stochastic differential equations, Journal of Computational and Applied
Mathematics 206, pp. 158 – 173 (2007)
DOI:10.1016/j.cam.2006.06.006

NON: High Weak Order Method
Fixed step weak order 2.0 for Stratonovich SDEs (deterministic order 4).
Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise.
"""
struct NON <: StochasticDiffEqAlgorithm end

"""
Komori, Y., Weak order stochastic Runge–Kutta methods for commutative stochastic
differential equations, Journal of Computational and Applied Mathematics 203,
pp. 57 – 79 (2007)
DOI:10.1016/j.cam.2006.03.010
"""
struct COM <: StochasticDiffEqAlgorithm end

"""
Komori, Y., & Burrage, K. (2011). Supplement: Efficient weak second order stochastic
Runge–Kutta methods for non-commutative Stratonovich stochastic differential equations.
Journal of computational and applied mathematics, 235(17), pp. 5326 - 5329 (2011)
DOI:10.1016/j.cam.2011.04.021
"""
struct NON2 <: StochasticDiffEqAlgorithm end

"""
Tocino, A. and Vigo-Aguiar, J., Weak Second Order Conditions for Stochastic Runge-
Kutta Methods, SIAM Journal on Scientific Computing 24, pp. 507 - 523 (2002)
DOI:10.1137/S1064827501387814

SIEA:High Weak Order Method
Fixed step weak order 2.0 for Ito SDEs (deterministic order 2).
Can handle diagonal and scalar additive noise.†
Stochastic generalization of the improved Euler method.
"""
struct SIEA <: StochasticDiffEqAlgorithm end

"""
Tocino, A. and Vigo-Aguiar, J., Weak Second Order Conditions for Stochastic Runge-
Kutta Methods, SIAM Journal on Scientific Computing 24, pp. 507 - 523 (2002)
DOI:10.1137/S1064827501387814

SMEA: High Weak Order Method
Fixed step weak order 2.0 for Ito SDEs (deterministic order 2).
Can handle diagonal and scalar additive noise.†
Stochastic generalization of the modified Euler method.
"""
struct SMEA <: StochasticDiffEqAlgorithm end

"""
Tocino, A. and Vigo-Aguiar, J., Weak Second Order Conditions for Stochastic Runge-
Kutta Methods, SIAM Journal on Scientific Computing 24, pp. 507 - 523 (2002)
DOI:10.1137/S1064827501387814

SIEB: High Weak Order Method
Fixed step weak order 2.0 for Ito SDEs (deterministic order 2).
Can handle diagonal and scalar additive noise.†
Stochastic generalization of the improved Euler method.
"""
struct SIEB <: StochasticDiffEqAlgorithm end

"""
Tocino, A. and Vigo-Aguiar, J., Weak Second Order Conditions for Stochastic Runge-
Kutta Methods, SIAM Journal on Scientific Computing 24, pp. 507 - 523 (2002)
DOI:10.1137/S1064827501387814

SMEB: High Order Weak Method
Fixed step weak order 2.0 for Ito SDEs (deterministic order 2).
Can handle diagonal and scalar additive noise.†
Stochastic generalization of the modified Euler method.
"""
struct SMEB <: StochasticDiffEqAlgorithm end

################################################################################

# IIF

struct IIF1M{F} <: StochasticDiffEqAlgorithm
  nlsolve::F
end
IIF1M(;nlsolve=NLSOLVEJL_SETUP()) = IIF1M{typeof(nlsolve)}(nlsolve)

struct IIF2M{F} <: StochasticDiffEqAlgorithm
  nlsolve::F
end
IIF2M(;nlsolve=NLSOLVEJL_SETUP()) = IIF2M{typeof(nlsolve)}(nlsolve)

struct IIF1Mil{F} <: StochasticDiffEqAlgorithm
  nlsolve::F
end
IIF1Mil(;nlsolve=NLSOLVEJL_SETUP()) = IIF1Mil{typeof(nlsolve)}(nlsolve)

################################################################################

# SDIRK
"""
ImplicitEM: Stiff Method
An order 0.5 Ito drift-implicit method.
This is a theta method which defaults to theta=1 or the Trapezoid method on the drift term.
This method defaults to symplectic=false, but when true and theta=1/2 this is the implicit Midpoint method on the drift term and is symplectic in distribution.
Can handle all forms of noise, including non-diagonal, scalar, and colored noise. Uses a 1.0/1.5 heuristic for adaptive time stepping.
"""
struct ImplicitEM{CS,AD,F,F2,FDT,ST,T2,Controller} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,FDT,ST,Controller}
  linsolve::F
  nlsolve::F2
  theta::T2
  extrapolant::Symbol
  new_jac_conv_bound::T2
  symplectic::Bool
end
ImplicitEM(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                          standardtag = Val{true}(),
                          linsolve=nothing,nlsolve=NLNewton(),
                          extrapolant=:constant,
                          theta = 1,symplectic=false,
                          new_jac_conv_bound = 1e-3,
                          controller = :Predictive) =
                          ImplicitEM{chunk_size,autodiff,
                          typeof(linsolve),typeof(nlsolve),diff_type,
                          OrdinaryDiffEq._unwrap_val(standardtag),
                          typeof(new_jac_conv_bound),controller}(
                          linsolve,nlsolve,
                          symplectic ? 1/2 : theta,
                          extrapolant,new_jac_conv_bound,symplectic)

STrapezoid(;kwargs...) = ImplicitEM(;theta=1/2,kwargs...)
SImplicitMidpoint(;kwargs...) = ImplicitEM(;theta=1/2,symplectic=true,kwargs...)
"""
ImplicitEulerHeun: Stiff Method
An order 0.5 Stratonovich drift-implicit method.
This is a theta method which defaults to theta=1/2 or the Trapezoid method on the drift term.
This method defaults to symplectic=false, but when true and theta=1 this is the implicit Midpoint method on the drift term and is symplectic in distribution.
Can handle all forms of noise, including non-diagonal, scalar, and colored noise. Uses a 1.0/1.5 heuristic for adaptive time stepping.
"""
struct ImplicitEulerHeun{CS,AD,F,FDT,ST,N,T2,Controller} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,FDT,ST,Controller}
  linsolve::F
  nlsolve::N
  theta::T2
  extrapolant::Symbol
  new_jac_conv_bound::T2
  symplectic::Bool
end
ImplicitEulerHeun(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                          standardtag = Val{true}(),
                          linsolve=nothing,nlsolve=NLNewton(),
                          extrapolant=:constant,
                          theta = 1,symplectic = false,
                          new_jac_conv_bound = 1e-3,
                          controller = :Predictive) =
                          ImplicitEulerHeun{chunk_size,autodiff,
                          typeof(linsolve),diff_type,
                          OrdinaryDiffEq._unwrap_val(standardtag),
                          typeof(nlsolve),
                          typeof(new_jac_conv_bound),controller}(
                          linsolve,nlsolve,
                          symplectic ? 1/2 : theta,
                          extrapolant,
                          new_jac_conv_bound,symplectic)

"""
ImplicitRKMil: Stiff Method
An order 1.0 drift-implicit method.
This is a theta method which defaults to theta=1 or the Trapezoid method on the drift term.
Defaults to solving the Ito problem, but ImplicitRKMil(interpretation=:Stratonovich) makes it solve the Stratonovich problem.
This method defaults to symplectic=false, but when true and theta=1/2 this is the implicit Midpoint method on the drift term and is symplectic in distribution.
Handles diagonal and scalar noise. Uses a 1.5/2.0 heuristic for adaptive time stepping.
"""
struct ImplicitRKMil{CS,AD,F,FDT,ST,N,T2,Controller,interpretation} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,FDT,ST,Controller}
  linsolve::F
  nlsolve::N
  theta::T2
  extrapolant::Symbol
  new_jac_conv_bound::T2
  symplectic::Bool
end
ImplicitRKMil(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                          standardtag = Val{true}(),
                          linsolve=nothing,nlsolve=NLNewton(),
                          extrapolant=:constant,
                          theta = 1,symplectic = false,
                          new_jac_conv_bound = 1e-3,
                          controller = :Predictive,interpretation=:Ito) =
                          ImplicitRKMil{chunk_size,autodiff,
                          typeof(linsolve),diff_type,
                          OrdinaryDiffEq._unwrap_val(standardtag),
                          typeof(nlsolve),typeof(new_jac_conv_bound),
                          controller,interpretation}(
                          linsolve,nlsolve,
                          symplectic ? 1/2 : theta,
                          extrapolant,
                          new_jac_conv_bound,symplectic)
"""
ISSEM: Stiff Method
An order 0.5 split-step Ito implicit method.
It is fully implicit, meaning it can handle stiffness in the noise term.
This is a theta method which defaults to theta=1 or the Trapezoid method on the drift term.
This method defaults to symplectic=false, but when true and theta=1/2 this is the implicit Midpoint method on the drift term and is symplectic in distribution.
Can handle all forms of noise, including non-diagonal, scalar, and colored noise. Uses a 1.0/1.5 heuristic for adaptive time stepping.
"""
struct ISSEM{CS,AD,F,FDT,ST,N,T2,Controller} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,FDT,ST,Controller}
  linsolve::F
  nlsolve::N
  theta::T2
  extrapolant::Symbol
  new_jac_conv_bound::T2
  symplectic::Bool
end
ISSEM(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                       standardtag = Val{true}(),
                       linsolve=nothing,nlsolve=NLNewton(),
                       extrapolant=:constant,
                       theta = 1,symplectic=false,
                       new_jac_conv_bound = 1e-3,
                       controller = :Predictive) =
                       ISSEM{chunk_size,autodiff,
                       typeof(linsolve),diff_type,
                       OrdinaryDiffEq._unwrap_val(standardtag),
                       typeof(nlsolve),
                       typeof(new_jac_conv_bound),controller}(
                       linsolve,nlsolve,
                       symplectic ? 1/2 : theta,
                       extrapolant,
                       new_jac_conv_bound,symplectic)
"""
ISSEulerHeun: Stiff Method
An order 0.5 split-step Stratonovich implicit method.
It is fully implicit, meaning it can handle stiffness in the noise term.
This is a theta method which defaults to theta=1 or the Trapezoid method on the drift term.
This method defaults to symplectic=false, but when true and theta=1/2 this is the implicit Midpoint method on the drift term and is symplectic in distribution.
Can handle all forms of noise, including non-diagonal,Q scalar, and colored noise. Uses a 1.0/1.5 heuristic for adaptive time stepping.
"""
struct ISSEulerHeun{CS,AD,F,FDT,ST,N,T2,Controller} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,FDT,ST,Controller}
 linsolve::F
 nlsolve::N
 theta::T2
 extrapolant::Symbol
 new_jac_conv_bound::T2
 symplectic::Bool
end
ISSEulerHeun(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                      standardtag = Val{true}(),
                      linsolve=nothing,nlsolve=NLNewton(),
                      extrapolant=:constant,
                      theta = 1,symplectic=false,
                      new_jac_conv_bound = 1e-3,
                      controller = :Predictive) =
                      ISSEulerHeun{chunk_size,autodiff,
                      typeof(linsolve),diff_type,
                      OrdinaryDiffEq._unwrap_val(standardtag),
                      typeof(nlsolve),typeof(new_jac_conv_bound),controller}(
                      linsolve,nlsolve,
                      symplectic ? 1/2 : theta,
                      extrapolant,
                      new_jac_conv_bound,symplectic)
"""
SKenCarp: Stiff Method
Adaptive L-stable drift-implicit strong order 1.5 for additive Ito and Stratonovich SDEs with weak order 2.
Can handle diagonal, non-diagonal and scalar additive noise.
"""
struct SKenCarp{CS,AD,F,FDT,ST,N,T2,Controller} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,FDT,ST,Controller}
  linsolve::F
  nlsolve::N
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
  ode_error_est::Bool
end

SKenCarp(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                   standardtag = Val{true}(),
                   linsolve=nothing,nlsolve=NLNewton(),
                   smooth_est=true,extrapolant=:min_correct,
                   new_jac_conv_bound = 1e-3,controller = :Predictive,
                   ode_error_est = true) =
 SKenCarp{chunk_size,autodiff,typeof(linsolve),diff_type,
        OrdinaryDiffEq._unwrap_val(standardtag),
        typeof(nlsolve),typeof(new_jac_conv_bound),controller}(
        linsolve,nlsolve,smooth_est,extrapolant,new_jac_conv_bound,
        ode_error_est)

################################################################################

# Jumps

struct TauLeaping <: StochasticDiffEqJumpAdaptiveAlgorithm end
struct CaoTauLeaping <: StochasticDiffEqJumpAdaptiveAlgorithm end

################################################################################

# Etc.

struct StochasticCompositeAlgorithm{T,F} <: StochasticDiffEqCompositeAlgorithm
  algs::T
  choice_function::F
end

struct RandomEM <: StochasticDiffEqRODEAlgorithm end

const SplitSDEAlgorithms = Union{IIF1M,IIF2M,IIF1Mil,SKenCarp,SplitEM}

@doc raw"""
Leimkuhler B., Matthews C., Robust and efficient configurational molecular sampling via
Langevin dynamics, J. Chem. Phys. 138, 174102 (2013)
DOI:10.1063/1.4802990

```math
du = vdt \\
dv = f(v,u) dt - \gamma v dt + g(u) \sqrt{2\gamma} dW
```
"""
struct BAOAB{T} <: StochasticDiffEqAlgorithm
  gamma::T
end
BAOAB(;gamma=1.0) = BAOAB(gamma)


text:
__precompile__(false)

module Process

using Distributed

workers = Future[]
channels=Dict{Symbol,Distributed.AbstractRemoteRef}()

channels[:JOBS]=RemoteChannel(()->Channel{Function}(Inf))
#channels[:TRIGGER]=RemoteChannel(()->Channel{Bool}(1))
#channels[:CAMERA]=RemoteChannel(()->Channel{Tuple{Array{Float32,1},Array{Float32,1}}}(Inf))
#channels[:BOOL]=RemoteChannel(()->Channel{Tuple{Symbol,Bool}}(Inf))
#channels[:SCENE]=RemoteChannel(()->Channel{Int}(1))
#channels[:UPLOAD]=RemoteChannel(()->Channel{Symbol}(Inf))

updateChannels(chs::Dict{Symbol,Distributed.AbstractRemoteRef}) = global channels=chs
getChannels() = channels
register(id, channel) = channels[id]=channel

function createProcess(pool)
  global channels
  chs = channels
  push!(workers, remotecall(function()
    println("LOOP START")
    while true
      #fetch(remoteChannels)  #wait until first channel
      #chs = Dict([c for c in remoteChannels])
      println("updateChannels")
      updateChannels(chs)

      has = haskey(chs, :JOBS)
      if has
        println("GET JOB")
        job=take!(chs[:JOBS])
        println("START JOB")
        job() # execute job
        println("END JOB")
      end

      sleep(has ? 0 : 1)
    end
    true
  end, pool))
end

waitToEnd() = for w in workers println(fetch(w)) end

addJob(job) = put!(channels[:JOBS],job)

#remoteChannels=RemoteChannel(()->Channel{Tuple{Symbol,Distributed.AbstractRemoteRef}}(0))
#put!(remoteChannels,(:JOBS,RemoteChannel(()->Channel{Function}(0))))
#put!(remoteChannels,(:TRIGGERUPDATE,RemoteChannel(()->Channel{Bool}(0))))
#channels = Dict([c for c in remoteChannels])

function getInput()
  println("\x1b[91m>\x1b[0m\x1b[96m ")
  input = readline()
  println(input,"\x1b[0m")
end

function inputLoop()
  input = ""
  while input != "q" && input != "quit"
    input = getInput()
  end
end

function outputLoop()
  #Main.App.run()
  while true
    println(".")
    sleep(3)
  end
end

function start(mainProcess::Function)
  p=procs()
  popfirst!(p)
  if length(p) > 0
    pool=WorkerPool(p)
    createProcess(pool)
    addJob(inputLoop)
  end
  mainProcess()
  waitToEnd()
end

#=
function thread_compute(this::ThreadManager.Thread)
  println("thread_compute")
  #timeRef = Ref(0.0)
  #i=0
  while true
    #take!(TRIGGERUPDATE)
    #if TimeManager.OnTime(0.25, timeRef)
    #  info(this, i ;title="CALC")
    #  i+=1
    #end
    #sleep(0.001)
    println("thread_compute: work")
    sleep(3)
  end
end

function thread_renderer(this::ThreadManager.Thread)
  println("thread_renderer")
  run()
end
=#

#=
function test_thread()
  id=Threads.threadid()
  i=1
  while true
    #if id==1
    #  print(".")
    #  Libc.systemsleep(0.1)
    #elseif id==2
    #  println("Compute", id)
    #  for i=1:99999999 b=i^i; end;
    #  tprintln("Compute", id, "end.")
    #  Libc.systemsleep(0.1)
    #else
      #println("Idle", id)
      #Libc.systemsleep(10)  #rand(1:10) #sleep(rand(0.1:3))
      i += 1
    #end
  end
end

# Does not work with OpenGL, because OpenGL has to be on main thread...
function start_threads()
  return false

  println("Start Threads")

  max=Threads.nthreads()
  Threads.@threads for i = 1:max
    if i == 1
      run()
    elseif i == 2
      test_thread()
    end
  end

  #pool = ThreadManager.thread_pool()

  #ThreadManager.set(pool, thread_compute, "Compute")
  #ThreadManager.set(pool, thread_renderer, "Renderer")
  #ThreadManager.set(pool, thread_sound, "Sound")

  #ThreadManager.start(pool) # anything below this line will be paused until threads are closed
end

function fetchlast!(c::Distributed.AbstractRemoteRef)
  v=nothing; while isready(c) v=take!(c); end; v
end

function waitAndUpdateChunks()
  global CAMERA, FRUSTUM_CULLING, HIDE_UNSEEN_CUBES

  println("Presets")

  cameradata=channels[:CAMERA]
  bools=channels[:BOOL]

  init()

  trigger=false

  println("start update")
  while true

    if (data=fetchlast!(cameradata)) != nothing
      CAMERA.position, CAMERA.rotation = data
      Update(CAMERA)
      #println("updated Camera")
    end

    while isready(bools)
      v=take!(bools)
      id, value = v
      if id == :FRUSTUM_CULLING FRUSTUM_CULLING = value;
      elseif id == :HIDE_UNSEEN_CUBES HIDE_UNSEEN_CUBES = value;
      end
      trigger = true
    end

    if trigger
      trigger = false
      setFrustumMode()
    end

    sleep(0.1)
  end
  nothing
end

updateChunk(this::Chunk) =   #println("send Chunk data")
   if myid() != 1 remotecall(setChunkInstances, 1, chunk_instances) end

 createChunk(this::Chunk)=

  if myid() != 1
    sz=channels[:SCENE]
    if isready(sz) SCENE=take!(sz) end
  end

function uploadChunk(s)
  global uploaded = s
  if myid() == 1 return end
  #println("upload Chunk")
  #@save "chunk.jld2" chunk_instances
  put!(channels[:UPLOAD],s)
end

function loadChunk()
  if length(procs()) <= 1 return uploaded end

  upload=channels[:UPLOAD]
  if !isready(upload) return :NOTHING end
  #println("load Chunk")
  m = take!(upload)
  #@load "chunk.jld2" chunk_instances
  m
end

uploadData()
  m=loadChunk()
  if m == :NOTHING return end

  #channels["KEYS"]
  #put!(channels[:BOOL], (:NOTHING, false))
  put!(channels[:BOOL], (:FRUSTUM_CULLING, FRUSTUM_CULLING))
  put!(channels[:BOOL], (:HIDE_UNSEEN_CUBES, HIDE_UNSEEN_CUBES))

checkCamera() =
      if length(procs()) > 1 put!(channels[:BOOL], (:NOTHING, false))
      else setFrustumMode()
      end

chooseRenderMethod(method=RENDER_METHOD) =   if myid() != 1 put!(channels[:SCENE],SCENE) end
=#

end #Process

ext: jl

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
module Day8

using StatsBase

inp = readline("day8.e") |> Vector{Char} |> x -> parse.(Int8, x)
inp = reshape(inp, 25, 6, :)

function solvea()
    maxlayer = 0
    layer = argmin(count.(==(0), eachslice(inp; dims = 3)))
    counts = countmap(inp[:, :, layer])
    counts[1] * counts[2]
end

function solveb()
    ans = zeros(Int, 25, 6)
    for i in 1:25, j in 1:6
        slice = inp[i, j, :]
        ans[i, j] = slice[findfirst(!=(2), slice)]
    end

    transpose(ans)
end

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
#=
ellipse_sino.jl
2019-07-12, Helena H
2019-07-13, Jeff Fessler, refactor to use sg.grid
=#

export ellipse_sino

#using MIRT: sino_geom, SinoGeom, downsample2

"""
`sino = ellipse_sino(sg, ells ; oversample=1, xscale=1, yscale=1)`

Create sinogram projection of one or more ellipses.
Works for any sinogram geometry.

in
- `sg::SinoGeom` sinogram geometry object from `sino_geom()`
- `ells::Matrix` `[ne 6]` ellipse parameters
   * `[centx centy radx rady angle_degrees amplitude]`

options
- `oversample::Int`	oversampling factor for emulating "strips"
   * default 1: just 1 ray per detector element
- `xscale::Int`		use -1 to flip in x (not recommended); default 1
- `yscale::Int`		use -1 to flip in y (not recommended); default 1

out
- `sino` `[nb na]` sinogram

To get the sample locations, use `(pos,ang) = sg.grid`
"""
function ellipse_sino(
    sg::SinoGeom,
    ells::AbstractMatrix{<:Real} ;
    oversample::Int = 1,
    kwargs...,
)

	sg = sg.over(oversample)
	(rg, ϕg) = sg.grid
	sino = ellipse_sino(rg, ϕg, ells ; kwargs...)
	if oversample > 1
		sino = downsample2(sino, (oversample, 1))
	end
	return sino
end

"""
`sino = ellipse_sino(rg::AbstractArray{<:Real}, ϕg::AbstractArray{<:Real},
		ells ; oversample=1, xscale=1, yscale=1)`

Create sinogram projection of one or more ellipses,
for arbitrary radial/angular sampling grid locations `(rg, ϕg)`

in
- `rg::AbstractArray{<:Real}`	radial sampling locations
- `ϕg::AbstractArray{<:Real}`	angular sampling locations (radians)
- `ells::Matrix` `[ne 6]` ellipse parameters
   * `[centx centy radx rady angle_degrees amplitude]`

options
- `xscale::Int` use -1 to flip in x (not recommended); default 1
- `yscale::Int` use -1 to flip in y (not recommended); default 1

out
- `sino::AbstractArray{Float32}` same size as `rg` and `ϕg`
"""
function ellipse_sino(
    rg::AbstractArray{<:Real},
    ϕg::AbstractArray{<:Real},
    ells::AbstractMatrix{<:Real} ;
    xscale::Int = 1,
    yscale::Int = 1,
)

	size(ells,2) != 6 && throw("6 parameters per ellipse")
	size(rg) != size(ϕg) && throw("rg and ϕg size mismatch")

	sino = zeros(Float32, size(rg))

	cangs = cos.(ϕg)
	sangs = sin.(ϕg)

	#loop over ellipses
	#ticker reset
	ne = size(ells, 1)
	for ie in 1:ne
		#ticker(mfilename, ie, ne)
		ell = ells[ie, :]

		cx = ell[1] * xscale
		cy = ell[2] * yscale
		rx = ell[3]
		ry = ell[4]
		(rx <= 0) || (ry <= 0) && throw("need positive radii")
		eang = deg2rad(ell[5])
		val = ell[6]

		if yscale == -1
			eang = -eang
		end
		if xscale == -1
			eang = pi - eang
		end
		scale = 2 * val * rx * ry

		# square of projected radius:
		rp2 = @. (rx * (cangs * cos(eang) + sangs * sin(eang)))^2 +
				(ry * (sangs * cos(eang) - cangs * sin(eang)))^2
		sp = cx * cangs + cy * sangs # radial shift
		dis2 = (rg - sp).^2 # square of distances from center
		@. sino += scale / rp2 * sqrt(max(rp2 - dis2, 0))
	end

	return sino
end


text:
const wall_rebound = [
    v -> v[1] ==  1 ? - rot * v :   rot * v,
    v -> v[2] ==  1 ? - rot * v :   rot * v,
    v -> v[2] == -1 ?   rot * v : - rot * v,
    v -> v[2] == -1 ? - rot * v :   rot * v ]

function step!(sq :: SquareParticles, pc, bc)

    n = sq.n
    ϵ = sq.ϵ
    tempsp, i1, i2 = dt_min_position(pc)
    tempsm, j1, j2 = dt_min_position(bc)

    if tempsp < tempsm

        for i in 1:n
            sq.q[i] = sq.q[i] + tempsp * sq.v[i]
        end

        if sq.v[i1]'sq.v[i2] == 0

            sq.v[[i1, i2]] = sq.v[[i2, i1]]  # swap velocities

        elseif sq.v[i1]'sq.v[i2] == -1

            rot = [0 -1; 1 0]

            if (rot * sq.v[i1])'*(sq.q[i2] - sq.q[i1]) < 0
                sq.v[i1] = rot * sq.v[i1]
                sq.v[i2] = rot * sq.v[i2]
            else
                sq.v[i1] = -rot * sq.v[i1]
                sq.v[i2] = -rot * sq.v[i2]
            end
        end

        pc.dt .= pc.dt .- tempsp
        pc.dt[:, i1] .= Inf
        pc.dt[:, i2] .= Inf

        bc.dt .= bc.dt .- tempsp
        bc.dt[i1,:] .= Inf
        bc.dt[i2,:] .= Inf

        compute_dt!(pc, i1, sq)
        compute_dt!(pc, i2, sq)

        t, i = compute_dt(sq, i1)
        bc.dt[i1, i] = t

        t, i = compute_dt(sq, i2)
        bc.dt[i2, i] = t

    else

        for i in 1:n
            sq.q[i] = sq.q[i] + tempsm * sq.v[i]
        end

        sq.v[j1] = wall_rebound[j2](sq.v[j1])

        pc.dt .-= tempsm
        reset!(pc, j1)

        bc.dt .-= tempsm
        bc.dt[j1,:] .= Inf

        compute_dt!(pc, j1, sq)

        t, i = compute_dt(sq, j1)
        bc.dt[j1, i] = t

    end

    return min(tempsp, tempsm)

end


text:
module StructDualDynProg

using StochOptInterface
const SOI = StochOptInterface

# Generic implementation of Stochastic Program supporting the conversion from a StructJuMP model to a Stochastic Program
include("StructProg/StructProg.jl")
export StructProg

# SDDP algorithm on top of these abstract components
include("SDDP/SDDP.jl")
export SDDP

# Wait and see value computation for StructProg
include("WaitAndSee/waitandsee.jl")
export WaitAndSee

end # module


text:
# load an example data set
function example_data(name)
    basepath = @artifact_str name
    paths = readdir(basepath)
    if length(paths) == 1
        path = joinpath(basepath, paths[1])
        data, header = readdlm(path, ',', header = true)
        data_raw = DataFrame(data, vec(header))
        for col in eachcol(data_raw)
            replace!(col,"NA" => missing)
        end
        data = identity.(data_raw)
        return data
    else
        out = Dict()
        for path in paths
            key = Symbol(split(path, ".")[1])
            path = joinpath(basepath, path)
            data, header = readdlm(path, ',', header = true)
            data_raw = DataFrame(data, vec(header))
            for col in eachcol(data_raw)
                replace!(col,"NA" => missing)
            end
            data = identity.(data_raw)
            push!(out, key => data)
        end
        return out
    end
end


text:
## LICENSE
##   Copyright (c) 2018 GEOPHYSTECH LLC
##
##   Licensed under the Apache License, Version 2.0 (the "License");
##   you may not use this file except in compliance with the License.
##   You may obtain a copy of the License at
##
##       http://www.apache.org/licenses/LICENSE-2.0
##
##   Unless required by applicable law or agreed to in writing, software
##   distributed under the License is distributed on an "AS IS" BASIS,
##   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##   See the License for the specific language governing permissions and
##   limitations under the License.

## initial release by Andrey Stepnov, email: a.stepnov@geophsytech.ru

## read VS30 file
function read_vs30_file(filename::String)
  A = readdlm(filename)
  vs30_row_num = length(A[:,1])
  B = Array{Point_vs30}(undef,0)
  for i=1:vs30_row_num
    B = push!(B,Point_vs30(A[i,1],A[i,2],A[i,3]))
  end
  return B
end
## convert Array{Float64,2} array to Array{GroundMotion.Point_vs30,1}
function convert_to_point_vs30(A::Array{Float64,2})
  vs30_row_num = length(A[:,1])
  B = Array{Point_vs30}(undef,0)
  for i=1:vs30_row_num
    B = push!(B,Point_vs30(A[i,1],A[i,2],A[i,3]))
  end
  return B
end
## read VS30 file with Dl column
function read_vs30_dl_file(filename::String)
  A = readdlm(filename)
  vs30_row_num = length(A[:,1])
  B = Array{Point_vs30_dl}(undef,0)
  for i=1:vs30_row_num
    B = push!(B,Point_vs30_dl(A[i,1],A[i,2],A[i,3],A[i,4]))
  end
  return B
end
## convert Array{Point_pga_out}(N,1) array to Array{Float64}(N,4)
"""
convert_to_float_array(B::Array{T,N})
  where T is Point_{pga,pgv,pgd}_out, Point_vs30

  Convert 1-d array of custom type to Array{Float64}(N,X)
"""
function convert_to_float_array(B::Array{Point_pga_out})
  A = Array{Float64}(undef,length(B),3)
  for i=1:length(B)
   A[i,1] = B[i].lon
   A[i,2] = B[i].lat
   A[i,3] = B[i].pga
  end
  return A
end
function convert_to_float_array(B::Array{Point_vs30})
  A = Array{Float64}(undef,length(B),3)
  for i=1:length(B)
   A[i,1] = B[i].lon
   A[i,2] = B[i].lat
   A[i,3] = B[i].vs30
  end
return A
end
function convert_to_float_array(B::Array{Point_vs30_dl})
  A = Array{Float64}(undef,length(B),4)
  for i=1:length(B)
   A[i,1] = B[i].lon
   A[i,2] = B[i].lat
   A[i,3] = B[i].vs30
   A[i,4] = B[i].dl
  end
return A
end


text:
include("VecUtils.jl")

struct NURBS
    number::Array{Int64,1}
    coefs::Array{Float64}
    knots::Array{Array{Float64,1}}
    order::Array{Int64,1}
end

function nrbmak(coefs, knots)
    np = [i for i in size(coefs)]
    dim = np[1]

    if length(knots)==3
        if length(np)==3
            np = push!(np,1)
        elseif length(np)==2
            np = append!(np,[1,1])
        end
        #constructing a volume
        number = np[2:4]
        if dim<4
            temp_coefs = coefs
            coefs = repeat([0., 0., 0., 1.], 1, np[2], np[3], np[4])
            coefs[1:dim,:,:,:] = temp_coefs
        end
        uorder = length(knots[1])-np[2]
        vorder = length(knots[2])-np[3]
        worder = length(knots[3])-np[4]
        uknots = sort(knots[1])
        vknots = sort(knots[2])
        wknots = sort(knots[3])
        knots = [uknots, vknots, wknots]
        order = [uorder, vorder, worder]
    elseif length(knots)==2
        if length(np)==2
            np = push!(np,1)
        end
        #constructing a surface
        number = np[2:3]
        if dim<4
            temp_coefs = coefs
            coefs = repeat([0., 0., 0., 1.], 1, np[2], np[3])
            coefs[1:dim,:,:,:] = temp_coefs
        end
        uorder = length(knots[1])-np[2]
        vorder = length(knots[2])-np[3]
        uknots = sort(knots[1])
        vknots = sort(knots[2])
        knots = [uknots, vknots]
        order = [uorder, vorder]
    else
        #constructing a curve
        if length(np)==1
            number = [np[1]]
        else
            number = [np[2]]
        end
        if dim<4
            temp_coefs = coefs
            coefs = repeat([0., 0., 0., 1.], 1, number[1])
            coefs[1:dim,:,:,:] = temp_coefs
        end
        order = [length(knots[1])-number[1]]
        uknots = sort(knots[1])
        knots = [uknots]
    end

    nrb = NURBS(number, coefs, knots, order)
    return nrb
end

function findspan(n::Int64, u, knot::Array{Float64,1})
    if (minimum(u)<knot[1]) || (maximum(u)>knot[end])
        error("Some value is outside the knot span")
    end
    if isa(u, Number)
        s = [0]
    else
        s = similar(u, Int64)
    end
    for j=1:length(u)
        if u[j]==knot[n+2]
            s[j]=n
            continue
        end
        s[j]=findlast(knot.<=u[j])-1
    end
    return s
end

function basisfun(iv::Array{Int64}, uv, p::Integer, U)
    B = zeros(length(uv), p+1)
    N = zeros(p+1)
    for jj = 1:length(uv)
        i = iv[jj]+1 #findspan uses 0-based numbering
        u = uv[jj]
        left = zeros(p+1)
        right = zeros(p+1)
        N[1] = 1

        for j=1:p
            left[j+1] = u - U[i+1-j]
            right[j+1] = U[i+j] - u;
            saved = 0
            for r = 0:j-1
                temp = N[r+1]/(right[r+2]+left[j-r+1])
                N[r+1] = saved + right[r+2]*temp
                saved = left[j-r+1]*temp
            end
            N[j+1] = saved
        end
        B[jj,:] = N
    end
    return B
end

function bspeval(d, c, k, u)
    nu = length(u)
    mc, nc = size(c)
    #readline(stdin)
    s = findspan(nc-1, u, k)
    N = basisfun(s, u, d, k)
    tmp1 = s .- d .+ 1
    p = zeros(mc, nu)
    for i=0:d
        p = p + repeat(N[:,i+1]',mc,1).*c[:,tmp1.+i]
    end
    return p
end

function nrbeval(nurbs::NURBS, tt::AbstractArray)
    foption = true
    degree = nurbs.order .- 1
    if length(nurbs.knots)==3
        #NURBS structure represents a volume
        num1 = nurbs.number[1]
        num2 = nurbs.number[2]
        num3 = nurbs.number[3]

        if isa(tt[1], Array)
            nt1 = length(tt[1])
            nt2 = length(tt[2])
            nt3 = length(tt[3])

            #Evaluate along the w direction
            val = reshape(nurbs.coefs, 4*num1*num2, num3)
            val = bspeval(degree[3], val, nurbs.knots[3], tt[3])
            val = reshape(val, (4, num1, num2, nt3))

            #Evaluate along the v direction
            val = permutedims(val, [1, 2, 4, 3])
            val = reshape(val, 4*num1*nt3, num2)
            val = bspeval(degree[2], val, nurbs.knots[2], tt[2])
            val = reshape(val, 4, num1, nt3, nt2)
            val = permutedims(val, [1, 2, 4, 3])

            #Evaluate along the u direction
            val = permutedims(val, [1, 3, 4, 2])
            val = reshape(val, 4*nt2*nt3, num1)
            val = bspeval(degree[1], val, nurbs.knots[1], tt[1])
            val = reshape(val, 4, nt2, nt3, nt1)
            val = permutedims(val, [1, 4, 2, 3])
            pnts = val
            p = pnts[1:3,:,:,:]
            w = pnts[4,:,:,:,:]'
            if foption
                p = p./repeat(w, 3, 1, 1, 1)
            end
        else
            #Evaluate at scattered points
            #tt[1,:] represents the u direction
            #tt[2,:] represents the v direction
            #tt[3,:] represents the w direction
            st = size(tt)
            if (st[1]!=3 && st[2]==3 && length(st)==2)
                tt = tt'
                st = size(tt)
            end
            nt = prod(st[2:end])
            tt = reshape(tt, (3, nt))
            # Evaluate along the w direction
            val = reshape(nurbs.coefs, 4*num1*num2, num3)
            val = bspeval(degree[3], val, nurbs.knots[3], tt[3,:])
            val = reshape(val, (4, num1, num2, nt))
            # Evaluate along the v direction
            val2 = zeros(4*num1, nt)
            for v=1:nt
                coefs = reshape(val[:,:,:,v], (4*num1, num2))
                val2[:,v] = bspeval(degree[2], coefs, nurbs.knots[2], tt[2,v])
            end
            val2 = reshape(val2, (4, num1, nt))
            # Evaluate along the u direction
            pnts = zeros(4, nt)
            for v = 1:nt
                coefs = reshape(val2[:,:,v], (4, num1))
                pnts[:,v] = bspeval(degree[1], coefs, nurbs.knots[1], tt[1,v])
            end
            w = pnts[4,:]'
            p = pnts[1:3,:]
            if (foption)
                p = p./repeat(w,3,1)
            end

            if length(st)!=2
                w = reshape(w, st[2:end])
                p = reshape(p, tuple((3,)...,st[2:end]...))
            end
        end
    elseif length(nurbs.knots)==2
        # NURBS structure represents a surface
        num1 = nurbs.number[1]
        num2 = nurbs.number[2]
        if isa(tt[1],Array)
            #Evaluate over a [u,v] grid
            # tt[1] represents the u direction
            # tt[2] represents the v direction
            nt1 = length(tt[1])
            nt2 = length(tt[2])

            #Evaluate along the v direction
            val = reshape(nurbs.coefs, 4*num1, num2)
            val = bspeval(degree[2], val, nurbs.knots[2], tt[2])
            val = reshape(val, 4, num1, nt2)

            #Evaluate along the u direction
            val = permutedims(val, [1, 3, 2])
            val = reshape(val, 4*nt2, num1)
            val = bspeval(degree[1], val, nurbs.knots[1], tt[1])
            val = reshape(val, 4, nt2, nt1)
            val = permutedims(val, [1, 3, 2])

            w = val[[4],:,:]
            p = val[1:3, :, :]
            if (foption)
                p = p./repeat(w, 3, 1, 1)
            end
        else
            # Evaluate at scattered points
            # tt[1,:] represents the u direction
            # tt[2,:] represents the v direction
            st = size(tt)
            if (st[1]!=2 && st[2]==2 && length(st)==2)
                tt = tt'
                st = size(tt)
            end
            nt = prod(st[2:end])
            tt = reshape(tt, (3, nt))
            # Evaluate along the v direction
            val = reshape(nurbs.coefs, 4*num1, num2)
            val = bspeval(degree[2], val, nurbs.knots[2], tt[2,:])
            val = reshape(val, (4, num1, nt))
            #Evaluate along the u direction
            pnts = zeros(4, nt)
            for v = 1:nt
                coefs = reshape(val[:,:,v], (4, num1))
                pnts[:,v] = bspeval(degree[1], coefs, nurbs.knots[1], tt[1,v])
            end
            w = pnts[4,:]'
            p = pnts[1:3,:]
            if (foption)
                p = p./repeat(w,3,1)
            end
            if (length(st)!=2)
                w = reshape(w, st[2:end])
                p = reshape(p, tuple((3,)...,st[2:end]...))
            end
        end
    else
        if isa(tt[1], Array) && length(tt)==1
            tt =tt[1]
        end
        st = size(tt)
        val = bspeval(nurbs.order[1]-1, nurbs.coefs, nurbs.knots[1], [tt...])
        w = val[4,:]'
        p = val[1:3,:]
        if foption
            p = p./repeat(w,3,1)
        end
        if (st[1]!=1 || length(st)!=2)
            w = reshape(w, st)
            p = reshape(p, tuple((3,)...,st...))
        end
    end
    return p
end

function nrbextract(srf::NURBS)
    for idim = 1:length(srf.knots)
        ord = srf.order[idim]
        if (srf.knots[idim][1]!=srf.knots[idim][ord] ||
            srf.knots[idim][end]!=srf.knots[idim][end-ord+1])
            error("nrbextract: only working with open knot vectors")
        end
    end

    if (length(srf.knots)==2)
        crvs = Array{NURBS}(undef, 4)
        for ind=1:2
            ind2 = mod(ind,2)+1
            bnd1 = (ind-1)*2+1
            bnd2 = (ind-1)*2+2
            if ind==1
                coefs1 = srf.coefs[:,1,:]
                coefs2 = srf.coefs[:,end,:]
            elseif ind==2
                coefs1 = srf.coefs[:,:,1]
                coefs2 = srf.coefs[:,:,end]
            end
            crvs[bnd1] = nrbmak(coefs1, srf.knots[ind2])
            crvs[bnd2] = nrbmak(coefs2, srf.knots[ind2])
        end
    elseif (length(srf.knots)==3)
        crvs = Array{NURBS}(undef, 6)
        for ind=1:3
            inds = setdiff(1:3, ind)
            bnd1 = (ind - 1)*2 + 1
            bnd2 = (ind - 1)*2 + 2
            if ind==1
                coefs1 = srf.coefs[:,1,:,:]
                coefs2 = srf.coefs[:,end,:,:]
            elseif ind==2
                coefs1 = srf.coefs[:,:,1,:]
                coefs2 = srf.coefs[:,:,end,:]
            elseif ind==3
                coefs1 = srf.coefs[:,:,:,1]
                coefs2 = srf.coefs[:,:,:,end]
            end
            crvs[bnd1] = nrbmak(coefs1, [srf.knots[inds[1]], srf.knots[inds[2]]])
            crvs[bnd2] = nrbmak(coefs2, [srf.knots[inds[1]], srf.knots[inds[2]]])
        end
    else
        error("The entity is not a surface nor a volume")
    end
    return crvs
end

function nrbline(points...)
    numPts = length(points)
    if numPts==1
        error("More than one point should be specfified")
    end
    coefs = [zeros(3,numPts); ones(1,numPts)]
    for iPt = 1:numPts
        coefs[1:length(points[iPt]),iPt] .= points[iPt]
    end
    knots = vcat(0, LinRange(0,1,numPts), 1)
    line = nrbmak(coefs, [knots])
    return line
end

function nrbline()
    p1 = [0,0]
    p2 = [1,0]
    line = nrbline(p1, p2)
    return line
end

function nrbrect(width, height)
    coefs =[0. width width width width 0. 0. 0.;
          0. 0. 0. height height height height 0.;
          0. 0. 0. 0. 0. 0. 0. 0.;
          1. 1. 1. 1. 1. 1. 1. 1.]
    knots = [0, 0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1, 1]
    curve = nrbmak(coefs, [knots])
    return curve
end

function nrbrect(size)
    return nrbrect(size, size)
end

function nrbrect()
    return nrbrect(1.)
end

function nrbcirc(radius::Real, center::Vector, sang::Real, eang::Real)
    sweep = eang - sang
    if sweep < 0
        sweep = 2*pi + sweep
    end

    if abs(sweep) <= pi/2
        narcs = 1   # number of arc segments
        knots = [0, 0, 0, 1, 1, 1]
    elseif abs(sweep) <= pi
        narcs = 2
        knots = [0, 0, 0, 0.5, 0.5, 1, 1, 1]
    elseif abs(sweep) <= 3*pi/2
        narcs = 3
        knots = [0, 0, 0, 1/3, 1/3, 2/3, 2/3, 1, 1, 1]
    else
        narcs = 4
        knots = [0, 0, 0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1, 1, 1]
    end

    dsweep = sweep/(2*narcs)  #arc segment sweep angle/2

    # determine middle control point and weight
    wm = cos(dsweep)
    x  = radius*wm
    y  = radius*sin(dsweep)
    xm = x+y*tan(dsweep)

    # arc segment control points
    ctrlpt = [ x wm*xm x;    # w*x - coordinate
              -y 0     y;    # w*y - coordinate
               0 0     0;    # w*z - coordinate
               1 wm    1]    # w   - coordinate

    # build up complete arc from rotated segments
    coefs = zeros(4,2*narcs+1)   # nurbs control points of arc
    xx = vecrotz(sang + dsweep)
    coefs[:,1:3] = xx*ctrlpt;     # rotate to start angle
    xx = vecrotz(2*dsweep);
    for n = 2:narcs
       m = 2*n.+[0, 1]
       coefs[:,m] = xx*coefs[:,m.-2]
    end

    # vectrans arc if necessary
    if !isempty(center)
      xx = vectrans(center)
      coefs = xx*coefs
    end

    curve = nrbmak(coefs,[knots])
    return curve
end

function nrbcirc(radius::Real, center::Vector)
    sang = 0.
    eang = 2*pi
    curve = nrbcirc(radius, center, sang, eang)
    return curve
end

function nrbcirc(radius::Real)
    center = []
    curve = nrbcirc(radius, center)
    return curve
end

function nrbcirc()
    radius = 1
    curve = nrbcirc(radius)
    return curve
end

function bspkntins(d, c, k, u)
    tolEq = 1e-10
    mc, nc = size(c)
    sort!(u)
    nu = length(u)
    nk = length(k)

    ic = zeros(mc, nc+nu)
    ik = zeros(nk+nu)

    n = nc - 1
    r = nu - 1

    m = n + d + 1
    a = findspan(n, u[1], k)[1]
    b = findspan(n, u[r+1], k)[1]
    b += 1

    ic[:, 1:a-d+1] = c[:, 1:a-d+1]
    ic[:, b+nu:nc+nu] = c[:, b:nc]

    ik[1:a+1] = k[1:a+1]
    ik[b+d+nu+1:m+nu+1] = k[b+d+1:m+1]

    ii = b + d - 1
    ss = ii + nu

    for jj = r:-1:0
        ind = (a+1):ii
        ind = ind[u[jj+1].<=k[ind.+1]]
        ic[:, ind.+ss.-ii.-d] = c[:,ind.-d]
        ik[ind.+ss.-ii.+1] = k[ind.+1]
        ii = ii - length(ind)
        ss = ss - length(ind)

        ic[:,ss-d] = ic[:, ss-d+1]
        for l=1:d
            ind = ss - d + l
            alfa = ik[ss+l+1] - u[jj+1]
            if abs(alfa) < tolEq
                ic[:,ind] = ic[:, ind+1]
            else
                alfa = alfa/(ik[ss+l+1]-k[ii-d+l+1])
                tmp = (1-alfa) * ic[:, ind+1]
                ic[:,ind] = alfa*ic[:,ind] + tmp
            end
        end
        ik[ss+1] = u[jj+1]
        ss = ss-1
    end
    return ic, ik
end

function nrbkntins(nurbs::NURBS, iknots::AbstractArray)
    degree = nurbs.order .- 1
    fmax(x,y) = any(y.>maximum(x))
    fmin(x,y) = any(y.<minimum(x))
    for i=1:length(nurbs.knots)
        if any(fmax(nurbs.knots[i], iknots[i])) || any(fmin(nurbs.knots[i], iknots[i]))
            error("Trying to insert a knot outside the interval of definition")
        end
    end
    knots = Array{Array}(undef, length(nurbs.knots))
    if length(nurbs.knots)==3
        # NURBS represents a volume
        num1 = nurbs.number[1]
        num2 = nurbs.number[2]
        num3 = nurbs.number[3]

        # Insert knots along the w direction
        if isempty(iknots[3])
            coefs = nurbs.coefs
            knots[3] = nurbs.knots[3]
        else
            coefs = reshape(nurbs.coefs, 4*num1*num2, num3)
            coefs, knots[3] = bspkntins(degree[3], coefs, nurbs.knots[3], iknots[3])
            num3 = size(coefs, 2)
            coefs = reshape(coefs, 4, num1, num2, num3)
        end

        # Insert knots along the v direction
        if isempty(iknots[2])
            knots[2] = nurbs.knots[2]
        else
            coefs = permutedims(coefs, [1, 2, 4, 3])
            coefs = reshape(coefs, 4*num1*num3, num2)
            coefs, knots[2] = bspkntins(degree[2], coefs, nurbs.knots[2], iknots[2])
            num2 = size(coefs, 2)
            coefs = reshape(coefs, 4, num1, num3, num2)
            coefs = permutedims(coefs, [1, 2, 4, 3])
        end

        #Insert knots along the u direction
        if isempty(iknots[1])
            knots[1] = nurbs.knots[1]
        else
            coefs = permutedims(coefs, [1, 3, 4, 2])
            coefs = reshape(coefs, 4*num2*num3, num1)
            coefs, knots[1] = bspkntins(degree[1], coefs, nurbs.knots[1], iknots[1])
            coefs = reshape(coefs, 4, num2, num3, size(coefs, 2))
            coefs = permutedims(coefs, [1, 4, 2, 3])
        end
    elseif length(nurbs.knots)==2
        # NURBS rrepresents a surface
        num1 = nurbs.number[1]
        num2 = nurbs.number[2]

        # Insert knots along the v direction
        if isempty(iknots[2])
            coefs = nurbs.coefs
            knots[2] = nurbs.knots[2]
        else
            coefs = reshape(nurbs.coefs, 4*num1, num2)
            coefs, knots[2] = bspkntins(degree[2], coefs, nurbs.knots[2], iknots[2])
            num2 = size(coefs, 2)
            coefs = reshape(coefs, 4, num1, num2)
        end

        # Insert knots along the u directions
        if isempty(iknots[1])
            knots[1] = nurbs.knots[1]
        else
            coefs = permutedims(coefs, [1, 3, 2])
            coefs = reshape(coefs, 4*num2, num1)
            coefs, knots[1] = bspkntins(degree[1], coefs, nurbs.knots[1], iknots[1])
            coefs = reshape(coefs, 4, num2, size(coefs,2))
            coefs = permutedims(coefs, [1, 3, 2])
        end
    else
        # NURBS represents a curve
        if isempty(iknots[1])
            coefs = nurbs.coefs
            knots = nurbs.knots
        else
            coefs, knots[1] = bspkntins(degree[1], nurbs.coefs, nurbs.knots[1], iknots[1])
        end
    end

    #construct the new NURBS
    inurbs = nrbmak(coefs, knots)
    return inurbs
end

"""
Joins two 1D NURBS objects (NURBS curves), nurbsA and nurbsB, with matching end control points. The
output is a single NURBS curve with scaled and normalized knot vector
"""
function nrbjoin(nurbsA::NURBS, nurbsB::NURBS)
    coefsA = nurbsA.coefs
    coefsB = nurbsB.coefs
    knotsA = nurbsA.knots[1]
    knotsB = nurbsB.knots[1]
    tol_eq = 1e-10
    if norm(coefsA[:,end]-coefsB[:,1])>tol_eq
        error("The two NURBS objects do not appear to have a common endpoint")
    end
    coefs_new = hcat(coefsA, coefsB[:,2:end])
    knotsB .+= nurbsA.knots[1][end]
    knots_new = vcat(knotsA, knotsB[2:end])
    #linearly map knots_new to the interval [0,1]
    knots_new = (knots_new .- knots_new[1])./(knots_new[end]-knots_new[1])
    nrb_new = nrbmak(coefs_new, [knots_new])
    return nrb_new
end

function bspdegelev(d, c, k, t)
    mc, nc = size(c)
    ic = zeros(mc, nc*(t+1))
    n = nc - 1
    bezalfs = zeros(d+1, d+t+1)
    bpts = zeros(mc, d+1)
    ebpts = zeros(mc, d+t+1)
    Nextbpts = zeros(mc, d+1)
    alfs = zeros(d)

    m = n + d + 1
    ph = d + t
    ph2 = floor(Int, ph/2)
    # Compute Bezier degree elevation coefficients
    bezalfs[1,1] = 1
    bezalfs[d+1, ph+1] = 1

    for i=1:ph2
        inv = 1/binomial(ph, i)
        mpi = min(d,i)

        for j=max(0, i-t):mpi
            bezalfs[j+1, i+1] = inv*binomial(d,j)*binomial(t, i-j)
        end
    end

    for i=ph2+1:ph-1
        mpi = min(d,i)
        for j=max(0,i-t):mpi
            bezalfs[j+1, i+1] = bezalfs[d-j+1, ph-i+1]
        end
    end

    mh = ph
    kind = ph+1
    r = -1
    a = d
    b = d + 1
    cind = 1
    ua = k[1]

    ic[1:mc,1] = c[1:mc,1]
    ik = ua.*ones(ph+1)

    # Initialize the first Bezier segment
    #bpts[1:mc, 1:d+1] = c[1:mc, 1:d+1]
    bpts = copy(c)

    # Big loop through knot vector
    while b<m
        i=b
        while b<m && k[b+1] == k[b+2]
            b += 1
        end
        mul = b - i + 1
        mh = mh + mul + t
        ub = k[b+1]
        oldr = r
        r = d - mul

        # Insert knot u[b] r times
        if oldr > 0
            lbz = floor(Int, (oldr+2)/2)
        else
            lbz = 1
        end

        if r>0
            rbz = ph - floor(Int, (r+1)/2)
        else
            rbz = ph
        end

        if r>0
            # Insert knot to get Bezier segment
            numer = ub - ua
            for q = d:-1:mul+1
                alfs[q-mul] = numer / (k[a+q+1]-ua)
            end

            for j = 1:r
                save = r-j
                s = mul + j

                for q = d:-1:s
                    for ii = 0:mc-1
                        tmp1 = alfs[q-s+1]*bpts[ii+1, q+1]
                        tmp2 = (1-alfs[q-s+1])*bpts[ii+1, q]
                        bpts[ii+1, q+1] = tmp1 + tmp2
                    end
                end

                Nextbpts[:, save+1] = bpts[:, d+1]
            end
        end
        # End of insert knot

        #Degree elevate Bezier
        for i=lbz:ph
            ebpts[:,i+1] = zeros(mc)
            mpi = min(d,i)
            for j=max(0,i-t):mpi
                for ii=0:mc-1
                    tmp1 = ebpts[ii+1, i+1]
                    tmp2 = bezalfs[j+1, i+1]*bpts[ii+1, j+1]
                    ebpts[ii+1, i+1] = tmp1 + tmp2
                end
            end
        end
        # End of degree elevating Bezier

        if oldr > 1
            # Must remove knot u=k[a] oldr times
            first = kind - 2
            last = kind
            den = ub - ua
            bet = floor(Int, (ub-ik[kind])/den)

            # Knot removal loop
            for tr = 1:oldr-1
                i = first
                j = last
                kj = j - kind + 1
                while j-i > tr
                    # Loop and compute the new control points for one removal step
                    if i < cind
                        alf = (ub - ik[i+1])/(ua - ik[i+1])
                        tmp1 = alf.*ic[:,i+1]
                        tmp2 = (1-alf).*ic[:, i]
                        ic[:,i+1] = tmp1 + tmp2
                    end
                    if j >= lbz
                        if j-tr <= kind - ph + oldr
                            gam = (ub-ik[j-tr+1])/den
                            tmp1 = gam.*ebpts[:, kj+1]
                            tmp2 = (1-gam).*ebpts[:, kj+2]
                            ebpts[:, kj+1] = tmp1 + tmp2
                        else
                            tmp1 = bet.*ebpts[:, kj+1]
                            tmp2 = (1-bet).*ebpts[:, kj+2]
                            ebpts[:, kj+1] = tmp1 + tmp2
                        end
                    end
                    i += 1
                    j -= 1
                    kj -= 1
                end
                first -= 1
                last += 1
            end
        end
        # End of removing knot n=k[a]

        # Load the knot ua
        if a!=d
            for i=0:ph-oldr-1
                push!(ik, ua)
                #ik[kind+1] = ua
                kind += 1
            end
        end

        for j = lbz:rbz
            for ii = 0:mc-1
                ic[ii+1, cind+1] = ebpts[ii+1, j+1]
            end
            cind += 1
        end

        if b < m
            # Setup for next pass through loop
            bpts[:,1:r] = Nextbpts[:, 1:r]
            bpts[:,r+1:d+1] = c[:, b-d+r+1:b+1]
            a = b
            b += 1
            ua = ub
        else
            for i=0:ph
                #ik[kind+i+1] = ub
                push!(ik, ub)
            end
        end
    end
    #End big while loop
    ic = ic[:, 1:cind]
    return ic, ik
end

function nrbdegelev(nurbs::NURBS, ntimes::Array{Int,1})
    degree = nurbs.order .- 1
    knots = Array{Array}(undef, length(nurbs.knots))
    if length(nurbs.knots)==3
        # NURBS represents a volume
        _, num1, num2, num3 = size(nurbs.coefs)

        # Degree elevate along the w direction
        if ntimes[3] == 0
            coefs = nurbs.coefs
            knots[3] = nurbs.knots[3]
        else
            coefs = reshape(nurbs.coefs, 4*num1*num2,num3)
            coefs, knots[3] = bspdegelev(degree[3], coefs, nurbs.knots[3], ntimes[3])
            num3 = size(coefs, 2)
            coefs = reshape(coefs, 4, num1, num2, num3)
        end

        # Degree elevate along the v direction
        if ntimes[2] == 0
            knots[2] = nurbs.knots[2]
        else
            coefs = permutedims(coefs, [1, 2, 4, 3])
            coefs = reshape(coefs, 4*num1*num3, num2)
            coefs, knots[2] = bspdegelev(degree[2], coefs, nurbs.knots[1], ntimes[1])
            num2 = size(coefs, 2)
            coefs = reshape(coefs, 4, num1, num3, num2)
            coefs = permutedims(coefs, [1, 2, 4, 3])
        end

        # Degree elevate along the u direction
        if ntimes[1] == 0
            knots[1] = nurbs.knots[1]
        else
            coefs = permutedims(coefs, [1, 3, 4, 2])
            coefs = reshape(coefs, 4*num2*num3, num1)
            coefs, knots[1] = bspdegelev(degree[1], coefs, nurbs.knots[1], ntimes[1])
            coefs = reshape(coefs, 4, num2, num3, size(coefs, 2))
            coefs = permutedims(coefs, [1, 4, 2, 3])
        end
    elseif length(nurbs.knots)==2
        _, num1, num2 = size(nurbs.coefs)

        # Degree elevate along the v direction
        if ntimes[2] == 0
            coefs = nurbs.coefs
            knots[2] = nurbs.knots[2]
        else
            coefs = reshape(nurbs.coefs, 4*num1, num2)
            coefs, knots[2] = bspdegelev(degree[2], coefs, nurbs.knots[2], ntimes[2])
            num2 = size(coefs, 2)
            coefs = reshape(coefs, 4, num1, num2)
        end

        # Degree elevate along the u direction
        if ntimes[1] == 0
            knots[1] = nurbs.knots[1]
        else
            coefs = permutedims(coefs, [1, 3, 2])
            coefs = reshape(coefs, 4*num2, num1)
            coefs, knots[1] = bspdegelev(degree[1], coefs, nurbs.knots[1], ntimes[1])
            coefs = reshape(coefs, 4, num2, size(coefs, 2))
            coefs = permutedims(coefs, [1, 3, 2])
        end
    else
        # NURBS represents a curve
        if isempty(ntimes)
            coefs = nurbs.coefs
            knots = nurbs.knots
        else
            coefs, knots[1] = bspdegelev(degree[1], nurbs.coefs, nurbs.knots[1], ntimes[1])
        end
    end

    #construct new NURBS
    inurbs = nrbmak(coefs, knots)
    return inurbs
end

function nrbruled(crv1::NURBS, crv2::NURBS)
    # Ensure both curves have a common degree
    d = max(crv1.order, crv2.order)[1]
    crv1 = nrbdegelev(crv1, [d - crv1.order[1]])
    crv2 = nrbdegelev(crv2, [d - crv2.order[1]])

    # Merge the knot vectors to obtain a common knot vector
    k1 = crv1.knots[1]
    k2 = crv2.knots[1]

    ku = unique(vcat(k1, k2))
    n = length(ku)

    ka = Array{Float64,1}(undef,0)
    kb = Array{Float64,1}(undef,0)

    for i = 1:n
        i1 = length(findall(k1.==ku[i]))
        i2 = length(findall(k2.==ku[i]))
        m = max(i1, i2)
        ka = vcat(ka, ku[i].*ones(m-i1))
        kb = vcat(kb, ku[i].*ones(m-i2))
    end
    crv1 = nrbkntins(crv1, [ka])
    crv2 = nrbkntins(crv2, [kb])

    coefs = cat(crv1.coefs, crv2.coefs, dims=3)
    srf = nrbmak(coefs, [crv1.knots[1], [0,0,1,1]])
    return srf
end

function nrbtform(nurbs::NURBS, tmat)
    if length(nurbs.knots)==2
        # NURBS is a surface
        dim, nu, nv = size(nurbs.coefs)
        coefs = reshape(tmat*reshape(nurbs.coefs, dim, nu*nv), dim, nu, nv)
    elseif length(nurbs.knots)==3
        # NURBS is a volume
        dim, nu, nv, nw = size(nurbs.coefs)
        coefs = reshape(tmat*reshape(nurbs.coefs, dim, nu*nv*nw), dim, nu, nv, nw)
    else
        # NURBS is a curve
        coefs = tmat*nurbs.coefs
    end
    tnurbs = nrbmak(coefs, nurbs.knots)
    return tnurbs
end

function nrbtransp(srf)
    if length(srf.knots)==1
        error("A NURBS curve cannot be transposed")
    elseif length(srf.knots)==3
        error("The transposition of NURBS volumes has not been implemented")
    end
    tsrf = nrbmak(permutedims(srf.coefs, [1,3,2]), reverse(srf.knots))
    return tsrf
end

function nrb4surf(p11, p12, p21, p22)
    coefs = cat(zeros(3,2,2), ones(1,2,2), dims=1)
    coefs[1:length(p11),1,1] = p11[:]
    coefs[1:length(p12),2,1] = p12[:]
    coefs[1:length(p21),1,2] = p21[:]
    coefs[1:length(p22),2,2] = p22[:]

    knots = [[0,0,1,1], [0,0,1,1]]
    srf = nrbmak(coefs, knots)
    return srf
end

function nrbcoons(u1::NURBS, u2::NURBS, v1::NURBS, v2::NURBS)
    tolEq = 1e-10
    if maximum(abs.(nrbeval(u1,u1.knots[1][[1]]) - nrbeval(v1,v1.knots[1][[1]])))>tolEq ||
            maximum(abs.(nrbeval(u1,u1.knots[1][[end]]) - nrbeval(v2,v2.knots[1][[1]])))>tolEq ||
            maximum(abs.(nrbeval(u2,u2.knots[1][[1]]) - nrbeval(v1,v1.knots[1][[end]])))>tolEq ||
            maximum(abs.(nrbeval(u2,u2.knots[1][[end]]) - nrbeval(v2,v2.knots[1][[end]])))>tolEq
        error("The four curves do not define a closed boundary")
    end

    r1 = nrbruled(u1, u2)
    r2 = nrbtransp(nrbruled(v1, v2))
    t = nrb4surf(u1.coefs[:,1], u1.coefs[:,end], u2.coefs[:,1], u2.coefs[:,end])

    # Raise all surfaces to a common degree
    du = max(r1.order[1], r2.order[1], t.order[1])
    dv = max(r1.order[2], r2.order[2], t.order[2])
    r1 = nrbdegelev(r1, [du - r1.order[1], dv - r1.order[2]])
    r2 = nrbdegelev(r2, [du - r2.order[1], dv - r2.order[2]])
    t = nrbdegelev(t, [du - t.order[1], dv - t.order[2]])

    # Merge the knots vectors, to obtain a commmon knot vector

    #Uknots
    k1 = r1.knots[1]
    k2 = r2.knots[1]
    k3 = t.knots[1]
    k = unique(vcat(k1, k2, k3))
    n = length(k)
    kua = Array{Float64,1}(undef,0)
    kub = Array{Float64,1}(undef,0)
    kuc = Array{Float64,1}(undef,0)

    for i=1:n
        i1 = length(findall(k1.==k[i]))
        i2 = length(findall(k2.==k[i]))
        i3 = length(findall(k3.==k[i]))
        m = max(i1, i2, i3)
        kua = vcat(kua, k[i].*ones(m-i1))
        kub = vcat(kub, k[i].*ones(m-i2))
        kuc = vcat(kuc, k[i].*ones(m-i3))
    end

    # Vknots
    k1 = r1.knots[2]
    k2 = r2.knots[2]
    k3 = t.knots[2]

    k = unique(vcat(k1, k2, k3))
    n = length(k)

    kva = Array{Float64,1}(undef,0)
    kvb = Array{Float64,1}(undef,0)
    kvc = Array{Float64,1}(undef,0)

    for i=1:n
        i1 = length(findall(k1.==k[i]))
        i2 = length(findall(k2.==k[i]))
        i3 = length(findall(k3.==k[i]))
        m = max(i1, i2, i3)
        kva = vcat(kva, k[i].*ones(m-i1))
        kvb = vcat(kvb, k[i].*ones(m-i2))
        kvc = vcat(kvc, k[i].*ones(m-i3))
    end

    r1 = nrbkntins(r1, [kua, kva])
    r2 = nrbkntins(r2, [kub, kvb])
    t = nrbkntins(t, [kuc, kvc])

    # combine coefficients to construct the Coons surface
    coefs = zero(r1.coefs)
    coefs[1,:,:] = r1.coefs[1,:,:] + r2.coefs[1,:,:] - t.coefs[1,:,:]
    coefs[2,:,:] = r1.coefs[2,:,:] + r2.coefs[2,:,:] - t.coefs[2,:,:]
    coefs[3,:,:] = r1.coefs[3,:,:] + r2.coefs[3,:,:] - t.coefs[3,:,:]
    coefs[4,:,:] = r1.coefs[4,:,:] + r2.coefs[4,:,:] - t.coefs[4,:,:]
    srf = nrbmak(coefs, r1.knots)
end

function nrbrevolve(curve::NURBS, pnt, vec, theta)
    if length(curve.knots)==3
        error("The function nrbrevolve is not yet ready to create volumes")
    end

    # Translate the center point to the origin
    if isempty(pnt)
        pnt = zeros(3)
    end

    if length(pnt)!=3
        error("All point and vecto coordinates must be 3D")
    end

    # Translate and rotate the original curve or surface into alignment with the z-axis
    T = vectrans(-pnt)
    angx = vecangle(vec[1], vec[3])
    RY = vecroty(-angx[1])
    vectmp = RY*vcat(vecnorm(vec), 1.)
    angy = vecangle(vectmp[2], vectmp[3])
    RX = vecrotx(angy[1])
    curve = nrbtform(curve, RX*RY*T)

    # Construct an arc
    arc = nrbcirc(1.0, [], 0.0, theta)
    if length(curve.knots)==2
        # Construct the revolved volume
        coefs = zeros(4, arc.number[1], curve.number[1], curve.number[2])
        angle = vecangle(curve.coefs[2,:,:], curve.coefs[1,:,:])
        radius = dropdims(vecmag(curve.coefs[1:2, :, :]), dims=1)
        for i=1:curve.number[1]
            for j=1:curve.number[2]
                coefs[:,:,i,j] = vecrotz(angle[i,j])*vectrans([0., 0., curve.coefs[3,i,j]])*
                            vecscale([radius[i,j], radius[i,j]])*arc.coefs
                coefs[4,:,i,j] = coefs[4, :, i, j]*curve.coefs[4,i,j]
            end
        end
        surf = nrbmak(coefs, vcat(arc.knots, curve.knots))
    else
        # Construct the revolved surface
        coefs = zeros(4, arc.number[1], curve.number[1])
        angle = vecangle(curve.coefs[2,:], curve.coefs[1,:])
        radius = vecmag(curve.coefs[1:2, :])
        for i = 1:curve.number[1]
            coefs[:, :, i] = vecrotz(angle[i])*vectrans([0., 0., curve.coefs[3,i]])*
                        vecscale([radius[i], radius[i]])*arc.coefs
            coefs[4, :, i] = coefs[4, :, i]*curve.coefs[4, i]
        end
        surf = nrbmak(coefs, vcat(arc.knots, curve.knots))
    end

    # Rotate and vectrans the surface back into position
    T = vectrans(pnt)
    RX = vecrotx(-angy[1])
    RY = vecroty(angx[1])
    surf = nrbtform(surf, T*RY*RX)
    return surf
end

function nrbrevolve(curve::NURBS, pnt, vec)
    surf = nrbrevolve(curve, pnt, vec, 2*pi)
    return surf
end

"""
Refine a given knot vector by dividing each interval uniformly,
             maintaining the continuity in previously existing knots.
"""
function kntrefine(knots, n_sub, degree, regularity)
    if length(n_sub)!=length(degree) || length(n_sub)!=length(regularity) ||
        length(n_sub)!=length(knots)
        error("kntrefine: n_sub, degree and regularity must have the same length
            as the number of knot vectors")
    end
    rknots = Array{Array}(undef, length(knots))
    new_knots = Array{Array}(undef, length(knots))
    zeta = Array{Array}(undef, length(knots))
    for idim = 1:length(n_sub)
        min_mult = degree[idim] - regularity[idim]
        z = unique(knots[idim])
        nz = length(z)
        deg = sum(knots[idim].==z[1])-1
        rknots[idim] = z[ones(Int, degree[idim]+1)]
        new_knots[idim] = [];
        for ik = 2:nz
            insk = LinRange(z[ik-1], z[ik], n_sub[idim]+2)
            insk = vec(repeat(insk[2:end-1], min_mult,1))
            old_mult = sum(knots[idim].==z[ik])
            mult = max(min_mult, degree[idim] - deg + old_mult)

            rknots[idim] = vcat(rknots[idim], insk, z[ik*ones(Int, mult)])
            if mult>=old_mult
                new_knots[idim] = vcat(new_knots[idim], insk, z[ik*ones(Int, mult-old_mult)])
            else
                new_knots[idim] = vcat(new_knots[idim], insk)
            end
        end
        zeta[idim] = unique(rknots[idim])
    end
    return rknots, zeta, new_knots
end

function nrbsquare(corner, lengthx, lengthy, degree, nsubdiv)
    if isempty(corner)
        corner = [0,0]
    end

    if length(degree)==1
        deg = [degree, degree]
    elseif length(degree)==2
        deg = degree
    else
        error("The degree vector should provide the degree in each direction (one or two values)")
    end

    if length(nsubdiv)==1
        nsub = [nsubdiv, nsubdiv]
    elseif length(nsubdiv)==2
        nsub = nsubdiv
    else
        error("The nsubdiv vector should provide the number of intervals in each direction (one or two values)")
    end

    srf = nrb4surf(corner, corner+[lengthx, 0], corner+[0, lengthy], corner+[lengthx, lengthy])
    srf = nrbdegelev(srf, deg-[1, 1])
    _, _, new_knots = kntrefine(srf.knots, nsub.-1, deg, deg - [1, 1] )
    srf = nrbkntins(srf, new_knots)
    return srf
end

"""
Reverse the evaluation direction of a NURBS geometry
Only implemented for curves for now
"""
function nrbreverse(nrb)
    if length(nrb.knots)==1
        #reverse a NURBS curve
        new_knots = sort(nrb.knots[1][end].-nrb.knots[1])
        new_coefs = reverse(nrb.coefs, dims=2)
        new_nrb = nrbmak(new_coefs, [new_knots])
    end
    return new_nrb
end


text:
using SoftSquishyMatter # load simulation package
using Random # useful if you want to shuffle arrays
cd(dirname(@__FILE__)) # set directory to this file's folder

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Initialize empty simulation with a description for future reference
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
simulation = Simulation()
simulation.descriptor = "Passive colloids with Lennard-Jones interactions"

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Initialize particles
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
R = 1e-6 # radius of particle
γ_trans = 6 * pi * 8.9e-4 * R # friction coefficient of particle
D_trans = 1.38e-23 * 300 / γ_trans # diffusivity of particle

# initialize positions of particles on a triangular lattice
# use lattice dimensions for the simulation dimensions (L_x, L_y)
positions, simulation.L_x, simulation.L_y = triangular_lattice(s = 3e-6, M_x = 30, M_y = 20)

# create particles with properties
for (x, y) in positions
    push!(simulation.particles, Particle(ptype = :passivecolloid, x = x, y = y, R = R, γ_trans = γ_trans, D_trans = D_trans))
end
# particle groups are useful for interactions and integrators
# in this case, it doesn't really matter
pgroup_all = simulation.particles

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Initialize cell list and Lennard-Jones interactions
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
cell_list = CellList(particles = pgroup_all, L_x = simulation.L_x, L_y = simulation.L_y, cutoff = 2^(1.0 / 6.0) * 2 * R)
lj = LennardJones(particles = pgroup_all, cell_list = cell_list, ϵ = D_trans * γ_trans, multithreaded = true)
push!(simulation.cell_lists, cell_list)
push!(simulation.interactions, lj)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Initialize Brownian integrator for dynamics
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
simulation.dt = 0.001
brownian = Brownian(particles = pgroup_all, dt = simulation.dt, multithreaded = true)
push!(simulation.integrators, brownian)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Set the number of steps and which particles to save periodically
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
simulation.num_steps = trunc(Int64, 10 / simulation.dt)
simulation.save_interval = trunc(Int64, 1 / simulation.dt)
simulation.particles_to_save = pgroup_all

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Run simulation and save data
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
run_simulation!(simulation; save_as = "out/Example_LennardJonesFluid_data.out")

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# OPTIONAL: Load simulation data and animate the simulation.
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
simulation = load_simulation(file = "out/Example_LennardJonesFluid_data.out")
visualize!(simulation; save_as = (:gif, "frames/Example_LennardJonesFluid.gif"), fps = 1, particle_colors = Dict(:passivecolloid => "black"))


text:
#
include("objects.jl")
import Base.show

mutable struct Collector
    health
    player  # who owns this unit
    pos::WorldPos
    Collector(player) = new(max_health(Collector), player, WorldPos(0,0))
    Collector(player, pos) = new(max_health(Collector), player, pos)
end
mutable struct Fighter
    health
    player  # who owns this unit
    pos::WorldPos
    attackTargetUnit  # who currently attacking
    Fighter(player) = new(max_health(Fighter), player, WorldPos(0,0), nothing)
    Fighter(player, pos) = new(max_health(Fighter), player, pos, nothing)
end

Base.show(io::IO, u::Collector) = print(io, "Collector($(u.health), $(Ptr{PlayerUnits}(pointer_from_objref(u.player))), $(u.pos)")
health_percent(u) = health(u) / max_health(typeof(u))

health(x::Collector) = x.health
set_health!(x::Collector, h) = (x.health = h)

health(x::Fighter) = x.health
set_health!(x::Fighter, h) = (x.health = h)

# Unit attributes
@noinline max_health(::Type{Collector}) = (global kMaxHealth_Collector; kMaxHealth_Collector)
@noinline max_health(::Type{Fighter}) = (global kMaxHealth_Collector; kMaxHealth_Fighter)

@noinline attack_damage(::Collector) = (global kMaxHealth_Collector; kAttackDamage_Collector)
@noinline attack_damage(::Fighter) = (global kMaxHealth_Collector; kAttackDamage_Fighter)

@noinline build_time(::Type{Collector}) = (global kMaxHealth_Collector; kBuildTime_Collector)
@noinline build_time(::Type{Fighter}) = (global kMaxHealth_Collector; kBuildTime_Fighter)
@noinline build_cost(::Type{Collector}) = (global kMaxHealth_Collector; kBuildCost_Collector)
@noinline build_cost(::Type{Fighter}) = (global kMaxHealth_Collector; kBuildCost_Fighter)

unit_types = (Collector, Fighter)
UnitTypes = Union{unit_types...}

# ---------------------
# PLAYER
# --------------------
mutable struct PlayerUnits
    units::Array
    collectors::Array{Collector}
    fighters::Array{Fighter}
end
PlayerUnits() = PlayerUnits([],[],[])
Base.show(io::IO, pu::PlayerUnits) = print(io,
                     "PlayerUnits(\n"
                    *"            units[$(length(pu.units))]\n"
                    *"            collectors[$(length(pu.collectors))]\n"
                    *"            fighters[$(length(pu.fighters))]\n"
                    *"           )")

function add_unit_helper(p::PlayerUnits, u)
    @assert u.player == p
    push!(p.units, u)
    return u
end
function add_unit!(p::PlayerUnits, u)
    return add_unit_helper(p,u)
end
function add_unit!(p::PlayerUnits, u::Collector)
    add_unit_helper(p,u)
    push!(p.collectors, u)
    return u
end
function add_unit!(p::PlayerUnits, u::Fighter)
    add_unit_helper(p,u)
    push!(p.fighters, u)
    return u
end
add_collector!(p::PlayerUnits) = add_unit!(p, Collector(p))
add_fighter!(p::PlayerUnits) = add_unit!(p, Fighter(p))

remove_unit!(p, unit) = filter!(u->u≠unit, p.units)
function remove_unit!(p, unit::Collector)
    filter!(u->u≠unit, p.collectors)
    filter!(u->u≠unit, p.units)
end
function remove_unit!(p, unit::Fighter)
    filter!(u->u≠unit, p.fighters)
    filter!(u->u≠unit, p.units)
end

function clear_units!(p)
    empty!(p.units)
    empty!(p.collectors)
    empty!(p.fighters)
end
function clear_collectors!(p)
    filter!(u->!isa(u, Collector), p.units)
    empty!(p.collectors)
end
function clear_fighters!(p)
    filter!(u->!isa(u, Fighter), p.units)
    empty!(p.fighters)
end

owning_player(x)::PlayerUnits = x.player

# --------------

function attack!(target, attacker)
    set_health!(target, health(target) - attack_damage(attacker))
    if (health(target) <= 0)
        destroy_unit!(target)
    end
    return health(target)
end

function destroy_unit!(target)
    remove_unit!(owning_player(target), target)
end

# Tests
using Test
collectors = PlayerUnits()
add_collector!(collectors)
destroy_unit!(collectors.units[1])
@test 0 == length(collectors.units)
w = add_collector!(collectors)
w = add_collector!(collectors)
@test 2 == length(collectors.units)
attack!(w, Fighter(nothing))
attack!(w, Fighter(nothing))
@test 1 == length(collectors.units)


text:
function fTBGrd(x,Param)
if lower(Param.ProfTheta)== "baldaufcart"
    T=Param.T0;
end
return T
end


text:
"""
$(TYPEDEF)

Physical space with unstructured mesh

# Fields

$(FIELDS)
"""
struct UnstructPSpace{
    A,
    B<:AM{<:FN},
    C<:AM{<:Integer},
    D<:AV{<:Integer},
    E<:AV{<:FN},
    F<:AA{<:FN,3},
} <: AbstractPhysicalSpace
    cells::A # all information: cell, line, vertex
    points::B # locations of vertex points
    cellid::C # node indices of elements
    cellType::D # inner/boundary cell
    cellNeighbors::C # neighboring cells id
    cellFaces::C # cell edges id
    cellCenter::B # cell center location
    cellArea::E # cell size
    cellNormals::F # cell unit normal vectors
    facePoints::C # ids of two points at edge
    faceCells::C # ids of two cells around edge
    faceCenter::B # edge center location
    faceType::D # inner/boundary face
    faceArea::E # face area
end

function UnstructPSpace(file::T) where {T<:AbstractString}
    cells, points = read_mesh(file)
    p = mesh_connectivity_2D(cells, points)

    return UnstructPSpace((cells, points)..., p...)
end

function write_vtk(ks::T1, ctr) where {T1<:AbstractSolverSet}
    write_vtk(ks.ps, ctr)

    return nothing
end

function write_vtk(ps::UnstructPSpace, ctr)
    cdata = zeros(length(ctr), length(ctr[1].w))
    for i in eachindex(ctr)
        cdata[i, :] .= ctr[i].prim
        if size(cdata, 2) > 1
            cdata[i, end] = 1.0 / cdata[i, end]
        end
    end
    write_vtk(ps.points, ps.cellid, cdata)

    return nothing
end


text:
### Instruction:
# This file runs the main estimations in serial-processing mode to isolate
# the source of problem in parallel computation setting
# For the list of input parameters, see gmmOptim! function in GMMInterface.jl
# For the list of input files, check all *.csv files in Helper.jl
# The post estimation step is executed by main.jl

push!(LOAD_PATH,pwd())

using TypeDef, Helper, Estim, Gmm, GMMInterface, UserFileIO, UnitTest
using PostEstim

const ver   = "824"
const fCSV1 = "gmmStage1.csv"
const fCSV2 = "gmmStage2.csv"
const fJLD1 = "gmmStage1.jld"
const fJLD2 = "gmmStage2.jld"
const fINT1 = "gmmStage1Interim.csv"
const fINT2 = "gmmStage2Interim.csv"
const fileName = fileNameType(fCSV1,fCSV2,fJLD1,fJLD2,fINT1,fINT2)

### Continuous updating optimal IV for GMM (Tables 4, A.1 & A.15)
gmmOptim!(fileName, 0,0,"opt",ver,"base";update=true)
# gmmOptim!(fileName,15,0,"opt",ver,"base";update=true)
# gmmOptim!(fileName,27,0,"opt",ver,"base";update=true,nPerson=3000)
# gmmOptim!(fileName,15,0,"opt",ver,"ms15";update=true)
# gmmOptim!(fileName,15,0,"opt",ver,"noAllow";update=true)
# gmmOptim!(fileName,15,0,"opt",ver,"extra";update=true)
#
# ### Comparison of alternative IVs (Table A.17)
# for model=[0,15]
# 	# GMM with BLP IV
# 	gmmOptim!(fileName,model,1,"blp",ver,"base";update=false)
# 	gmmOptim!(fileName,model,2,"blp",ver,"base";update=false,mcid=1)
#
# 	# GMM with diff-quad IV
# 	gmmOptim!(fileName,model,1,"diff-quad",ver,"base";update=false)
# 	gmmOptim!(fileName,model,2,"diff-quad",ver,"base";update=false,mcid=1)
#
# 	# GMM with diff-local IV
# 	gmmOptim!(fileName,model,1,"diff-local",ver,"base";update=false)
# 	gmmOptim!(fileName,model,2,"diff-local",ver,"base";update=false,mcid=1)
# end

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
export RecentSubject, RecentSubjectFactory

import Base: show, similar

"""
    RecentSubject(::Type{D}) where D
    RecentSubject(::Type{D}, factory::F) where { D, F <: AbstractSubjectFactory }
    RecentSubject(::Type{D}, subject::S) where { D, S }

A variant of Subject that emits its recent value whenever it is subscribed to.

See also: [`RecentSubjectFactory`](@ref), [`Subject`](@ref), [`SubjectFactory`](@ref)
"""
function RecentSubject end

"""
    RecentSubjectFactory(factory::F) where { F <: AbstractSubjectFactory }
    RecentSubjectFactory(; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler }

A variant of SubjectFactory that creates an instance of RecentSubject.

See also: [`SubjectFactory`](@ref), [`AbstractSubjectFactory`](@ref), [`RecentSubject`](@ref), [`Subject`](@ref)
"""
function RecentSubjectFactory end

##

mutable struct RecentSubjectInstance{D, S} <: AbstractSubject{D}
    subject :: S
    recent  :: Union{D, Nothing}
end

Base.show(io::IO, ::RecentSubjectInstance{D, S}) where { D, S } = print(io, "RecentSubject($D, $S)")

Base.similar(subject::RecentSubjectInstance{D, S}) where { D, S } = RecentSubject(D, similar(subject.subject))

function RecentSubject(::Type{D}) where D
    return RecentSubject(D, SubjectFactory(AsapScheduler()))
end

function RecentSubject(::Type{D}, factory::F) where { D, F <: AbstractSubjectFactory }
    return RecentSubject(D, create_subject(D, factory))
end

function RecentSubject(::Type{D}, subject::S) where { D, S }
    return as_recent_subject(D, as_subject(S), subject)
end

as_recent_subject(::Type{D},  ::InvalidSubjectTrait,    subject)    where D          = throw(InvalidSubjectTraitUsageError(subject))
as_recent_subject(::Type{D1}, ::ValidSubjectTrait{D2},  subject)    where { D1, D2 } = throw(InconsistentSubjectDataTypesError{D1, D2}(subject))
as_recent_subject(::Type{D},  ::ValidSubjectTrait{D},   subject::S) where { D, S }   = RecentSubjectInstance{D, S}(subject, nothing)

##

getrecent(subject::RecentSubjectInstance)         = subject.recent
setrecent!(subject::RecentSubjectInstance, value) = subject.recent = value

##

function on_next!(subject::RecentSubjectInstance{D}, data::D) where D
    setrecent!(subject, data)
    next!(subject.subject, data)
end

function on_error!(subject::RecentSubjectInstance, err)
    error!(subject.subject, err)
end

function on_complete!(subject::RecentSubjectInstance)
    complete!(subject.subject)
end

##

function on_subscribe!(subject::RecentSubjectInstance, actor)
    recent = getrecent(subject)
    if recent !== nothing
        next!(actor, recent)
    end
    return subscribe!(subject.subject, actor)
end

##

struct RecentSubjectFactoryInstance{ F <: AbstractSubjectFactory } <: AbstractSubjectFactory
    factory :: F
end

Base.show(io::IO, ::RecentSubjectFactoryInstance{F}) where F = print(io, "RecentSubjectFactory($F)")

create_subject(::Type{L}, factory::RecentSubjectFactoryInstance) where L = RecentSubject(L, factory.factory)

function RecentSubjectFactory(factory::F) where { F <: AbstractSubjectFactory }
    return RecentSubjectFactoryInstance(factory)
end

function RecentSubjectFactory(; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler }
    return RecentSubjectFactoryInstance(SubjectFactory{H}(scheduler))
end

##