text:
using Latexify
using Test
using Markdown

str = "2*x^2 - y/c_2"
ex = :(2*x^2 - y/c_2)

desired_output = "2 \\cdot x^{2} - \\frac{y}{c_{2}}"

@test latexify(:(a = [x / y; 3; 4])) == replace(
raw"$a = \left[
\begin{array}{c}
\frac{x}{y} \\
3 \\
4 \\
\end{array}
\right]$", "\r\n"=>"\n")

@test latexify(:(a = [x / y 2 3 4])) == replace(
raw"$a = \left[
\begin{array}{cccc}
\frac{x}{y} & 2 & 3 & 4 \\
\end{array}
\right]$", "\r\n"=>"\n")

@test latexify(:(a = [x / y 2; 3 4])) == replace(
raw"$a = \left[
\begin{array}{cc}
\frac{x}{y} & 2 \\
3 & 4 \\
\end{array}
\right]$", "\r\n"=>"\n")

@test latexraw(str) == latexraw(ex)
@test latexraw(ex) == desired_output

@test latexify(:(u[1, 2]); index = :bracket) == raw"$u\left[1, 2\right]$"
@test latexify(:(u[1, 2]); index = :subscript) == raw"$u_{1,2}$"

array_test = [ex, str]
@test all(latexraw.(array_test) .== desired_output)

@test latexraw(:(@__dot__(x / y))) == raw"\frac{x}{y}"
@test latexraw(:(@. x / y)) == raw"\frac{x}{y}"
@test latexraw(:(eps())) == raw"\mathrm{eps}()"

@test latexraw(:y_c_a) == "y_{c\\_a}"
@test latexraw(1.0) == "1.0"
@test latexraw(1.00) == "1.0"
@test latexraw(1) == "1"
@test latexraw(:(log10(x))) == "\\log_{10}\\left( x \\right)"
@test latexraw(:(sin(x))) ==  "\\sin\\left( x \\right)"
@test latexraw(:(sin(x)^2)) ==  "\\sin^{2}\\left( x \\right)"
@test latexraw(:(asin(x))) ==  "\\arcsin\\left( x \\right)"
@test latexraw(:(asinh(x))) ==  "\\mathrm{arcsinh}\\left( x \\right)"
@test latexraw(:(sinc(x))) ==  "\\mathrm{sinc}\\left( x \\right)"
@test latexraw(:(acos(x))) ==  "\\arccos\\left( x \\right)"
@test latexraw(:(acosh(x))) ==  "\\mathrm{arccosh}\\left( x \\right)"
@test latexraw(:(cosc(x))) ==  "\\mathrm{cosc}\\left( x \\right)"
@test latexraw(:(atan(x))) ==  "\\arctan\\left( x \\right)"
@test latexraw(:(atan2(x))) ==  "\\arctan\\left( x \\right)"
@test latexraw(:(atanh(x))) ==  "\\mathrm{arctanh}\\left( x \\right)"
@test latexraw(:(acot(x))) ==  "\\mathrm{arccot}\\left( x \\right)"
@test latexraw(:(acoth(x))) ==  "\\mathrm{arccoth}\\left( x \\right)"
@test latexraw(:(asec(x))) ==  "\\mathrm{arcsec}\\left( x \\right)"
@test latexraw(:(sech(x))) ==  "\\mathrm{sech}\\left( x \\right)"
@test latexraw(:(asech(x))) ==  "\\mathrm{arcsech}\\left( x \\right)"
@test latexraw(:(acsc(x))) ==  "\\mathrm{arccsc}\\left( x \\right)"
@test latexraw(:(csch(x))) ==  "\\mathrm{csch}\\left( x \\right)"
@test latexraw(:(acsch(x))) ==  "\\mathrm{arccsch}\\left( x \\right)"
@test latexraw(:(x ± y)) == "x \\pm y"
@test latexraw(:(f(x))) ==  "f\\left( x \\right)"
@test latexraw("x = 4*y") == "x = 4 \\cdot y"
@test latexraw(:(sqrt(x))) == "\\sqrt{x}"
@test latexraw(complex(1,-1)) == "1-1\\textit{i}"
@test latexraw(1//2) == "\\frac{1}{2}"
@test latexraw(missing) == "\\textrm{NA}"
@test latexraw("x[2]") == raw"x\left[2\right]"
@test latexraw("x[2, 3]") == raw"x\left[2, 3\right]"
@test latexraw("α") == raw"\alpha"
@test latexraw("α + 1") == raw"\alpha + 1"
@test latexraw("α₁") == raw"\alpha{_1}"
@test latexraw("γ³") == raw"\gamma{^3}"
@test latexraw("β₃_hello") == raw"\beta{_3}_{hello}"
@test latexraw("β₃₄") == raw"\beta{_{34}}"
@test latexraw("β₃₄¹⁴") == raw"\beta{_{34}^{14}}"
@test latexraw("β₃¹⁴") == raw"\beta{_3^{14}}"
@test latexraw("β¹⁴₃") == raw"\beta{^{14}_3}"
@test latexraw("β¹⁴") == raw"\beta{^{14}}"
@test latexraw("β⁴") == raw"\beta{^4}"
@test latexraw(Inf) == raw"\infty"
@test latexraw(:Inf) == raw"\infty"
@test latexraw("Inf") == raw"\infty"
@test latexraw(-Inf) == raw"-\infty"
@test latexraw(:($(3+4im)*a)) == raw"\left( 3+4\textit{i} \right) \cdot a"
@test latexraw(:(a*$(3+4im))) == raw"a \cdot \left( 3+4\textit{i} \right)"

@test latexify(:((-1) ^ 2)) == replace(
raw"$\left( -1 \right)^{2}$", "\r\n"=>"\n")
@test latexify(:($(1 + 2im) ^ 2)) == replace(
raw"$\left( 1+2\textit{i} \right)^{2}$", "\r\n"=>"\n")
@test latexify(:($(3 // 2) ^ 2)) == replace(
raw"$\left( \frac{3}{2} \right)^{2}$", "\r\n"=>"\n")

### Test broadcasting
@test latexraw(:(sum.((a, b)))) == raw"\mathrm{sum}\left( a, b \right)"

### Test for correct signs in nested sums/differences.
@test latexraw("-(-1)") == raw" + 1"
@test latexraw("+(-1)") == raw"-1"
@test latexraw("-(+1)") == raw" - 1"
@test latexraw("-(1+1)") == raw" - \left( 1 + 1 \right)"
@test latexraw("1-(-2)") == raw"1 + 2"
@test latexraw("1 + (-(2))") == raw"1 - 2"
@test latexraw("1 + (-2 -3 -4)") == raw"1 -2 - 3 - 4"
@test latexraw("1 - 2 - (- 3 - 4)") == raw"1 - 2 - \left(  - 3 - 4 \right)"
@test latexraw("1 - 2 - (- 3 -(2) + 4)") == raw"1 - 2 - \left(  - 3 - 2 + 4 \right)"
@test latexraw("1 - 2 - (- 3 -(2 - 8) + 4)") == raw"1 - 2 - \left(  - 3 - \left( 2 - 8 \right) + 4 \right)"
@test latexraw(:(-$(3+5im))) == raw" - \left( 3+5\textit{i} \right)"
@test latexraw(:($(3+4im)-$(3+5im))) == raw"3+4\textit{i} - \left( 3+5\textit{i} \right)"

# @test_throws ErrorException latexify("x/y"; env=:raw, bad_kwarg="should error")

@test latexraw(:(3 * (a .< b .<= c < d <= e > f <= g .<= h .< i == j .== k != l .!= m))) ==
raw"3 \cdot \left( a < b \leq c < d \leq e > f \leq g \leq h < i = j = k \neq l \neq m \right)"

#### Test the fmt keyword option
@test latexify([32894823 1.232212 :P_1; :(x / y) 1.0e10 1289.1]; env=:align, fmt="%.2e") == replace(
raw"\begin{align}
3.29e+07 =& 1.23e+00 =& P_{1} \\
\frac{x}{y} =& 1.00e+10 =& 1.29e+03
\end{align}
", "\r\n"=>"\n")

@test latexify([32894823 1.232212 :P_1; :(x / y) 1.0e10 1289.1]; env=:equation, fmt="%.2e") == replace(
raw"\begin{equation}
\left[
\begin{array}{ccc}
3.29e+07 & 1.23e+00 & P_{1} \\
\frac{x}{y} & 1.00e+10 & 1.29e+03 \\
\end{array}
\right]
\end{equation}
", "\r\n"=>"\n")

@test latexify([32894823 1.232212 :P_1; :(x / y) 1.0e10 1289.1]; env=:table, fmt="%.2e") == replace(
raw"\begin{tabular}{ccc}
$3.29e+07$ & $1.23e+00$ & $P_{1}$\\
$\frac{x}{y}$ & $1.00e+10$ & $1.29e+03$\\
\end{tabular}
", "\r\n"=>"\n")

@test latexify([32894823 1.232212 :P_1; :(x / y) 1.0e10 1289.1]; env=:mdtable, fmt="%.2e") ==
Markdown.md"|    $3.29e+07$ | $1.23e+00$ |    $P_{1}$ |
| -------------:| ----------:| ----------:|
| $\frac{x}{y}$ | $1.00e+10$ | $1.29e+03$ |
"

@test latexify(1234.2234; env=:inline, fmt="%.2e") ==
raw"$1.23e+03$"

@test latexify(1234.2234; env=:raw, fmt="%.2e") ==
raw"1.23e+03"

test_functions = [:sinh, :alpha, :Theta, :cosc, :acoth, :acot, :asech, :lambda,
                  :asinh, :sinc, :eta, :kappa, :nu, :asin, :epsilon, :sigma,
                  :upsilon, :phi, :tanh, :iota, :Psi, :acosh, :log, :zeta, :mu,
                  :csc, :xi, :tau, :beta, :Lambda, :Xi, :Phi, :acsc, :atan,
                  :sech, :atanh, :gamma, :Delta, :rho, :sec, :log10, :delta,
                  :pi, :cot, :log2, :cos, :Omega, :psi, :atan2, :Gamma, :cosh,
                  :acos, :Pi, :Upsilon, :omega, :coth, :chi, :tan, :csch,
                  :acsch, :theta, :asec, :Sigma, :sin]

@test latexify(["3*$(func)(x)^2/4 -1" for func = test_functions]) == replace(
raw"\begin{equation}
\left[
\begin{array}{c}
\frac{3 \cdot \sinh^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \alpha\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \Theta\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \mathrm{cosc}^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \mathrm{arccoth}^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \mathrm{arccot}^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \mathrm{arcsech}^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \lambda\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \mathrm{arcsinh}^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \mathrm{sinc}^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \eta\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \kappa\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \nu\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \arcsin^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \epsilon\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \sigma\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \upsilon\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \phi\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \tanh^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \iota\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \Psi\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \mathrm{arccosh}^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \log\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \zeta\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \mu\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \csc^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \xi\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \tau\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \beta\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \Lambda\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \Xi\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \Phi\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \mathrm{arccsc}^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \arctan^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \mathrm{sech}^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \mathrm{arctanh}^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \Gamma\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \Delta\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \rho\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \sec^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \log_{10}\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \delta\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \pi\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \cot^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \log_{2}\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \cos^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \Omega\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \psi\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \arctan^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \Gamma\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \cosh^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \arccos^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \Pi\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \Upsilon\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \left( \omega\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \coth^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \chi\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \tan^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \mathrm{csch}^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \mathrm{arccsch}^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \theta\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \mathrm{arcsec}^{2}\left( x \right)}{4} - 1 \\
\frac{3 \cdot \left( \Sigma\left( x \right) \right)^{2}}{4} - 1 \\
\frac{3 \cdot \sin^{2}\left( x \right)}{4} - 1 \\
\end{array}
\right]
\end{equation}
", "\r\n"=>"\n")

## Test logical operators
@test latexraw(:(x && y)) == "x \\wedge y"
@test latexraw(:(x || y)) == "x \\vee y"
@test latexraw(:(x || !y)) == "x \\vee \\neg y"

## Test {cases} enviroment
@test latexraw(:(R(p,e,d) = e ? 0 : log(p) - d)) == replace(
raw"R\left( p, e, d \right) = \begin{cases}
0 & \text{if } e\\
\log\left( p \right) - d & \text{otherwise}
\end{cases}", "\r\n"=>"\n")

@test latexraw(:(R(p,e,d,t) = if (t && e); 0 elseif (t && !e); d else log(p) end)) == replace(
raw"R\left( p, e, d, t \right) = \begin{cases}
0 & \text{if } t \wedge e\\
d & \text{if } t \wedge \neg e\\
\log\left( p \right) & \text{otherwise}
\end{cases}", "\r\n"=>"\n")

@test latexraw(:(function reward(p,e,d,t)
    if t && e
        return 0
    elseif t && !e
        return -1d
    elseif 2t && e
        return -2d
    elseif 3t && e
        return -3d
    else
        return log(p)
    end
end)) == replace(
raw"\mathrm{reward}\left( p, e, d, t \right) = \begin{cases}
0 & \text{if } t \wedge e\\
-1 \cdot d & \text{if } t \wedge \neg e\\
-2 \cdot d & \text{if } 2 \cdot t \wedge e\\
-3 \cdot d & \text{if } 3 \cdot t \wedge e\\
\log\left( p \right) & \text{otherwise}
\end{cases}", "\r\n"=>"\n")


text:
module Controller

using Modia
using Modia.StaticArrays

# Single-Input-Single-Output continuous control block
SISO = Model(
    u = input,
    y = output
)

# Gain
Gain = SISO | Model(
    k = 1.0,
    equations = :[
        y = k*u ]
)

# First-order transfer function block (= 1 pole)
FirstOrder = SISO | Model(
    k = 1.0,
    T = 1.0,
    x = Var(init=0.0),
    equations = :[
        der(x) = (k * u - x) / T
        y = x
    ]
)

# Integrator
Integrator = SISO | Model(
    k = 1.0,
    x = Var(init=0.0),
    equations = :[
        der(x) = k * u
        y = x
    ]
)

# VoltagePin = Model(
#     vRe = potential,
#     vIm = potential,
# )

VoltagePin = Model(
    v = Var(potential=true, start= [0.0, 0.0]),
)

"""
function ComplexMulti(vRe::Float64,vIm::Float64, Thetap::Float64)::Float64

y = abs(V * exp(J Theta))

"""
function ComplexMulti(vRe::Float64,vIm::Float64, Kp::Float64, Thetap::Float64)::Float64
    return Kp*abs(Complex(vRe,vIm) * exp(Thetap*1.0im))
end

using LinearAlgebra

#FIXME
#FIXME
#FIXME
function Betrag1(A::Vector{Float64}):: Float64
    return norm(A)
end

function Betrag1(A::SVector{2}):: Float64
    return norm(A)
end

GENvolts = Model(   # für IEEE exciter model ; Xl = 0, Kii = 0 -> siehe Beispiel Parameterliste ST8C
    Kp = 1.0,
    Thetap = 0.0,
    vT = VoltagePin,
    y = output,
    equations = :[
#         y = ComplexMulti(vT.vRe, vT.vIm, Kp, Thetap)
        y = Kp* Controller.Betrag1(vT.v)
    ]
)

"""
# IEEEpart = Model()
Partial model with transducer and sum
interface: Verror

"""
IEEEpart = Model(
    vT = VoltagePin,    # Complex GEN volts

    TR= 0.0226,     # Regulator input filter time constant
    vS = input,     # PSS in
    vREF = input,   # setpoint in
    vGEN = input,   # GEN volts
    iFD = input,    # field current
    y = output,     # EXC out

    transducer = FirstOrder | Map(T=:TR, k=1.0, x=Var(init=1.0)),

    equations = :[
        transducer.u = vGEN
        Verror = vREF+vS-transducer.y
    ]
)

PInonwindup = SISO | Model(    # NO state events !!
    Ki = 1.0,
    Kp = 1.0,
    A = -5.0,   # min
    B = 5.0,    # max
    limitA = Boolean | Map(init=false),
    limitB = Boolean | Map(init=false),
    Integr = Integrator | Map(k=1.0),
    equations = :[
        addy = u*Kp + Integr.y
        limitA = (addy<A)
        limitB = (addy>B)
        Integr.u = if  (limitA || limitB) ; 0.0; else; u*Ki; end;
        y = if limitB; B; elseif limitA; A; else; addy; end;
    ]
)

PInonwindup0 = SISO | Model(    # NO state events !!
    Kp = 1.0,
    Ki = 1.0,
    A = -5.0,   # min
    B = 5.0,    # max
    Integr = Integrator | Map(k=1.0),
    equations = :[
        Integr.u =  u*Ki
        y = u*Kp + Integr.y
    ]
)

"""
# PT1nonwindup = SISO | Model()

Keine state machine möglich.
### temp. Lösung
wenn obere Grenze erreicht (B) dann `Integr.u=1e-10` , damit Grenze nicht unterschritten wird. `Integr.y` kann sich in der Simulation etwas ändern, und damit könnte Grenze unterschritten werden.
"""
PT1nonwindup = SISO | Model(    # NO state events !!
    K = 1.0,
    T = 0.5,
    A = -5.0,   # min
    B = 5.0,    # max
    limitA = Boolean | Map(init=false),
    limitB = Boolean | Map(init=false),

    Integr = Integrator | Map(k=:(1.0/T)),
    equations = :[
        Integr.u = if (limitA || limitB); null; else; feedback; end;
        null = limitA ? -1e-10 : 1e-10
        feedbackPOS = feedback>0.0
        limitA = (Integr.y<A) && !feedbackPOS
        limitB = (Integr.y>B) && feedbackPOS
        feedback = u - Integr.y
        y = Integr.y * K
    ]
)

PT1nonwindup0 = SISO | Model(    # NO state events !!
    K= 1.0,
    T = 0.5,
    A = -5.0,   # min
    B = 5.0,    # max

    Integr = Integrator | Map(k=:(1.0/T)),
    equations = :[
        Integr.u = u - Integr.y
        y = Integr.y * K
    ]
)

"""
# ExcCONST = IEEEpart | Model()
Konstante Erregung. `y = vREF + vS`

"""
ExcCONST = IEEEpart | Model(
    vS = Var(),
    equations = :[
        y = vREF + vS
    ]
)

"""
# ExcST8C = IEEEpart | Model()
IEEE ST8C model

"""
ExcST8C = IEEEpart | Model(

    TR= 0.0226,         # Regulator input filter time constant
    Kp = 6.3,           # Potential circuit (voltage) gain coefficient
    Thetap = 0.0,       # Potential circuit phase angle (degrees)

    VPImin = 0.0,       # Minimum voltage regulator output
    VPImax = 1.8,       # Maximum voltage regulator output
    Kpr = 10.9,         # Voltage regulator proportional gain
    Kir = 14.29,        # Voltage regulator integral gain

    Kf = 0.452,         # Exciter field current feedback gain	// airgap base
    Tf = 0.005,         # Field current feedback time constant

    VAmin = -0.866,     # Minimum field current regulator output
    VAmax =  0.996,     # Maximum field current regulator output
    Kpa = 4.0,          # Field current regulator proportional gain
    Kia = 0.0,          # Field current regulator integral gain

    VRmin = -0.866,     # Minimum field current regulator output ??
    VRmax =  0.966,     # Maximum field current regulator output ??
    Ka = 1.0,           # Field current regulator proportional gain ??
    Ta = 0.0033,        # Controlled rectifier bridge equivalent time constant

    SW1 = true,         # Power source selector, A if true (from Gen voltage)

#TODO KEINE state events ==> VIEL SCHNELLER
    PI1 = PInonwindup | Map(A = :VPImin, B = :VPImax, Ki = :Kir, Kp = :Kpr),
    PIa = PInonwindup | Map(A = :VAmin,  B = :VAmax,  Ki = :Kia, Kp = :Kpa),
    PI1f = FirstOrder | Map(T=:Tf, k=:Kf),
    PT1a = PT1nonwindup | Map(A = :VRmin, B = :VRmax, K = :Ka, T = :Ta),
    genVolts = GENvolts | Map(Kp = :Kp, Thetap = :Thetap),

    equations = :[
        PI1.u = Verror
        PI1f.u = iFD
        PIa.u = PI1.y - PI1f.y
        PT1a.u = PIa.y
        y = if SW1; genVolts.y *PT1a.y; else; Kp*PT1a.y; end;
        ],

    connect = :[
         (vT, genVolts.vT)
    ]

)

end

ext: jl

max_stars_repo_name: johhell/EPSphasor

max_issues_repo_name: johhell/EPSphasor

max_forks_repo_name: johhell/EPSphasor


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 12/15/2014

module AbstractEncounterDBNImpl

export AbstractEncounterDBN

abstract AbstractEncounterDBN

end

module AbstractEncounterDBNInterfaces

export
       # common
       get,
       getInitialState

# common

function get()
end

function getInitialState()
end

end


text:
#
# module file
#

module MaybeEither

export
    Maybe, Just, None, isjust, isnone, eltype,
    Either, Left, Right, left, right, isleft, isright, lefttype, righttype, value

include("maybe.jl")
include("either.jl")

end # module


text:
using ImageFiltering, ImageCore, OffsetArrays
using Base.Test

@testset "cascade" begin
    a = rand(15)
    kern = OffsetArray(ones(3), -1:1)
    kern2 = OffsetArray(Float64[1,2,3,2,1], -2:2)  # the convolution of the above with self
    for border in ("replicate", "circular", "symmetric", "reflect", Fill(zero(eltype(a))))
        afc = imfilter(a, (kern, kern), border)
        af2 = imfilter(a, kern2, border)
        @test_approx_eq afc af2
    end

    a = rand(15, 15)
    kernx, kerny = OffsetArray(ones(3,1), -1:1, 0:0), OffsetArray(ones(1,3), 0:0, -1:1)
    c = Float64[1,2,3,2,1]
    kern2 = OffsetArray(c.*c', -2:2, -2:2)
    kern2x = OffsetArray(c.*ones(1,3), -2:2, -1:1)
    kern2y = OffsetArray(ones(3).*c',  -1:1, -2:2)
    for border in ("replicate", "circular", "symmetric", "reflect", Fill(zero(eltype(a))))
        afc = imfilter(a, (kernx, kerny, kernx, kerny), border)
        af2 = imfilter(a, kern2, border)
        @test_approx_eq afc af2
        afc = imfilter(a, (kernx, kernx, kerny, kerny), border)
        af2 = imfilter(a, kern2, border)
        @test_approx_eq afc af2
        afc = imfilter(a, (kernx, kernx, kerny), border)
        af2 = imfilter(a, kern2x, border)
        @test_approx_eq afc af2
        afc = imfilter(a, (kerny, kernx, kerny), border)
        af2 = imfilter(a, kern2y, border)
        @test_approx_eq afc af2
    end
end


text:
# To fix the broken jld
using JLD, Knet, ArgParse
include("src/header.jl")
const mfile = "src/new_model_feat3.jl"
include(mfile)

# old-style input
function fixjld(jldin, jldout, dfile)
    bundle = load(jldin)
    if haskey(bundle, "fnav")
        println("cp $jldin $jldout | it is fixed")
        return
    end

    model, optims      = bundle["allmodel"], bundle["optims"]
    wmodel, v          = bundle["wordmodel"], bundle["vocab"]
    #featdict, xposdict = bundle["featdict"], bundle["xposdict"]
    odict              = bundle["odict"]
    c, xposdict2, featdict2 = load_conllu3(dfile, v);

    fnav = createfnav(featdict2); # you need to save that dictionary for test-time usage
    wmodel2 = map2cpu(wmodel) # idiot
    JLD.save(jldout,
             "allmodel", model,    "optims", optims,
             "featdict", featdict2, "xposdict", xposdict2,
             "wordmodel", wmodel2, "vocab", v,
             "odict",      odict,  "fnav", fnav
             )
end

function main(args=ARGS)
    s = ArgParseSettings()
    @add_arg_table s begin
        ("--input";  help="Input jld")
        ("--output"; help="Output jld")
        ("--tfile"; help="To initialize train file")
    end
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    fixjld(o[:input], o[:output], o[:tfile])
end

!isinteractive() && main(ARGS)


text:
######################################################################
# Slightly advanced logging utility
# -----
# Licensed under Apache License 2.0
module Logging
import ExtraFun: curry

struct LogMessage
    message::String
    time::Float64
end
LogMessage(msg::AbstractString) = LogMessage(string(msg), time())

export LogCategory
struct LogCategory
    name::String
    default_duplex::Union{IO, Nothing}
    messages::Vector{LogMessage}
end
LogCategory(name::AbstractString = "uncategorized", default_duplex::Union{IO, Nothing} = stdout) = LogCategory(string(name), default_duplex, LogMessage[])
(cat::LogCategory)(args...; kwargs...) = log(cat, args...; kwargs...)

const Uncategorized    = LogCategory()
const DebugCategory    = LogCategory("Debug")
const InfoCategory     = LogCategory("Info")
const WarnCategory     = LogCategory("Warning")
const ErrorCategory    = LogCategory("Error", stderr)
const CriticalCategory = LogCategory("Critical", stderr)

export logdebug, loginfo, logwarn, logerror, logcritical
function log(cat::LogCategory, xs...; duplex::Union{IO, Nothing} = cat.default_duplex, joiner = ' ')
    msg = join(string.(xs), joiner)
    push!(cat.messages, LogMessage(msg))

    if duplex !== nothing
        write(duplex, msg * '\n')
    end
    nothing
end
log(args...; kwargs...) = log(Uncategorized, args...; kwargs...)

const logdebug    = curry(log, DebugCategory)
const loginfo     = curry(log, InfoCategory)
const logwarn     = curry(log, WarnCategory)
const logerror    = curry(log, ErrorCategory)
const logcritical = curry(log, CriticalCategory)

function Base.empty!(log::LogCategory; older_than::Real = Inf)
    filter!(msg->msg.time ≥ older_than, log.messages)
    return log
end

end # module Logging


text:
#!/usr/bin/env julia
include("commandUtils.jl")
activatePkg()
using ArgParse
using JLD
using Distributed
s = ArgParseSettings(description="Resonant Mode Search, runs algorithms 1 & 2 from Heider paper.")
@add_arg_table s begin
    "--simFile", "-s"
        help = "Simulation file to include, must set up boundary named simBoundary, can also edit or set k or other params to avoid the need to
                set in command arguments."
        required = true
    "--fieldFile", "-f"
        help = "File with calculated k's and field approximation as initial guess for x0."
    "--fieldIdx", "-i"
        help = "Index in field file to use for k and x0."
        arg_type = Int
        default = 1
    "--k0"
        help = "Search for resonances close to this k (can also set k0 in simFile)."
        arg_type = Complex{Float64}
    "--eps"
        help = "Epsilon value, stop iteration when k changes less than this value "
        arg_type = Float64
        default = .000001
    "--dx"
        help = "Total field dx"
        arg_type = Float64
        default = 0.01
    "--dy"
        help = "Total field dy"
        arg_type = Float64
        default = 0.01
    "--xRange"
        help = "X range (tuple)"
        default = "(-1.5, 1.5)"
    "--yRange"
        help = "Y range (tuple)"
        default = "(-1.5, 1.5)"
    "--procs", "-p"
        help = "Number of processes to use in parallel"
        arg_type = Int
        default = convert(Int, Sys.CPU_THREADS/2)
end
arguments = parse_args(s)

@debug arguments

addprocs(arguments["procs"])
@everywhere using BoundaryIntegralMethod

ϵ = arguments["eps"]
fieldFile = arguments["fieldFile"]
fieldIdx = arguments["fieldIdx"]

xRange = eval(Meta.parse(arguments["xRange"]))
yRange = eval(Meta.parse(arguments["yRange"]))
dx = arguments["dx"]
dy = arguments["dy"]

simFile = arguments["simFile"]
simName = importSim(simFile)
b = simBoundary

ks = nothing
xs = nothing
if (fieldFile != nothing)
  (ks, xs) = jldopen(fieldFile, "r") do file
    (read(file, "ks/full"), read(file, "xs/full"))
  end
end

if (ks != nothing)
  k0 = ks[fieldIdx]
  x0 = xs[fieldIdx,:]
end

n = convert(Int,size(x0,1)/4)

@info "Looking for resonance at k=$k0"

(kf, xf) = residualInverseIteration(10, ϵ, b, n, k0, x0)

@info "Final k=$kf"

kfForFile = string(real(kf),"_",imag(kf))
outFile = "output/field-$simName-$n-$kfForFile.h5"
mkpath("output")
boundaryOutFile = "output/boundaryField-$simName-$n-$kfForFile.h5"
jldopen(boundaryOutFile, "w") do file
  write(file, "k/full", kf)
  write(file, "k/Re", real(kf))
  write(file, "k/Im", imag(kf))
  write(file, "boundary/full", xf)
  write(file, "boundary/Re", real(xf))
  write(file, "boundary/Im", imag(xf))
end

exportTotalField(b, xf, kf; outFile=outFile, xRange=xRange, yRange=yRange, dx=dx, dy=dy)


text:
using Memoize
using GenericSchur
using JLD

"""
    greyencode(n::Integer) = n ⊻ (n >> 1)
returns grey encoding of given integer
"""
greyencode(n::Integer) = n ⊻ (n >> 1)

"""
    inverseM(N)
finds inverse of M matrix of size 2^(N-1)×2^(N-1) with elements M_{i,j} equal to
inner product of binary rep of i (classical) and j (grey code), i,j ∈
[1:2^(N-1)] utilizes algorithm from arxiv:quant-ph/0404089
"""
@memoize function inverseM(N)
    m = 2^(N - 1)
    M = zeros((m, m))
    for i in 1:m
        for j in 1:m
            M[i,j] = (-1)^count_ones((i - 1) .& greyencode(j - 1))
        end
    end
    return inv(M)
end

"""
    stateprep(u, N)

creates Array of CircuitGate objects preparing state `u` over `N` qubits.. uses
algorithm from arxiv:quant-ph/0410066
"""
function stateprep(u, N, n=1)
    cg = CircuitGate[]
    k = 2
    grey = zeros(Int, 2^(N - n))
    for i in 0:N - 2
        grey[2^i:2^(i + 1):end] .= i + n + 1
    end
    grey[end] = N
    θ = real.(atan.(-u[2:2:end] ./ u[1:2:end]) .* 2)
    θ[isnan.(θ)] .= -π
    M = inverseM(N - n + 1)
    θ = M * θ
    for i in 1:2^(N - n)
        push!(cg, circuit_gate((n,), RyGate(θ[i])))
        push!(cg, circuit_gate((n), X, grey[i]))
    end
    return cg
end

"""
    unitary2circuit(m::Matrix{Complex64,2}) where {N}

compiles quantum circuit from given unitary matrix. Implements algorithms from
ff10.1016/j.cpc.2019.107001f, 10.1103/PhysRevA.69.032315, arxiv.org:1003.5760,
arxiv:quant-ph/0404089

# Examples
```@example; setup = :(using RandomMatrices)
M = Stewart(ComplexF64, 4);
cgc = unitary2circuit(M)
```
"""
function unitary2circuit(m::AbstractMatrix{ComplexF64}, N::Union{Int, Nothing}=nothing, wires=nothing)
    m * m' ≈ I || error("Only unitary matrices can be compiled into a quantum circuit")
    s = size(m)

    if isnothing(N)
        N = intlog2(s[1])
    end

    if Diagonal(diag(m)) ≈ m
        return compilediag(diag(m), N, wires)
    end

    if size(m)[1] == 2
        return compile1qubit(m, wires)[1]
    end

    if size(m)[1] == 4
        return compile2qubit(m, N, wires)
    end

    cgc = CircuitGate[]
    QR, τ = qr_unblocked(m)
    cgc = compilediag(diag(QR), N, wires)

    Dg = fill(1.0 + 0.0im, (2^N))
    Dg[1] = -1
    Dgm = compilediag(Dg, N)

    for i in 2^(N) - 1:-1:1
        prepd = CircuitGate[]
        u = zeros(ComplexF64, (2^N))
        u[i + 1:end] = QR[i + 1:end, i]
        u[i] = 1
        u = u ./ norm(u)
        angles = exp.(im .* angle.(u))
        u = u ./ angles
        D = compilediag(angles, N)

        for j in 1:N - 1
            cg = stateprep(u[1:2^(j - 1):end], N, j)
            append!(prepd, cg)
            u = apply(u, cg)
        end
        θ = real(atan(-u[2^(N - 1) + 1] ./ u[1]) .* 2)
        if !isnan(θ)
            cg = circuit_gate((N,), RyGate(θ))
        else
            cg = circuit_gate((N,), RyGate(-π))
        end
        push!(prepd, cg)

        prep = reverse(adjoint.(prepd))
        append!(prep, D)
        prepend!(prepd, adjoint.(D))

        append!(cgc, prepd)
        append!(cgc, Dgm)
        append!(cgc, prep)
    end
    return cgc
end

"""
    qr_unblocked(m::Matrix{Complex64,2})

performs serial qr decomposition on matrix `m`
"""
function qr_unblocked!(m::AbstractMatrix{ComplexF64}, n=1::Integer)
    N = size(m)[1]
    τ = ComplexF64[]
    for i in n:N - 1
        b = norm(m[i,i])
        a = angle(m[i,i])
        push!(τ, 1 + b)
        m[i,i] = -exp(im * a)
        m[i + 1:N, i] =  m[i + 1:N, i] ./ (exp(im * a) * (1 + b))
        m[i + 1:N,i + 1:N] -= (m[i + 1:N, i] * transpose(m[i,i + 1:N]))
    end
    return τ
end

"""
    qr_unblocked(m::Matrix{Complex64,2})

performs serial qr decomposition on matrix `m`
"""
function qr_unblocked(m::AbstractMatrix{ComplexF64}, n=1::Integer)
    m = deepcopy(m)
    N = size(m)[1]
    τ = ComplexF64[]
    for i in n:N - 1
        b = norm(m[i,i])
        a = angle(m[i,i])
        push!(τ, 1 + b)
        m[i,i] = -exp(im * a)
        m[i + 1:N, i] =  m[i + 1:N, i] ./ (exp(im * a) * (1 + b))
        m[i + 1:N,i + 1:N] -= (m[i + 1:N, i] * transpose(m[i,i + 1:N]))
    end
    return m, τ
end

"""
    compile1qubit(m::AbstractMatrix{ComplexF64}, N, wires=nothing)

compiles a U(2) matrix into quantum gates
"""
function compile1qubit(m::AbstractMatrix{ComplexF64}, wires=nothing)
    if isnothing(wires)
        wires = [1]
    end
    cg = CircuitGate[]
    phase = sqrt(m[1,1]*m[2,2] - m[2,1]*m[1,2])
    m  = m / phase
    b = norm.(m)
    if (norm(b[1,1]) + norm(b[2,2])) < 1e-5
        θ2 = -imag(log(m[2,1]/m[1,2]))
        θ1 = 0
        ϕ = (real(m[1,2]/exp(-θ1/2*im)) < 0 ? 3π/2 : π/2) * 2
    elseif (norm(b[1,2]) + norm(b[2,1])) < 1e-5
        θ2 = -imag(log(m[1,1]/m[2,2]))
        θ1 = 0
        ϕ = (real(m[1,1]/exp(-θ1/2*im)) < 0 ? π : 0) * 2
    else
        ϕ = real(acos(sqrt(m[1,1] * m[2,2])) * 2)
        θ1 = imag(log(2im * m[2,1] * m[2,2] / sin(ϕ)))
        θ2 = imag(-log(2im * m[1,1] * m[2,1] / sin(ϕ)))
    end

    if !(norm(θ2) < 1e-5)
        append!(cg, circuit_gate.((wires[1],), [RzGate(θ2), RxGate(ϕ)]))
    else
        append!(cg, circuit_gate.((wires[1],), [RxGate(ϕ)]))
    end

    if !(norm(θ1) < 1e-5)
        append!(cg, circuit_gate.((wires[1],), [RzGate(θ1)]))
    end
    return cg, phase
end

iso_permutations = [
                    [[1,1], [2,2], [3,3], [4,4]],
                    [[2,1], [1,2], [4,3], [3,4]],
                    [[3,1], [4,2], [1,3], [2,4]],
                    [[4,1], [3,2], [2,3], [1,4]],
                    [[2,1], [1,2], [4,3], [3,4]],
                    [[1,1], [2,2], [3,3], [4,4]],
                    [[4,1], [3,2], [2,3], [1,4]],
                    [[3,1], [4,2], [1,3], [2,4]],
                    [[3,1], [4,2], [1,3], [2,4]],
                    [[4,1], [3,2], [2,3], [1,4]],
                    [[1,1], [2,2], [3,3], [4,4]],
                    [[2,1], [1,2], [4,3], [3,4]],
                    [[4,1], [3,2], [2,3], [1,4]],
                    [[3,1], [4,2], [1,3], [2,4]],
                    [[2,1], [1,2], [4,3], [3,4]],
                    [[1,1], [2,2], [3,3], [4,4]],
                    ]

iso_signs = [
            [ 1, 1, 1, 1],
            [ 1,-1, 1,-1],
            [ 1,-1,-1, 1],
            [ 1, 1,-1,-1],
            [ 1,-1,-1, 1],
            [-1,-1, 1, 1],
            [-1,-1,-1,-1],
            [ 1,-1, 1,-1],
            [ 1, 1,-1,-1],
            [ 1,-1,-1, 1],
            [-1, 1,-1, 1],
            [-1,-1,-1,-1],
            [ 1,-1, 1,-1],
            [-1,-1,-1,-1],
            [ 1, 1,-1,-1],
            [-1, 1, 1,-1]
            ]

function _get_permutation(m::AbstractMatrix{Float64}, ::Type{Val{I}}, ::Type{Val{J}}) where {I,J}
    val = 0
    for i in 1:4
        val += m[iso_permutations[4*(J-1)+I][i]...] * iso_signs[4*(J-1)+I][i]
    end
    val/4
end

function _get_col!(m::AbstractMatrix{Float64}, col::Vector{Float64}, col_id::Int)
    for i in 1:4
        col[i] = _get_permutation(m, Val{i}, Val{col_id})
    end
end

function _get_row!(m::AbstractMatrix{Float64}, row::Vector{Float64}, row_id::Int)
    for i in 1:4
        row[i] = _get_permutation(m, Val{row_id}, Val{i})
    end
end

"""
    decomposeSO4(m::AbstractMatrix{ComplexF64})

decomposes a SO(4) matrix into 2 SU(2) matrices, such that `m`∈SO(4), A,B ∈
SU(2) via isoclinic decomposition. m ≊ A ⊗ B under the magic basis homomorphism.
"""
function decomposeSO4(m::AbstractMatrix{Float64})
    size(m)[1] == size(m)[2] || error("decomposeSO4 only works on square matrices")
    size(m)[1] == 4 || error("decomposeSO4 only works on 4x4 matrices")
    det(m) ≈ 1 || error("matrix `m` is not in SO(4), determinant is not 1")

    col = zeros(Float64, 4)
    row = zeros(Float64, 4)

    col_id = 0
    row_id = 0

    for i in 4:-1:1
        _get_col!(m, col, i)
        if sum(norm.(col)) > 0.5
            col_id = i
            break
        end
    end

    for i in 4:-1:1
        _get_row!(m, row, i)
        if sum(norm.(row)) > 0.5
            row_id = i
            break
        end
    end

    col_val = sqrt(sum(col.^2))
    row_val = sqrt(sum(row.^2))

    for i in 0:3

        col_val = col_val * (-1)^(i)
        row_val = row_val * (-1)^(i÷2)

        a, b, c, d = col./col_val
        p, q, r, s = row./row_val

        if col[row_id]/col_val ≈ row_val && row[col_id]/row_val ≈ col_val
            @assert p^2 + q^2 + r^2 + s^2 ≈ 1
            @assert a^2 + b^2 + c^2 + d^2 ≈ 1
            A = [p + q * im s + r * im; -s + r * im p - q * im]
            B = [a + b * im -d + c * im; d + c * im a - b * im]
            return A, B
        end
    end
    error("Algorithm does not work")
end

function decomposeSO4(m::AbstractMatrix{ComplexF64})
    decomposeSO4(real.(m))
end

"""
    compile2qubit(m::AbstractMatrix{ComplexF64}, N, wires=nothing)

compiles an arbitrary U(4) matrix into a quantum circuit. Algorithm taken from
arxiv:quant-ph/0308006, arxiv:quant-ph/0211002
"""
function compile2qubit(m::AbstractMatrix{ComplexF64}, N, wires=nothing)
    cg = CircuitGate[]
    E = 1 / sqrt(2) .* [1 im 0 0; 0 0 im 1; 0 0 im -1; 1 -im 0 0]
    U = E' * m * E

    P2 = U * transpose(U)

    Diag, K_2 = eigen(P2)
    K_2 = real.(K_2)

    # ensure that eigenvectors for degenerate eigenvalues are orthogonal
    if !(det(K_2) ≈ 1) || !(det(K_2) ≈ -1)
        for i in 1:4
            for j in i+1:4
                if norm(dot(K_2[:,j], K_2[:, i])) > 1e-12
                    K_2[:, [i,j]] .= gramm_schmidt!(K_2[:, [i,j]])
                end
            end
        end
    end

    Diag[1] = Diag[1] / det(K_2)^2
    K_2[:,1] = K_2[:,1] * det(K_2)

    Diag = diagm(sqrt.(Diag))
    Diag[1] = Diag[1] * det(U) / det(Diag)

    P = K_2 * Diag * inv(K_2)
    K_1 = inv(P) * U
    K = inv(K_2) * K_1

    C, D = decomposeSO4(inv(K_2) * K_1)

    # println(norm(E'*kron(C,D)*E - inv(K_2) * K_1))
    append!(cg, compile1qubit(C, [2])[1])
    append!(cg, compile1qubit(D, [1])[1])

    append!(cg, [circuit_gate(2, X, 1),
                    circuit_gate((1,), HadamardGate()),
                    circuit_gate((2,), SdagGate()),
                    circuit_gate((1,), SdagGate())])
    append!(cg, compilediag(diag(Diag), N))
    append!(cg, [circuit_gate((1,), SGate()),
                    circuit_gate((2,), SGate()),
                    circuit_gate((1,), HadamardGate()),
                    circuit_gate(2, X, 1)])

    A, B = decomposeSO4(K_2)

    # println(norm(E'*kron(A,B)*E - K_2))
    append!(cg, compile1qubit(A, [2])[1])
    append!(cg, compile1qubit(B, [1])[1])

    return cg
end

"""
    fillψ!(d::Vector{ComplexF64}, ψ::Vector{Float64}, l::Integer)

fills empty vector `ψ`, such that  ψ = −im*[logχ1(d) log χ2(d) ··· logχl−1(U)],
where χn(d) = d[2n-1]*d[2n+2]/(d[2n]*d[2n+1]) per arxiv:quant-ph/0303039
"""
function fillψ!(d::AbstractVector{ComplexF64}, ψ::Vector{Float64}, l::Integer)
    for i in StepRange(1, 1, l - 1)
        ψ[i] = imag(log(d[2i - 1] * d[2i + 2] / (d[2i] * d[2i + 1])))
    end
end

"""
    ηcol(l::Integer, n::Integer)

calculates a single flip state for given integer `n` and a column length of `l`,
    where 0 <= n <= l. returns a vector of length l, where vec[n] = 1, vec[n+1]
    = -1. if n == 0, vec[1] = -1, if n==l, vec[l] = 1
# Examples
```julia-repl
julia> l = 4
julia> ηcol(4, 0)
[-1, 0, 0, 0]
julia> ηcol(4, 1)
[1, -1, 0, 0]
julia> ηcol(4, 2)
[0, 1, -1, 0]
julia> ηcol(4, 3)
[0, 0, 1, -1]
julia> ηcol(4, 4)
[0, 0, 0, 1]
```
"""
@memoize Dict function ηcol(l::Integer, n::Integer)
    vec = zeros(Integer, l)
    if n == 0
        vec[1] = -1
        return vec
    elseif n == l
    vec[l] = 1
        return vec
    end

    vec[n] = 1
    vec[n + 1] = -1
    vec
end

"""
    flip_state(m::Integer, l::Integer)

calculates the `m`th column of the η⊗ matrix by summing all flip states
associated with the binary representation of `m` over the span of [1,l-1]
algorithm taken from arxiv:quant-ph/0303039
"""
@memoize Dict function flip_state(m::Integer, l::Integer)
    # calculates flip states of a given integer `x`
    f = filter(x -> count_ones(x & m) % 2 == 1, 1:l - 1)
    sum(ηcol.((l - 1,), f))
end

"""
    svalue(i::Integer)

returns tuple of position of 1s in the binary representation of integer `i`
starting from least significant digit.
# Examples
```julia-repl
julia> svalue(2)
(2,)
julia> svalue(4)
(3,)
julia> svalue(7)
(1,2,3)
```
"""
function svalue(i::Integer)
    b = Int[]
    count = 1
    while i != 0
        if i & 1 == 1
            push!(b, count)
        end
        count += 1
i = i >> 1
    end
    Tuple(b)
end

"""
    compilediag(d::Vector{ComplexF64}, N, cg=nothing, j=0)

recursively compiles an arbitrary diagonal unitary matrix D_{2^N} ∈ C^(2^N×2^N)
into a quantum circuit. algorithm taken from arxiv:quant-ph/0303039 D_{2^N} can
be decomposed into D_{2^(N-1)} ∈ C^(2^(N-1)×2^(N-1)) ⊗ exp(i*ϕ)*RzGate(θ)
"""
function compilediag(d::AbstractVector{ComplexF64}, N, cg=nothing, j=0)
    cgs = CircuitGate[]
    # single qubit diagonal matrix can be split into a phase shift + RzGate
    if length(d) == 2
        logd = imag.(log.(d))
        β = (logd[1] - logd[2])
        if isnothing(cg)
            cg = CircuitGate[]
        end
        if !(β ≈ 0)
            push!(cg, circuit_gate((j + 1,), RzGate(-β)))
        end
        return cg
    end
    l = 2^(N - j - 1)
    ψ = zeros(Float64, l - 1)
    α = zeros(Float64, l)

    fillψ!(d, ψ, l)

    ηplus = zeros(Float64, (l - 1, l - 1))
    for i in StepRange(1, 1, l - 1)
    ηplus[:,i] = flip_state(greyencode(i), l)
    end

    ηplus = inv(ηplus)
    α[2:end] = -0.5 * ηplus * ψ

    gatewires = svalue(greyencode(1))
    oldwires = ()

    for i in StepRange(2, 1, l)
        rz = [exp(-im * α[i] / 2), exp(im * α[i] / 2)]
        for _ in 1:N - j - 1
            rz = append!(rz, rz)
        end
        s = filter(x -> count_ones(x & greyencode(i - 1)) % 2 == 1, 1:l - 1)
        s = append!(2 .* s .+ 1, 2 .* s .+ 2)

        for k in symdiff(gatewires, oldwires)
            push!(cgs, circuit_gate((j + 1), X, k + 1 + j))
        end
        oldwires = gatewires
        gatewires = svalue(greyencode(i))
        if !(α[i]≈0)
            push!(cgs, circuit_gate((j + 1,), RzGate(α[i])))
        end
        rz[s] .= 1 ./ rz[s]
        d = d ./ rz
    end

    for k in reverse(oldwires)
        push!(cgs, circuit_gate(j + 1, X, k + 1 + j))
    end

    if !(d[1] ≈ d[2])
        logd = imag.(log.(d))
        α[1] = logd[1] - logd[2]
        if !(α[1] ≈ 0)
            pushfirst!(cgs, circuit_gate((j + 1,), RzGate(-α[1])))
        end
        d[1:2:end] = d[1:2:end] ./ exp(im * α[1] / 2)
    end
    if isnothing(cg)
        return compilediag(d[1:2:end], N, cgs, j + 1)
    end
    return compilediag(d[1:2:end], N, append!(cg, cgs), j + 1)
end


text:
module QuantumAlgebra

export scal,param,a,adag,f,fdag,OpSumAnalytic,ExpVal,Corr
export σx,σy,σz,σp,σm,comm,latex
export Avac,vacA,vacExpVal
export CorrOrExp,ascorr
export @Pr_str, @Pc_str, ∑
#export σ,preftuple,exptuple,optuple,prodtuples,sumtuples

using Printf

# we will want to overload these operators and functions for our custom types
import Base: ==, ≈, *, +, -, isless, length, adjoint, print, zero, one

# define for σx, σy, σz
@enum SpatialIndex x=1 y=2 z=3
SpatialIndex(a::SpatialIndex) = a

# we sometimes want to sort integers and symbols together, with integers coming first
sortsentinel(x::Integer) = (1,x)
sortsentinel(x::Symbol) = (2,x)

const OpIndex = Union{Int,Symbol}
const OpIndices = NTuple{N,OpIndex} where N

abstract type Operator end
abstract type Scalar <: Operator; end
"Represent a scalar value (i.e., a number)"
struct scal{T<:Number} <: Scalar; v::T; end
"`param(g,(:i,:j),'n')`: represent a scalar named parameter ``g_{i,j}``. state can be purely real (`'r'`), not conjugated (`'n'`), or conjugated (`'c'`)"
struct param{T<:OpIndices} <: Scalar
    name::Symbol
    state::Char
    inds::T
    function param(name,state::Char,inds::OpIndex...)
        state in ('n','r','c') || error("state has to be one of n,r,c")
        new{typeof(inds)}(name,state,inds)
    end
    param(name,inds::OpIndex...) = param(name,'n',inds...)
    param(name,state::Char,inds::Tuple) = param(name,state,inds...)
    param(name,inds::Tuple) = param(name,'n',inds...)
end
struct δ <: Scalar
    inds::Tuple{OpIndex,OpIndex}
    δ(inds) = δ(inds...)
    function δ(iA::OpIndex,iB::OpIndex)
        # sort indices
        if iA == iB
            scal(1)
        elseif iA isa Integer && iB isa Integer
            # e.g., δ_1,3 = 0 (integers are not symbolic!)
            scal(0)
        elseif sortsentinel(iB) < sortsentinel(iA)
            new((iB,iA))
        else
            new((iA,iB))
        end
    end
end
function δ(Ainds::OpIndices,Binds::OpIndices)::Operator
    if length(Ainds) == length(Binds)
        prod(collect(Operator,δ.(Ainds,Binds)))
    else
        scal(0)
    end
end

abstract type BaseOperator <: Operator; end
for (op,desc,sym) in (
    (:a,   "bosonic annihilation","a"),
    (:adag,"bosonic creation","a^†"),
    (:f,   "fermionic annihilation","f"),
    (:fdag,"fermionic creation","f^†"),
    (:σminus,"TLS annihilation","σ^-"),
    (:σplus,"TLS creation","σ^+"))
    @eval begin
        "`$($op)(inds)`: represent $($desc) operator ``$($sym)_{inds}``"
        struct $op{T<:OpIndices} <: BaseOperator
            inds::T
            $op(inds::OpIndex...) = new{typeof(inds)}(inds)
            $op(inds::Tuple) = $op(inds...)
        end
    end
end

"`σ(a,n)`: represent Pauli matrix ``σ_{a,n}`` for two-level system (TLS) ``n``, where ``a ∈ \\{x,y,z\\}`` or ``\\{1,2,3\\}`` is the type of Pauli matrix."
struct σ{T<:OpIndices} <: BaseOperator
    a::SpatialIndex
    inds::T
    σ(a,inds::OpIndex...) = new{typeof(inds)}(SpatialIndex(a),inds)
    σ(a,inds::Tuple) = σ(a,inds...)
end

zero(::Type{<:Operator}) = scal(0)
zero(::Operator) = scal(0)
one(::Type{<:Operator}) = scal(1)
one(::Operator) = scal(1)

using_σpm = false
use_σpm(t::Bool=true) = eval(:( using_σpm = $t; nothing ))
use_σxyz() = use_σpm(false)

"`σp(n)`: construct ``σ^+_n = \\frac12 σ_{x,n} + \\frac{i}{2} σ_{y,n}``"
σp(n...) = using_σpm ? σplus(n...) : scal(1//2)*σ(x,n...) + scal(1im//2)*σ(y,n...)
"`σm(n)`: construct ``σ^-_n = \\frac12 σ_{x,n} - \\frac{i}{2} σ_{y,n}``"
σm(n...) = using_σpm ? σminus(n...) : scal(1//2)*σ(x,n...) - scal(1im//2)*σ(y,n...)
"`σx(n)`: construct ``σ_{x,n}``"
σx(n...) = using_σpm ? σminus(n...) + σplus(n...) : σ(x,n...)
"`σy(n)`: construct ``σ_{y,n}``"
σy(n...) = using_σpm ? scal(1im) * (σminus(n...) - σplus(n...)) : σ(y,n...)
"`σz(n)`: construct ``σ_{z,n}``"
σz(n...) = using_σpm ? scal(2)*σplus(n...)*σminus(n...) - scal(1) : σ(z,n...)

struct OpProd <: Operator; A::Operator; B::Operator; end
struct OpSum  <: Operator; A::Operator; B::Operator; end

"`OpSumAnalytic(i::Symbol,A::Operator)` or `∑(i,A)`: represent ``\\sum_{i} A``, with all possible values of ``i`` assumed to be included"
struct OpSumAnalytic <: Operator
    ind::Symbol
    A::Operator
    OpSumAnalytic(ind::Symbol,A::scal) = A.v == 0 ? A : new(ind,A)
    OpSumAnalytic(ind::Symbol,A::Operator) = begin
        if A isa OpProd && A.A isa scal
            A.A*OpSumAnalytic(ind,A.B)
        elseif A isa OpSumAnalytic && A.ind < ind
            # reorder nested sums by indices
            OpSumAnalytic(A.ind,OpSumAnalytic(ind,A.A))
        else
            for t in preftuple(A)
                # if there is a δ in the expression and it has the sum index, the sum disappears
                if t isa δ && ind in t.inds
                    irepl = ind==t.inds[1] ? t.inds[2] : t.inds[1]
                    return replace_index(A,ind,irepl)
                end
            end
            return new(ind,A)
        end
    end
    OpSumAnalytic(ind::Symbol,A::OpSum) = OpSumAnalytic(ind,A.A) + OpSumAnalytic(ind,A.B)
end

const ∑ = OpSumAnalytic

"`ExpVal(A::Operator)`: represent expectation value ``⟨A⟩``"
struct ExpVal <: Scalar
    A::Operator
    ExpVal(A::Scalar) = A
    ExpVal(A::OpSum) = ExpVal(A.A) + ExpVal(A.B)
    ExpVal(A::OpProd) = A.A isa Scalar ? A.A*ExpVal(A.B) : new(A)
    ExpVal(A::OpSumAnalytic) = OpSumAnalytic(A.ind,ExpVal(A.A))
    ExpVal(A::Operator) = new(A)
end

"""
    Corr(A::Operator)

Represents correlations of operator ``A``. ``A`` should be a product for this to
make sense, in which case ``⟨AB⟩_c = ⟨AB⟩ - ⟨A⟩⟨B⟩``, with corresponding
extensions for products of more operators.

See also: [`ascorr`](@ref)
"""
struct Corr <: Scalar
    A::Operator
    Corr(A::Scalar) = A
    Corr(A::OpSum) = Corr(A.A) + Corr(A.B)
    Corr(A::OpProd) = A.A isa Scalar ? A.A*Corr(A.B) : new(A)
    Corr(A::Operator) = new(A)
end

==(A::scal,B::scal) = A.v == B.v
==(A::OpProd,B::OpProd) = A.A == B.A && A.B == B.B
==(A::OpSum, B::OpSum)  = A.A == B.A && A.B == B.B
==(A::OpSumAnalytic, B::OpSumAnalytic) = A.A == B.A && A.ind == B.ind
==(A::ExpVal, B::ExpVal) = A.A == B.A
==(A::Corr, B::Corr) = A.A == B.A
==(A::param, B::param) = (A.name,A.inds,A.state) == (B.name,B.inds,B.state)

≈(A::Operator,B::Operator) = A == B
≈(A::scal,B::scal) = A.v ≈ B.v
≈(A::OpProd,B::OpProd) = A.A ≈ B.A && A.B ≈ B.B
≈(A::OpSum, B::OpSum)  = A.A ≈ B.A && A.B ≈ B.B
≈(A::OpSumAnalytic, B::OpSumAnalytic) = A.A ≈ B.A && A.ind == B.ind
≈(A::ExpVal, B::ExpVal) = A.A ≈ B.A
≈(A::Corr, B::Corr) = A.A ≈ B.A

prodtuple(A::Operator) = (A,)
prodtuple(A::OpProd) = (prodtuple(A.A)...,prodtuple(A.B)...)

# make a tuple from a product, containing only either prefactor types, expectation value types, or operators
for (name,types) in [(:pref,(scal,param,δ)),(:exp,(ExpVal,Corr)),(:op,(adag,a,f,fdag,σ,σplus,σminus))]
    name = Symbol(name,:tuple)
    types = Union{types...}
    @eval $name(A::$types) = (A,)
    @eval $name(A::Operator) = ()
    @eval $name(A::OpSum) = error("cannot get $($name) for OpSum!")
    @eval $name(A::OpSumAnalytic) = ($name(A.A)...,)
    @eval $name(A::OpProd) = ($name(A.A)...,$name(A.B)...)
end
prodtuples(A::Operator) = (preftuple(A), exptuple(A), optuple(A))

sumtuple(A::Operator) = (A,)
sumtuple(A::OpSum) = (sumtuple(A.A)...,sumtuple(A.B)...)

OpOrder = (scal,δ,param,ExpVal,Corr,adag,a,fdag,f,σplus,σminus,σ,OpProd,OpSumAnalytic,OpSum)
for (ii,op1) in enumerate(OpOrder)
    for op2 in OpOrder[ii+1:end]
        @eval isless(::$op1,::$op2) = true
        @eval isless(::$op2,::$op1) = false
    end
end
isless(A::scal,B::scal) = (real(A.v),imag(A.v)) < (real(B.v),imag(B.v))

# do not order parameters by whether they are purely real, or conjugated or not
isless(A::param,B::param) = (A.name,sortsentinel.(A.inds)) < (B.name,sortsentinel.(B.inds))
isless(A::σ,B::σ) = (sortsentinel.(A.inds),A.a) < (sortsentinel.(B.inds),B.a)
for op in (δ,a,adag,f,fdag,σminus,σplus)
    @eval isless(A::$op,B::$op) = sortsentinel.(A.inds) < sortsentinel.(B.inds)
end
for op in (ExpVal,Corr,OpSumAnalytic)
    @eval isless(A::$op,B::$op) = A.A < B.A
end
function isless(A::OpProd,B::OpProd)
    # order operator products first by number of operators (also within expectation values),
    # then by operators, expectation values and correlations, then by reversed prefactors (to order params, not scalars). Again use tuples to write this easily
    Aprefs,Aexps,Aops = prodtuples(A)
    Bprefs,Bexps,Bops = prodtuples(B)
    (length(A),Aops,Aexps,reverse(Aprefs)) < (length(B),Bops,Bexps,reverse(Bprefs))
end

# prefactors do not count for length calculation
for op in [scal,δ,param]
    @eval length(::$op) = 0
end
for op in [adag,a,fdag,f,σ,σminus,σplus]
    @eval length(::$op) = 1
end
for op in [ExpVal,Corr,OpSumAnalytic]
    @eval length(A::$op) = length(A.A)
end
length(A::OpProd) = length(A.A) + length(A.B)

*(A::Operator) = A # needed to make prod((A,)) work
*(A::OpProd,B::Operator) = A.A*(A.B*B)
*(A::OpSum,B::Operator)      = A.A*B + A.B*B
*(A::OpSum,B::OpSum)         = A.A*B + A.B*B # resolve ambiguity
*(A::OpSum,B::OpSumAnalytic) = A.A*B + A.B*B # resolve ambiguity
*(A::Operator,     B::OpSum) = A*B.A + A*B.B
*(A::OpProd,       B::OpSum) = A*B.A + A*B.B # resolve ambiguity
*(A::OpSumAnalytic,B::OpSum) = A*B.A + A*B.B # resolve ambiguity

# allow addition, substraction, and multiplication with a number x by promoting it to scal(x) operator
*(x::Number,A::Operator) = scal(x)*A
*(A::Operator,x::Number) = scal(x)*A
+(x::Number,A::Operator) = scal(x)+A
+(A::Operator,x::Number) = scal(x)+A
-(x::Number,A::Operator) = scal(x)-A
-(A::Operator,x::Number) = A-scal(x)

function *(A::Operator,B::Operator)::Operator
    if A isa scal && A.v==0
        A
    elseif A isa scal && A.v==1
        B
    elseif A isa scal && B isa scal
        scal(A.v*B.v)
    elseif A isa scal && B isa OpProd && B.A isa scal
        scal(A.v*B.A.v)*B.B
    elseif B isa OpProd && (A*B.A) != OpProd(A,B.A)
        # if A*B.A is not ordered as we want, evaluate A*B.A first
        (A*B.A)*B.B
    elseif A isa δ && A.inds[2] in indextuple(B)
        # if second index of δ_iA,iB shows up on RHS, replace by iA
        # (relies on indices in δ being sorted)
        A * replace_index(B,A.inds[2],A.inds[1])
    elseif A isa σ && B isa σ && A.inds == B.inds
        # σa σb = δab + i ϵabc σc = δab + 1/2 [σa,σb]
        A.a==B.a ? scal(1) : scal(1//2)*comm(A,B)
    elseif A isa Union{σplus,σminus,f,fdag} && typeof(B)==typeof(A) && A.inds == B.inds
        scal(0)
    elseif A>B
        B*A + comm(A,B)
    else
        OpProd(A,B)
    end
end

# these are the only overrides for minus, everything else falls back to +
-(A::Operator) = scal(-1)*A
-(A::Operator,B::Operator) = A + scal(-1)*B

separate_prefac(A::Operator) = (1,A)
separate_prefac(A::OpProd) = A.A isa scal ? (A.A.v,A.B) : (1,A)
function combinable(A::Operator,B::Operator)
    sA, oA = separate_prefac(A)
    sB, oB = separate_prefac(B)
    if oA==oB
        true, scal(sA+sB)*oA
    else
        false, scal(0)
    end
end

+(A::OpSum,B::Operator) = A.A + (A.B + B)
function +(A::Operator,B::Operator)::Operator
    A isa scal && A.v == 0 && return B
    A isa scal && B isa scal && return scal(A.v+B.v)
    A isa scal && B isa OpSum && B.A isa scal && return scal(A.v+B.A.v) + B.B

    # ignore scalars in ordering
    sA, oA = separate_prefac(A)
    if B isa OpSum
        sB, oB = separate_prefac(B.A)
        oA>oB && return B.A + (A + B.B)
    else
        sB, oB = separate_prefac(B)
        oA>oB && return B + A
    end

    flag, S = combinable(A,B)
    flag && return S

    if B isa OpSum
        flag, S = combinable(A,B.A)
        flag && return S + B.B
    end
    return OpSum(A,B)
end

adjoint(A::scal) = scal(conj(A.v))
adjoint(A::param) = A.state=='r' ? A : param(A.name,A.state=='n' ? 'c' : 'n',A.inds)
adjoint(A::a) = adag(A.inds)
adjoint(A::adag) = a(A.inds)
adjoint(A::f) = fdag(A.inds)
adjoint(A::fdag) = f(A.inds)
adjoint(A::σminus) = σplus(A.inds)
adjoint(A::σplus) = σminus(A.inds)
adjoint(A::σ) = A
adjoint(A::OpSum) = A.A' + A.B'
adjoint(A::OpProd) = A.B' * A.A'
adjoint(A::ExpVal) = ExpVal(A.A')
adjoint(A::Corr) = Corr(A.A')
adjoint(A::OpSumAnalytic) = OpSumAnalytic(A.ind,A.A')

comm(A::Scalar,B::Operator) = scal(0)
comm(A::Scalar,B::OpProd)   = scal(0)
comm(A::Scalar,B::OpSum)    = scal(0)
comm(A::Scalar,B::OpSumAnalytic) = scal(0)
comm(A::Operator,B::Scalar) = scal(0)
comm(A::OpSum,   B::Scalar) = scal(0)
comm(A::OpProd,  B::Scalar) = scal(0)
comm(A::OpSumAnalytic,B::Scalar) = scal(0)
comm(A::Scalar,B::Scalar)   = scal(0)
comm(A::OpProd,B::Operator) = comm(A.A,B)*A.B + A.A*comm(A.B,B)
comm(A::OpProd,B::OpSumAnalytic) = comm(A.A,B)*A.B + A.A*comm(A.B,B)
comm(A::OpProd,B::OpProd)   = comm(A.A,B)*A.B + A.A*comm(A.B,B)
comm(A::Operator,B::OpProd) = comm(A,B.A)*B.B + B.A*comm(A,B.B)
comm(A::OpSumAnalytic, B::OpProd) = comm(A,B.A)*B.B + B.A*comm(A,B.B)
comm(A::OpSum,B::Operator)  = comm(A.A,B) + comm(A.B,B)
comm(A::OpSum,B::OpSum)     = comm(A.A,B) + comm(A.B,B)
comm(A::OpSum,B::OpProd)    = comm(A.A,B) + comm(A.B,B)
comm(A::Operator,B::OpSum)  = comm(A,B.A) + comm(A,B.B)
comm(A::OpProd,  B::OpSum)  = comm(A,B.A) + comm(A,B.B)
comm(A::OpSumAnalytic, B::OpSum) = comm(A,B.A) + comm(A,B.B)
# different types of operators commute
commgroups = (Union{a,adag},Union{f,fdag},Union{σ,σminus,σplus})
for (ii,op) in enumerate(commgroups)
    for op2 in commgroups[ii+1:end]
        @eval comm(A::$op,B::$op2) = scal(0)
        @eval comm(A::$op2,B::$op) = scal(0)
    end
end
# these operators commute with themselves
for op in (a,adag,σplus,σminus)
    @eval comm(A::$op,B::$op) = scal(0)
end
comm(A::a,B::adag) = δ(A.inds,B.inds)
comm(A::adag,B::a) = -comm(B,A)
comm(A::σplus,B::σminus) = δ(A.inds,B.inds)*σz(A.inds)
comm(A::σminus,B::σplus) = -comm(B,A)

# {f_i, fdag_j} = f_i fdag_j + fdag_j f_i = δ_{i,j}
# f_i fdag_j = δ_{i,j} - fdag_j f_i
# [f_i, fdag_j] = f_i fdag_j - fdag_j f_i = δ_{i,j} - fdag_j f_i - fdag_j f_i = δ_{i,j} - 2 fdag_j f_i
comm(A::f,B::fdag) = δ(A.inds,B.inds) - 2*B*A
comm(A::fdag,B::f) = -comm(B,A)
# {f_i,f_j} = f_i f_j + f_j f_i = 0
# we assume that if we need the commutator, it's because of exchanging order
# [f_i,f_j] = f_i f_j - f_j f_i = -f_j f_i - f_j f_i
comm(A::f,B::f) = -2*B*A
comm(A::fdag,B::fdag) = -2*B*A

# levicivita_lut[a,b] contains the Levi-Cevita symbol ϵ_abc
# for c=6-a-b, i.e, when a,b,c is a permutation of 1,2,3
const levicivita_lut = [0 1 -1; -1 0 1; 1 -1 0]

function comm(A::σ,B::σ)
    if A.a == B.a
        scal(0)
    else
        a = Int(A.a)
        b = Int(B.a)
        # a+b+c == 6 (since a,b,c is a permutation of 1,2,3)
        c = 6 - a - b
        s = levicivita_lut[a,b]
        δ(A.inds,B.inds)*scal(2im*s)*σ(c,A.inds)
    end
end

replace_index(A::scal,iold,inew) = A
replace_index(A::param,iold,inew) = param(A.name,A.state,(n-> n==iold ? inew : n).(A.inds))
replace_index(A::ExpVal,iold,inew) = ExpVal(replace_index(A.A,iold,inew))
replace_index(A::Corr,iold,inew) = Corr(replace_index(A.A,iold,inew))
for op in (δ,a,adag,f,fdag,σminus,σplus)
    @eval replace_index(A::$op,iold,inew) = $op((n-> n==iold ? inew : n).(A.inds))
end
replace_index(A::σ,iold,inew) = σ(A.a,(n-> n==iold ? inew : n).(A.inds))
replace_index(A::OpProd,iold,inew) = replace_index(A.A,iold,inew)*replace_index(A.B,iold,inew)
replace_index(A::OpSum,iold,inew) = replace_index(A.A,iold,inew) + replace_index(A.B,iold,inew)
replace_index(A::OpSumAnalytic,iold,inew) = begin
    (A.ind==iold || A.ind==inew) && error("replace_index in OpSumAnalytic cannot have iold ($iold) or inew ($inew) be the same as the sum index ($(A.ind))!")
    tmp = replace_index(A.A,iold,inew)
    # we have to be careful here - when replacing an index inside the expression, there might be reorderings that
    # have be done with the same index explicitly to make sure that the commutator shows up
    # so take out the term with the same index explicitly, and redo that one term
    OpSumAnalytic(A.ind,tmp) - replace_index(tmp,A.ind,inew) + replace_index(replace_index(A.A,A.ind,iold),iold,inew)
end
replace_index(A::Operator,reps) = begin
    Anew = A
    for (iold,inew) in reps
        Anew = replace_index(Anew,iold,inew)
    end
    Anew
end

function exchange_inds(A::Operator,i1,i2)
    ss = gensym()
    A = replace_index(A,i1,ss)
    A = replace_index(A,i2,i1)
    A = replace_index(A,ss,i2)
    A
end

distribute_indices!(inds,A::scal) = A
distribute_indices!(inds,A::param) = param(A.name,A.state,(popfirst!(inds) for _ in A.inds)...)
distribute_indices!(inds,A::ExpVal) = ExpVal(distribute_indices!(inds,A.A))
distribute_indices!(inds,A::Corr) = Corr(distribute_indices!(inds,A.A))
for op in (a,adag,f,fdag,σminus,σplus)
    @eval distribute_indices!(inds,A::$op) = $op((popfirst!(inds) for _ in A.inds)...)
end
distribute_indices!(inds,A::σ) = σ(A.a,(popfirst!(inds) for _ in A.inds)...)
distribute_indices!(inds,A::OpProd) = distribute_indices!(inds,A.A)*distribute_indices!(inds,A.B)
distribute_indices!(inds,A::OpSum) = distribute_indices!(inds,A.A) + distribute_indices!(inds,A.B)
# on purpose do not define this for OpSumAnalytic or δ

"sum indices have no semantic meaning, so rename them in case they happen to occur in the other expression"
function ensure_compatible_sumind(S::OpSumAnalytic,A::Operator)
    Ainds = indextuple(A)
    if S.ind in Ainds
        oldinds = Set{OpIndex}((Ainds...,indextuple(S)...))
        m = match(r"(.*)_([0-9]+)",string(S.ind))
        indstem, ii = (m === nothing) ? (string(S.ind), 1) : (m.captures[1], 1+parse(Int,m.captures[2]))
        while (newind = Symbol(indstem,:_,ii)) in oldinds
            ii += 1
        end
        OpSumAnalytic(newind,replace_index(S.A,S.ind,newind))
    else
        S
    end
end

# when multiplying (or commuting) with an operator with an index, take into account that the term in the sum with equal index has to be treated specially
*(A::Union{param,ExpVal,Corr,BaseOperator},B::OpSumAnalytic) = (B = ensure_compatible_sumind(B,A); OpSumAnalytic(B.ind,A*B.A))
# no need to check indices here since we just dispatch to another routine
*(A::OpSumAnalytic,B::OpProd) = (A*B.A)*B.B
*(A::OpSumAnalytic,B::Operator) = (A = ensure_compatible_sumind(A,B); OpSumAnalytic(A.ind,A.A*B))

comm(A::Union{param,ExpVal,Corr,BaseOperator},B::OpSumAnalytic) = (B = ensure_compatible_sumind(B,A); OpSumAnalytic(B.ind,comm(A,B.A)))
comm(A::OpSumAnalytic,B::OpSumAnalytic) = (B = ensure_compatible_sumind(B,A); OpSumAnalytic(B.ind,comm(A,B.A)))
comm(A::OpSumAnalytic,B::Operator) = -comm(B,A)

print(io::IO,A::a) = print(io,"a($(A.inds...))")
print(io::IO,A::adag) = print(io,"a†($(A.inds...))")
print(io::IO,A::f) = print(io,"f($(A.inds...))")
print(io::IO,A::fdag) = print(io,"f†($(A.inds...))")
print(io::IO,A::σ) = print(io,"σ$(A.a)($(A.inds...))")
print(io::IO,A::σminus) = print(io,"σ-($(A.inds...))")
print(io::IO,A::σplus) = print(io,"σ+($(A.inds...))")
print(io::IO,A::scal) = print(io,"$(A.v)")
print(io::IO,A::param) = print(io,string(A.name, A.state=='c' ? "'" : "",length(A.inds)==0 ? "" : "($(A.inds...))"))
print(io::IO,A::OpProd) = print(io,"$(A.A) $(A.B)")
print(io::IO,A::OpSum) = print(io,"$(A.A) + $(A.B)")
print(io::IO,A::ExpVal) = print(io,"⟨$(A.A)⟩")
print(io::IO,A::Corr) = print(io,"⟨$(A.A)⟩c")
print(io::IO,A::OpSumAnalytic) = print(io,"∑_$(A.ind) $(A.A)")

mystring(x::Number) = @sprintf "%g" x
mystring(x::Rational) = denominator(x)==1 ? "$(numerator(x))" : "\\frac{$(numerator(x))}{$(denominator(x))}"
mystring(x::Complex) = @sprintf "(%g%+gi)" real(x) imag(x)
mystring(x::Complex{Rational{T}}) where T = @sprintf "\\left(%s%s%si\\right)" mystring(real(x)) (imag(x)>=0 ? "+" : "-") mystring(abs(imag(x)))

Base.show(io::IO, ::MIME"text/latex", A::Operator) = print(io,"\$",latex(A),"\$")
latex(A::σ) = string("\\sigma_{$(A.a)",length(A.inds)>0 ? ",$(A.inds...)}" : "}")
latexindstr(inds::OpIndices) = length(inds)==0 ? "" : "_{$(inds...)}"
latex(A::δ) = "δ" * latexindstr(A.inds)
latex(A::a) = "a" * latexindstr(A.inds)
latex(A::adag) = "a$(latexindstr(A.inds))^\\dagger"
latex(A::f) = "f" * latexindstr(A.inds)
latex(A::fdag) = "f$(latexindstr(A.inds))^\\dagger"
latex(A::σminus) = "\\sigma^-" * latexindstr(A.inds)
latex(A::σplus) = "\\sigma^+" * latexindstr(A.inds)
latex(A::scal) = imag(A.v)==0 ? mystring(real(A.v)) : (real(A.v)==0 ? mystring(imag(A.v))*"i" : mystring(A.v))
latex(A::param)  = string(A.name, latexindstr(A.inds), A.state=='c' ? "^*" : "")
latex(A::OpProd) = string(latex(A.A)," ",latex(A.B))
latex(A::OpSum) = string(latex(A.A)," + ",latex(A.B))
latex(A::ExpVal) = "\\langle $(latex(A.A)) \\rangle"
latex(A::Corr) = "\\langle $(latex(A.A)) \\rangle_{c}"
latex(A::OpSumAnalytic) = string("\\sum_{$(A.ind)}",latex(A.A))

indextuple(A::scal)::OpIndices = ()
indextuple(A::Union{param,δ,a,adag,f,fdag,σ,σminus,σplus})::OpIndices = A.inds
indextuple(A::Union{OpProd,OpSum})::OpIndices = (indextuple(A.A)...,indextuple(A.B)...)
indextuple(A::Union{ExpVal,Corr})::OpIndices = indextuple(A.A)
indextuple(A::OpSumAnalytic)::OpIndices = (A.ind,indextuple(A.A)...)
indexset(A) = Set{OpIndex}(indextuple(A))
sumindextuple(A::Operator)::OpIndices = ()
sumindextuple(A::Union{OpProd,OpSum})::OpIndices = (sumindextuple(A.A)...,sumindextuple(A.B)...)
sumindextuple(A::Union{ExpVal,Corr})::OpIndices = sumindextuple(A.A)
sumindextuple(A::OpSumAnalytic)::OpIndices = (A.ind,sumindextuple(A.A)...)
sumindexset(A) = Set{OpIndex}(sumindextuple(A))

"`extindices(A::Operator)` return externally visible indices of an expression"
extindices(A::Operator) = [ind for ind in indextuple(A) if !in(ind,sumindextuple(A))]

"`symmetric_index_nums(A::Operator)` return sequence of numbers of exchange-symmetric indices"
function symmetric_index_nums(A::Operator)
    inds = extindices(A)
    Nsyms = [1]
    for ii=2:length(inds)
        if A == exchange_inds(A,inds[ii-1],inds[ii])
            Nsyms[end] += 1
        else
            push!(Nsyms,1)
        end
    end
    Nsyms
end

"""
    ascorr(expr::Operator)

Take an expression `expr=A B C + D E...` and write its expectation value in
terms of single-body expectation values ``⟨A⟩, ⟨B⟩, \\ldots``, and many-body
correlations ``⟨AB⟩_c, ⟨ABC⟩_c``, etc. Currently, up to fourth-order
correlations (i.e., products of four operators) are supported.

E.g., `ascorr(adag(:n)*a(:n))` returns ``⟨a^\\dagger_n a_n⟩_c + ⟨a^\\dagger_n⟩
⟨a_n⟩`` (which is equal to ``⟨a^\\dagger_n a_n⟩``), while
`ascorr(adag(:n)*a(:m)*a(:n))` returns ``\\langle a_{n}^\\dagger a_{m} a_{n}
\\rangle_{c} + \\langle a_{n}^\\dagger \\rangle \\langle a_{m} \\rangle \\langle
a_{n} \\rangle + \\langle a_{n}^\\dagger \\rangle \\langle a_{m} a_{n}
\\rangle_{c} + \\langle a_{m} \\rangle \\langle a_{n}^\\dagger a_{n}
\\rangle_{c} + \\langle a_{n} \\rangle \\langle a_{n}^\\dagger a_{m}
\\rangle_{c}``.

See also: [`ExpVal`](@ref), [`Corr`](@ref)"""
function ascorr end

ascorr(A::Scalar) = A
for op in (adag,a,f,fdag,σ,σminus,σplus)
    @eval ascorr(A::$op) = ExpVal(A)
end
ascorr(A::OpSum) = ascorr(A.A) + ascorr(A.B)
ascorr(A::OpSumAnalytic) = begin
    # first calculate the correlation for the term in the sum with the "bare" indices, which means that the sum index
    # is assumed to be distinct from the indices of the expressions
    # then, calculate the correction expression for the term where the sum index is the same as any of the non-summed indices,
    # and add the correction between that and the "incorrect" term where index identity is not taken into account
    tmp = ascorr(A.A)
    res = OpSumAnalytic(A.ind,tmp)
    for ii in setdiff(indexset(A),sumindexset(A))
        res += ascorr(replace_index(A.A,A.ind,ii)) - replace_index(tmp,A.ind,ii)
    end
    res
end
function ascorr(A::OpProd)::Operator
    A.A isa OpSumAnalytic && error("should not occurr!")
    A.A isa scal && A.B isa OpSumAnalytic && return A.A*ascorr(A.B)
    preftup, exptup, optup = prodtuples(A)
    pref = prod((scal(1),preftup...,exptup...))
    if length(optup)==0
        pref
    elseif length(optup)==1
        pref*ExpVal(optup[1])
    elseif length(optup)==2
        A,B = optup
        pref*(Corr(A*B) + ExpVal(A)*ExpVal(B))
    elseif length(optup)==3
        A,B,C = optup
        pref*(Corr(A*B*C) + ExpVal(A)*Corr(B*C) + ExpVal(B)*Corr(A*C) + ExpVal(C)*Corr(A*B) + ExpVal(A)*ExpVal(B)*ExpVal(C))
    elseif length(optup)==4
        A,B,C,D = optup
        pref*(Corr(A*B*C*D) + ExpVal(A)*ExpVal(B)*ExpVal(C)*ExpVal(D)
            + ExpVal(A)*Corr(B*C*D) + ExpVal(B)*Corr(A*C*D) + ExpVal(C)*Corr(A*B*D) + ExpVal(D)*Corr(A*B*C)
            + ExpVal(A)*ExpVal(B)*Corr(C*D) + Corr(A*B)*ExpVal(C)*ExpVal(D) + Corr(A*B)*Corr(C*D)
            + ExpVal(A)*ExpVal(C)*Corr(B*D) + Corr(A*C)*ExpVal(B)*ExpVal(D) + Corr(A*C)*Corr(B*D)
            + ExpVal(A)*ExpVal(D)*Corr(B*C) + Corr(A*D)*ExpVal(B)*ExpVal(C) + Corr(A*D)*Corr(B*C))
    else
        error("ERROR: Only correlations up to fourth order are implemented for now")
    end
end
CorrOrExp(A::Operator) = length(A)==1 ? ExpVal(A) : Corr(A)

Avac(A::a) = scal(0)
Avac(A::adag) = A
Avac(A::f) = scal(0)
Avac(A::fdag) = A
Avac(A::σminus) = scal(0)
Avac(A::σplus) = A
# vacuum is an eigenstate of σz
Avac(A::σ) = (A.a == z) ? scal(-1) : A
Avac(A::OpSum) = Avac(A.A) + Avac(A.B)
Avac(A::OpSumAnalytic) = OpSumAnalytic(A.ind,Avac(A.A))
Avac(A::OpProd) = begin
    Bv = Avac(A.B)
    Bv != A.B && return Avac(A.A*Bv)
    # if A.A and A.B commute, try what happens with Avac(A.B*A.A)
    if comm(A.A,A.B) == scal(0)
        Av = Avac(A.A)
        Av != A.A && return Avac(A.B*Av)
    end
    # did not find a way to simplify
    A
end
Avac(A::Scalar) = A

vacA(A::a) = A
vacA(A::adag) = scal(0)
vacA(A::f) = A
vacA(A::fdag) = scal(0)
vacA(A::σminus) = A
vacA(A::σplus) = scal(0)
# vacuum is an eigenstate of σz
vacA(A::σ) = (A.a == z) ? scal(-1) : A
vacA(A::OpSum) = vacA(A.A) + vacA(A.B)
vacA(A::OpSumAnalytic) = OpSumAnalytic(A.ind,vacA(A.A))
vacA(A::OpProd) = begin
    if A.A isa Scalar
        A.A*vacA(A.B)
    else
        # for applying ⟨0|AB, see if ⟨0|A changes A, and if so,
        # keep going with the new operator
        vA = vacA(A.A)
        vA==A.A ?  A : vacA(vA * A.B)
    end
end
vacA(A::Scalar) = A

"""
    Avac(A::Operator), vacA(A::Operator)

Simplify operator by assuming it is applied to the vacuum from the left or
right, respectively. To be precise, `Avac(A)` returns ``A'`` such that ``A'|0⟩ =
A|0⟩``, while `vacA(A)` does the same for ``⟨0|A``."""
Avac, vacA

"""
    vacExpVal(A::Operator,S::Operator=scal(1))

Calculate the vacuum expectation value ``⟨0|S^\\dagger A S|0⟩``, i.e., the
expectation value ``⟨ψ|A|ψ⟩`` for the state defined by ``|ψ⟩= S|0⟩```.
"""
function vacExpVal(A::Operator,stateop::Operator=scal(1))
    # simplify down as much as possible by applying vacuum from left and right
    vsAsv = vacA(Avac(stateop' * A * stateop))
    # only operators that should survive here as operators are σs
    _vacExpVal(vsAsv)
end
_vacExpVal(A::Scalar) = A
# the terms that survive until here have at most a single σ of any particle
# so it does not matter if we are inside a product
_vacExpVal(A::σ) = A.a == z ? scal(-1) : scal(0)
_vacExpVal(A::OpSum) = _vacExpVal(A.A) + _vacExpVal(A.B)
_vacExpVal(A::OpSumAnalytic) = OpSumAnalytic(A.ind,_vacExpVal(A.A))
# we know that the operators here commute (all a and a† have disappeared, and at most a single σ remaining for each particle)
_vacExpVal(A::OpProd) = _vacExpVal(A.A) * _vacExpVal(A.B)

function parse_paramstr(s)
    s = split(s,"_")
    par = Symbol(s[1])
    @assert length(s) <= 2
    inds = length(s)==1 ? () : Meta.parse.(split(s[2],","))
    par,inds
end

macro Pc_str(s)
    par, inds = parse_paramstr(s)
    param(par,'n',inds...)
end
macro Pr_str(s)
    par, inds = parse_paramstr(s)
    param(par,'r',inds...)
end

include("precompile.jl")

end # module


text:
module GraphsExt
  using Graphs

  export read_edges_from_file,save_graph_to_file,
         create_subgraph,all_neighbors,all_degree,occurrence_of_degree

  include("I_O_graphs.jl")
  include("graphs_ext.jl")
end

max_stars_repo_name: johnybx/GraphsExt

max_issues_repo_name: johnybx/GraphsExt

max_forks_repo_name: johnybx/GraphsExt


text:
struct ogUNIFACParam <: EoSParam
    A::PairParam{Float64}
    R::SingleParam{Float64}
    Q::SingleParam{Float64}
end

abstract type ogUNIFACModel <: UNIFACModel end

struct ogUNIFAC{c<:EoSModel} <: ogUNIFACModel
    components::Array{String,1}
    icomponents::UnitRange{Int}
    groups::GroupParam
    params::ogUNIFACParam
    puremodel::Vector{c}
    absolutetolerance::Float64
    references::Array{String,1}
end

@registermodel ogUNIFAC
export ogUNIFAC

function ogUNIFAC(components; puremodel=PR,
    userlocations=String[],
     verbose=false)
    groups = GroupParam(components, ["Activity/UNIFAC/ogUNIFAC/ogUNIFAC_groups.csv"]; verbose=verbose)

    params = getparams(groups, ["Activity/UNIFAC/ogUNIFAC/ogUNIFAC_like.csv", "Activity/UNIFAC/ogUNIFAC/ogUNIFAC_unlike.csv"]; userlocations=userlocations, asymmetricparams=["A"], ignore_missing_singleparams=["A"], verbose=verbose)
    A  = params["A"]
    R  = params["R"]
    Q  = params["Q"]
    icomponents = 1:length(components)

    init_puremodel = [puremodel([groups.components[i]]) for i in icomponents]
    packagedparams = ogUNIFACParam(A,R,Q)
    references = String[]
    model = ogUNIFAC(components,icomponents,groups,packagedparams,init_puremodel,1e-12,references)
    return model
end

function lnγ_comb(model::ogUNIFACModel,V,T,z)
    Q = model.params.Q.values
    R = model.params.R.values

    v  = model.groups.n_flattenedgroups

    x = z ./ sum(z)

    r =[sum(v[i][k]*R[k] for k in @groups) for i in @comps]
    q =[sum(v[i][k]*Q[k] for k in @groups) for i in @comps]

    Φ = r/sum(x[i]*r[i] for i ∈ @comps)
    θ = q/sum(x[i]*q[i] for i ∈ @comps)
    lnγ_comb = @. log(Φ)+(1-Φ)-5*q*(log(Φ/θ)+(1-Φ/θ))
    return lnγ_comb
end

function Ψ(model::ogUNIFACModel,V,T,z)
    A = model.params.A.values
    return @. exp(-A/T)
end


text:
################################################################
# Base fixed datatypes
################################################################
const FILE         = Cvoid # File descriptor
if P4EST_ENABLE_MPI
    const MPI_Comm = MPI.MPI_Comm
else
    const MPI_Comm = Cint
end
const MPI_Datatype = MPI.MPI_Datatype
const MPI_File     = MPI.MPI_File
const sc_MPI_Comm  = MPI_Comm
const ssize_t      = Cssize_t

################################################################
# Non-wrapped but required data types
# p4est_base.h header not included as it depends on p4est_config.h and sc_config.h.
# This files are automatically generated during configuration.
################################################################
const p4est_qcoord_t = Cint         # int32_t
const p4est_locidx_t = Cint         # int32_t
const p4est_topidx_t = Cint         # int32_t
const p4est_gloidx_t = Clong        # int64_t

################################################################
# fix sc_array struct.
# Wrong element_size and elem_cont datatypes (Cint -> Cssize_t)
################################################################
struct sc_array
    elem_size::Csize_t
    elem_count::Csize_t
    byte_alloc::ssize_t
    array::Ptr{Cchar}
end

const sc_array_t = sc_array

################################################################
# C-Unions: Non supported
#
# A union is a special data type available in C.
# C-Unions allows to store different data types in the same memory location.
# The memory occupied by a C-Union will be equal to the largest member of the union.
#
# From the Julia point of view, we treat C-Unions as fixed size plain memory regions.
# These memory regions can be later reinterpreted as any of the original members of the C-Union.
################################################################

function _pointer_to_c_memory(obj::S) where {S}
    Base.unsafe_convert(Ptr{Cvoid}, [obj])
end

function _unsafe_reinterpret(::Type{T}, obj::S) where {T,S}
# Use unsafe_wrap to avoid copy from unsafe_load
# unsafe_load(Ptr{T}(Base.unsafe_convert(Ptr{Cvoid}, Ref{S}(obj))))
    unsafe_wrap(Array, Ptr{T}(_pointer_to_c_memory(obj)), 1)[]
end

    ############################################
    # Primitive type: quadrant_data
    # represents the following p4est C-Unions:
    #   - p4est_quadrant_data
    #   - p6est_quadrant_data
    #   - p8est_quadrant_data
    ############################################
struct piggy1 which_tree::p4est_topidx_t; owner_rank::Cint; end
struct piggy2 which_tree::p4est_topidx_t; from_tree::p4est_topidx_t; end
struct piggy3 which_tree::p4est_topidx_t; local_num::p4est_locidx_t; end

quadrant_data_union = Union{Ptr{Cvoid}, Clong, Cint, p4est_topidx_t, piggy1, piggy2, piggy3}
quadrant_data_properties = (:user_data, :user_long, :user_int, :which_tree, :piggy1, :piggy2, :piggy3)
primitive type quadrant_data 8*sizeof(quadrant_data_union) end # Biggest data size in union (8 bytes - 64 bits)

function Base.getproperty(x::quadrant_data, name::Symbol)
    if name == :user_data
        _unsafe_reinterpret(Ptr{Cvoid}, x)
    elseif name == :user_long
        _unsafe_reinterpret(Clong, x)
    elseif name == :user_int || name == :which_tree
        _unsafe_reinterpret(Cint, x)
    elseif name == :piggy1
        _unsafe_reinterpret(piggy1, x)
    elseif name == :piggy2
        _unsafe_reinterpret(piggy2, x)
    elseif name == :piggy3
        _unsafe_reinterpret(piggy3, x)
    else
        getfield(x, name)
    end
end

Base.propertynames(x::quadrant_data) = (quadrant_data_properties..., fieldnames(DataType)...)

const p4est_quadrant_data = quadrant_data
const p6est_quadrant_data = quadrant_data
const p8est_quadrant_data = quadrant_data

    ############################################
    # Primitive type: iter_side_data_2
    # represents the following p4est C-Unions:
    #   - p4est_iter_face_side_data
    #   - p8est_iter_edge_side_data
    ############################################

struct full is_ghost::Int8; quad::Ptr{Cvoid}; quadid::p4est_locidx_t; end
struct hanging{S} is_ghost::NTuple{S,Int8}; quad::NTuple{S,Ptr{Cvoid}}; quadid::NTuple{S,p4est_locidx_t}; end
iter_side_data_properties = (:full, :hanging)

# Biggest data size in union (24 bytes - 192 bits)
iter_side_data_2_union = Union{full, hanging{2}}
primitive type iter_side_data_2 8*sizeof(iter_side_data_2_union) end

function Base.getproperty(x::iter_side_data_2, name::Symbol)
    if name == :full
        _unsafe_reinterpret(full, x)
    elseif name == :hanging
        _unsafe_reinterpret(hanging{2}, x)
    else
        getfield(x, name)
    end
end

Base.propertynames(x::iter_side_data_2) = (iter_side_data_properties..., fieldnames(DataType)...)

const p4est_iter_face_side_data = iter_side_data_2
const p8est_iter_edge_side_data = iter_side_data_2

    ############################################
    # Primitive type: iter_side_data_4
    # represents the following p4est C-Unions:
    #   - p8est_iter_face_side_data
    ############################################

# Biggest data size in union (32 bytes - 256 bits)
iter_side_data_4_union = Union{full, hanging{4}}
primitive type iter_side_data_4 8*sizeof(iter_side_data_4_union) end

function Base.getproperty(x::iter_side_data_4, name::Symbol)
    if name == :full
        _unsafe_reinterpret(full, x)
    elseif name == :hanging
        _unsafe_reinterpret(hanging{4}, x)
    else
        getfield(x, name)
    end
end

Base.propertynames(x::iter_side_data_4) = (iter_side_data_properties..., fieldnames(DataType)...)

const p8est_iter_face_side_data = iter_side_data_4


text:
using Cytof5
using RCall, Distributions, Random
using JLD2, FileIO

R"""
library(xtable)
library(FlowSOM)
library(flowCore)
library(rcommon)
library(cytof3)
source("est_Z_from_clusters.R")
set.seed(3)
"""

include("../sim_study/util.jl")

if length(ARGS) == 0
  SIMDIR = "../sim_study/results/I3_J32_N_factor100_K8_L4_K_MCMC10_L_MCMC5_b0PriorSd0.1_b1PriorScale0.1_SEED0/"
  RESULTS_DIR = "results/"
else
  # Where to get data
  SIMDIR = ARGS[1]
  # Where to put results from FlowSOM analysis
  RESULTS_DIR = ARGS[2]
end

# Where in results dir to put output
expName(dir) = filter(x -> length(x) > 0, split(dir, "/"))[end]
OUTPUT_DIR = "$(RESULTS_DIR)/$(expName(SIMDIR))/"
mkpath(OUTPUT_DIR)

# READ DATA
println("Loading Data ...")
@load "$(SIMDIR)/output.jld2" out dat ll lastState c y_dat metrics

# MAIN
function replaceMissing(yi, x)
  out = deepcopy(yi)
  out[isnan.(out)] .= x
  return out
end

# Preimpute missing values
y_tilde = [ replaceMissing(yi, -20) .+ 0 for yi in y_dat.y ]

# Convert data to flowsom formats
Y_tilde = R"Reduce(rbind, $y_tilde)"
@rput Y_tilde
ff_Y = R"colnames(Y_tilde) <- 1:NCOL(Y_tilde); flowFrame(Y_tilde)"

# FlowSOM
println("Running FlowSOM...")
J = y_dat.J
I = y_dat.I
N = y_dat.N
FlowSOM = R"FlowSOM::FlowSOM"
@time fSOM = FlowSOM(ff_Y,
                     # Input options:
                     colsToUse = 1:J,
                     # Metaclustering options:
                     #nClus = 20,
                     maxMeta=20,
                     # Seed for reproducible results:
                     seed = 42)

@rput fSOM I J N
R"""
idx_upper = cumsum(N)
idx_lower = c(1,idx_upper[-I]+1)
idx = cbind(idx_lower, idx_upper)

fSOMClus = fSOM$meta[fSOM$FlowSOM$map$mapping[,1]]
fsEst = as.list(1:$I)
fsClus = as.numeric(fSOMClus)

mult=1
png($OUTPUT_DIR %+% 'YZ%03d_FlowSOM.png', height=500*mult, width=500*mult)
for (i in 1:$I) {
  clus = fsClus[idx[i,1]:idx[i,2]]
  print(length(unique(clus))) # Number of clusters learned
  clus = relabel_clusters(clus)
  my.image($(y_dat.y)[[i]][order(clus),], col=blueToRed(11), zlim=c(-5,5), addL=TRUE,
           na.color='black', cex.y.leg=2, xlab='cell types',  ylab='cells',
           cex.lab=1.5, cex.axis=1.5, xaxt='n',
           f=function(z) {
             add.cut(clus)
             axis(1, at=1:$(J), fg='grey', las=2, cex.axis=1.5)
           })
}
dev.off()

println("Computing ARI ...")
true.clus.ls <- $(dat[:lam])
fs.clus.ls = lapply(1:I, function(i) fsClus[idx_lower[i]:idx_upper[i]])

ARI = sapply(1:I, function(i) {
  ari(true.clus.ls[[i]], fs.clus.ls[[i]])
})

ARI_all = 'FlowSOM'=ari(unlist(true.clus.ls), fsClus)

sink($OUTPUT_DIR %+% "ari.txt")
  println("ARI:")
  print(ARI)
  println("ARI all:")
  print(ARI_all)
sink()
"""


text:
function get_h(obj::HessenbergTransformer)
    return jcall(obj, "getH", RealMatrix, ())
end

function get_p(obj::HessenbergTransformer)
    return jcall(obj, "getP", RealMatrix, ())
end

function get_pt(obj::HessenbergTransformer)
    return jcall(obj, "getPT", RealMatrix, ())
end


text:
# An error thrown by libhdf5
mutable struct H5Error <: Exception
    msg::String
    id::hid_t
end

macro h5error(msg)
    # Check if the is actually any errors on the stack. This is necessary as there are a
    # small number of functions which return `0` in case of an error, but `0` is also a
    # valid return value, e.g. `h5t_get_member_offset`

    # This needs to be a macro as we need to call `h5e_get_current_stack()` _before_
    # evaluating the message expression, as some message expressions can call API
    # functions, which would clear the error stack.
    quote
        err_id = h5e_get_current_stack()
        if h5e_get_num(err_id) > 0
            throw(H5Error($(esc(msg)), err_id))
        else
            h5e_close_stack(err_id)
        end
    end
end

Base.cconvert(::Type{hid_t}, err::H5Error) = err
Base.unsafe_convert(::Type{hid_t}, err::H5Error) = err.id

function Base.close(err::H5Error)
    if err.id != -1 && isvalid(err)
        h5e_close_stack(err)
        err.id = -1
    end
    return nothing
end
Base.isvalid(err::H5Error) = err.id != -1 && h5i_is_valid(err)

Base.length(err::H5Error) = h5e_get_num(err)
Base.isempty(err::H5Error) = length(err) == 0

function H5Error(msg::AbstractString)
    id = h5e_get_current_stack()
    err = H5Error(msg, id)
    finalizer(close, err)
    return err
end

const SHORT_ERROR = Ref(true)

function Base.showerror(io::IO, err::H5Error)
    n_total = length(err)
    print(io, "$(typeof(err)): ", err.msg)
    print(io, "\nlibhdf5 Stacktrace:")
    h5e_walk(err, H5E_WALK_UPWARD) do n, errptr
        n += 1
        if SHORT_ERROR[] && 1 < n < n_total
            return nothing
        end
        errval = unsafe_load(errptr)
        print(io, "\n", lpad("[$n] ", 4 + ndigits(n_total)))
        if errval.func_name != C_NULL
            printstyled(io, unsafe_string(errval.func_name); bold=true)
            print(io, ": ")
        end
        major = h5e_get_msg(errval.maj_num)[2]
        minor = h5e_get_msg(errval.min_num)[2]
        print(io, major, "/", minor)
        if errval.desc != C_NULL
            printstyled(io, "\n", " "^(4 + ndigits(n_total)), unsafe_string(errval.desc), color=:light_black)
        end
        if SHORT_ERROR[] && n == 1 && n_total > 2
            print(io, "\n", lpad("⋮", 2 + ndigits(n_total)))
        end
    end
    return nothing
end


text:
"""
  simulate_one_predictor(iter, m, p)

  It asks for number of iterations `iter` to be run, number of observations `m` in treatment
  and control groups `p` for the simulation of Semiparametric AUC regression adjusting for one discrete
  covariate. In this simulation, true model parameters are as follows: β0 = 0.15, β1 = 0.50, β2 = 1.
"""
simulate_one_predictor(iter = 500, m = 100, p = 120)

function simulate_one_predictor(;iter = 500, m = 100, p = 120)
    iter = iter
    finvhat = gamma1 = Array(Float64, 3)
    AUChat = Array(Float64, 3)
    Vhat_auchat = Array(Float64, 3)
    lo = up = Array(Float64, iter, 3)
    d  = Array(Float64, m, 3)
    nd = Array(Float64, p, 3)

    m_betas = Array(Float64, iter, 3)
    sd_betas = Array(Float64, iter, 3)
    lower = upper = score =  cov_b = Array(Float64, iter, 3)
    v_finv_auchat = gamma = Array(Float64, 3)
    all_var = Array(Float64, 3)
    var_finv_auchat = Array(Float64, iter, 3)

    for z in range(1,iter)
      d1 = [0.0, 0.0, 0.0]
      d2 = [0, 0.50, 1.00]
      d0 = 0.15
      y = 1:3
      for k in y
        result = Array(Float64, k, 3)
        u1 = randexp(p) # rand(Exponential(1), p)
        u2 = randexp(m) # rand(Exponential(1), m)
        d[:,k]=-log(u2) + d0 +(d1[k] + d2[k])
        nd[:,k]=-log(u1) + d1[k]

        result= calculate_auc_simulation(ya = d[:,k], yb = nd[:,k])
        AUChat[k]=result[1]
        finvhat[k]=result[2]
        Vhat_auchat[k] = result[3]
        v_finv_auchat[k] = Vhat_auchat[k]/(((AUChat[k])^2)*(1-(AUChat[k]))^2)  #Variance of F inverse
      end
      gamma1 = finvhat
      Z = reshape([1,0,0,1,1,0,1,0,1], 3, 3)'

      tau  =  diagm([1/i for i in v_finv_auchat])

      ztauz = inv(Z' * tau * Z)
      var_betas = diag(ztauz)
      std_error = sqrt(var_betas)
      betas = ztauz * Z' * tau * gamma1

      m_betas[z,:]  =  betas
      var_finv_auchat[z,:] = var_betas
    end
    lo = m_betas .- 1.96*sqrt(var_finv_auchat)
    up = m_betas .+ 1.96*sqrt(var_finv_auchat)
    ci = hcat(lo,up)
    ci_betas = ci[:,[1,4,2,5,3,6]]
    return(m_betas, var_finv_auchat, ci_betas, iter)
end


text:
export inf_sites
#=
A simplified version of the infinite sites model where sites evolve independently.
Called by nearly-neutral/src/run_infsites.jl
Author:  Alden H. Wright, Department of Computer Science, Univeristy of Montana, Missoula, MT 59812 USA
    alden.wright@umontana.edu
Example run:
[src]$ julia run.jl examples/is_example1
=#

@doc """ function inf_sites( )
Do a infinite sites simulation with popsize N, L loci, mutation rate mu per locus, ngens generations.
The site update is done in function update_sites() in src/site_collection.jl.
Returns a infs_result_type object, but also prints relevant information.
"""
function inf_sites( N::Int64, N_mu::Float64, ngens::Int64; dfe::Function=dfe_neutral)
  fitness_table = Dict{Int64,Float64}()
  # Create a new empty site collection
  sc = site_collection( N, 1.0 )   # fix_minimum = 1.0 so fixes only if count is N
  g_limit = 100000   # Maximum of the number of extra generations to run while waiting for all sites to fix
  mu = N_mu/N
  new_id = 1
  i = 1
  N_sum = 0
  g_count = 0
  g = 1
  num_gens = 0
  done = false
  while !done && g < g_limit+ngens
    num_mutations = rand(Binomial(N,mu))
    if g <= ngens
      for j in 1:num_mutations
        fit = dfe_fitness( new_id, dfe, fitness_table )  # Set fitness of new_id
        sc_push!(sc,site(new_id,N,g,fit))
        #println("new site: id: ",new_id,"  start gen: ",sc.list[new_id].start_gen,"  sum_h: ",sc.list[new_id].sum_heteroz)
        new_id += 1
      end
      g_count += 1
    end
    update_sites!( sc, g, N )
    g += 1
    done = (g > ngens) && length(sc.active) == 0
  end
  #print_summary( sc )   # print a summary of the values stored in site_collection sc
  count_adv_del_fixed( sc )   # sets sc.count_fixed_adv and sc.count_fixed_del
  check_globals_against_site_collection( sc )    # Check that global statistics agree with those in sc.
  sc
end

@doc """ function dfe_fitness( p::Int64, dfe::Function, fitness_table::Dict{Int64,Float64} )
Fitness function as saved in dictionary fitness_table.
"""
function dfe_fitness( p::Int64, dfe::Function, fitness_table::Dict{Int64,Float64} )
  val = get(fitness_table,p,-1.0)
  if val == -1.0   # p is not in fitness table
    val = dfe( p )
    fitness_table[p] = val
  end
  val
end

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
### A segment of a Generalized Metropolis-Hastings run, to be executed on a single worker
immutable GMHSegment <: AbstractJobSegment
    model::AbstractModel
    samplerstate::AbstractSamplerState
end

### Factory function
function _segment(p::MHRuntimePolicy,m::AbstractModel,s::AbstractSampler,nproposals::Integer)
    @assert traitvalue(p.runner) == :generalized
    state = samplerstate(s,nproposals,p.sampletype,p.calculationtype,true)
    GMHSegment(m,state)
end

### Perform one iteration
function iterate!(seg::GMHSegment,indicator::AbstractSamplerState)
    prepareauxiliary!(indicator,seg.samplerstate)
    propose!(seg.samplerstate)
    geometry!(seg.model,proposals(seg.samplerstate))
    acceptance!(seg.samplerstate)
end

### Update the indicatorstate with the sample of this auxiliary state
prepare!(seg::GMHSegment,indicator::AbstractSamplerState,i::Int) = prepareindicator!(indicator,seg.samplerstate,i)

### Get the samples
getsamples(seg::GMHSegment,sampleindex) = copy(proposals(seg.samplerstate),sampleindex)

### Tune the samplerstate
tune!(seg::GMHSegment,tvals...) = tune!(seg.samplerstate,tvals...)

###
function show(io::IO,s::GMHSegment)
    println(io,"GMHSegment with: ")
    println(io,"model:")
    show(io,s.model)
    println(io)
    println(io,"samplerstate:")
    show(io,s.samplerstate)
    println(io)
    nothing
end


text:
using Test
using LinearAlgebra
using PDMats
using Distributions
using DistributionsAD
using ConditionalDists
using Flux

using DistributionsAD: TuringMvNormal
using ConditionalDists: BatchMvNormal, SplitLayer

include("cond_dist.jl")
include("cond_mvnormal.jl")
include("categorical.jl")
include("utils.jl")

# for the BatchMvNormal tests to work BatchMvNormals have to be functors!
include("batch_mvnormal.jl")


text:
function FieldKeplerianPropagator(arg0::FieldOrbit)
    return FieldKeplerianPropagator((FieldOrbit,), arg0)
end

function FieldKeplerianPropagator(arg0::FieldOrbit, arg1::AttitudeProvider)
    return FieldKeplerianPropagator((FieldOrbit, AttitudeProvider), arg0, arg1)
end

function FieldKeplerianPropagator(arg0::FieldOrbit, arg1::AttitudeProvider, arg2::RealFieldElement)
    return FieldKeplerianPropagator((FieldOrbit, AttitudeProvider, RealFieldElement), arg0, arg1, arg2)
end

function FieldKeplerianPropagator(arg0::FieldOrbit, arg1::AttitudeProvider, arg2::RealFieldElement, arg3::RealFieldElement)
    return FieldKeplerianPropagator((FieldOrbit, AttitudeProvider, RealFieldElement, RealFieldElement), arg0, arg1, arg2, arg3)
end

function FieldKeplerianPropagator(arg0::FieldOrbit, arg1::RealFieldElement)
    return FieldKeplerianPropagator((FieldOrbit, RealFieldElement), arg0, arg1)
end

function reset_initial_state(obj::FieldKeplerianPropagator, arg0::FieldSpacecraftState)
    return jcall(obj, "resetInitialState", void, (FieldSpacecraftState,), arg0)
end


text:
###############################################################################
#
# Download the artifacts manually to avoid slow downloading from julia server
#
###############################################################################
"""
Artifacts stored on FTP cannot be installed directly for Windows using Julia 1.5. This problem may be fixed in future Julia versions. To address the
    problem for lower Julia versions, `PkgUtility` provides function `predownload_artifact!` to download and unpack the artifact manually for Windows:

$(METHODLIST)

"""
function predownload_artifact! end

"""
What `predownload_artifact!` does are
- determine if the artifact exists already
- if `true`, skip the artifact installation
- if `false`
  - install the artifact using the `download_artifact` from `Pkg.Artifacts`
  - if the installation fails (e.g., for Windows)
    - download the `.tar.gz` manually
    - unpack the  `.tar.gz` manually

Method for this function is

    predownload_artifact!(art_name::String, artifact_toml::String)

Download the artifact from given server if it does not exist, given
- `art_name` Artifact name to install
- `artifact_toml` `Artifacts.toml` file location

---
Examples
```julia
# download the artifact manually if default methods fail
predownload_artifact!("test_artifact", "Artifacts.toml");
```
"""
predownload_artifact!(art_name::String, artifact_toml::String) =
(
    _meta = artifact_meta(art_name, artifact_toml);
    _hash = Base.SHA1(_meta["git-tree-sha1"]);

    # if artifact exists already, skip it
    if artifact_exists(_hash)
        return nothing;
    end;

    # try to download the artifact from all entries if it does not exist
    for _entry in _meta["download"]
        _url          = _entry["url"];
        _tarball_hash = _entry["sha256"];

        @info tinfo("Downloading artifact from $(_url)...");
        _succeeded    = download_artifact(_hash, _url, _tarball_hash);
        if _succeeded
            break;
        end;
    end;

    # if artifact is installed successfully, skip it
    if artifact_exists(_hash)
        return nothing;
    end;

    # download and unpack the artifact manually if above fails
    @info tinfo("Install the artifact manually...");
    for _entry in _meta["download"]
        _dest_dir  = artifact_path(_hash; honor_overrides=false);
        _url       = _entry["url"];

        @info tinfo("Downloading artifact manually from $(_url)...");
        _succeeded = download(_url, "temp_artifact.tar.gz");
        if _succeeded == "temp_artifact.tar.gz"
            unpack("temp_artifact.tar.gz", _dest_dir);
            rm("temp_artifact.tar.gz");
            break;
        end;
    end;

    return nothing;
)


text:
using Revise
using AdFem
using PyCall
using LinearAlgebra
using PyPlot
using SparseArrays
using MAT
np = pyimport("numpy")

# mode = "data" generate data
# mode != "data" train
mode = "train"

β = 1/4; γ = 1/2
a = b = 0.1
m = 40
n = 20
h = 0.01
NT = 50
Δt = 5/NT
bdedge = []
for j = 1:n
  push!(bdedge, [(j-1)*(m+1)+m+1 j*(m+1)+m+1])
end
bdedge = vcat(bdedge...)

bdnode = Int64[]
for j = 1:n+1
  push!(bdnode, (j-1)*(m+1)+1)
end

# λ = Variable(1.0)
# μ = Variable(1.0)
# invη = Variable(1.0)

λ = constant(2.0)
μ = constant(0.5)
if mode=="data"
  global invη = placeholder(12.0)
else
  global invη = Variable(1.0)
end

G = tensor([1/Δt+μ*invη -μ/3*invη 0.0
  -μ/3*invη 1/Δt+μ*invη-μ/3*invη 0.0
  0.0 0.0 1/Δt+μ*invη])
S = tensor([2μ/Δt+λ/Δt λ/Δt 0.0
    λ/Δt 2μ/Δt+λ/Δt 0.0
    0.0 0.0 μ/Δt])

invG = inv(G)
# invG = Variable([1.0 0.0 0.0
#                 0.0 1.0 0.0
#                 0.0 0.0 1.0])
H = invG*S

M = compute_fem_mass_matrix1(m, n, h)
Zero = spzeros((m+1)*(n+1), (m+1)*(n+1))
M = SparseTensor([M Zero;Zero M])

K = compute_fem_stiffness_matrix(H, m, n, h)
C = a*M + b*K # damping matrix

L = M + γ*Δt*C + β*Δt^2*K
L, Lbd = fem_impose_Dirichlet_boundary_condition(L, bdnode, m, n, h)
# error()

a = TensorArray(NT+1); a = write(a, 1, zeros(2(m+1)*(n+1))|>constant)
v = TensorArray(NT+1); v = write(v, 1, zeros(2(m+1)*(n+1))|>constant)
d = TensorArray(NT+1); d = write(d, 1, zeros(2(m+1)*(n+1))|>constant)
U = TensorArray(NT+1); U = write(U, 1, zeros(2(m+1)*(n+1))|>constant)
Sigma = TensorArray(NT+1); Sigma = write(Sigma, 1, zeros(4*m*n, 3)|>constant)
Varepsilon = TensorArray(NT+1); Varepsilon = write(Varepsilon, 1,zeros(4*m*n, 3)|>constant)

Forces = zeros(NT, 2(m+1)*(n+1))
for i = 1:NT
  T = eval_f_on_boundary_edge((x,y)->0.1, bdedge, m, n, h)
  T = [T zeros(length(T))]
  rhs = compute_fem_traction_term(T, bdedge, m, n, h)

  if i*Δt>3.0
    rhs = zero(rhs)
  end
  Forces[i, :] = rhs
end
Forces = constant(Forces)

function condition(i, tas...)
  i <= NT
end

function body(i, tas...)
  a_, v_, d_, U_, Sigma_, Varepsilon_ = tas
  a = read(a_, i)
  v = read(v_, i)
  d = read(d_, i)
  U = read(U_, i)
  Sigma = read(Sigma_, i)
  Varepsilon = read(Varepsilon_, i)

  F = compute_strain_energy_term(Sigma*invG/Δt, m, n, h) - K * U
  rhs = Forces[i] - F
  td = d + Δt*v + Δt^2/2*(1-2β)*a
  tv = v + (1-γ)*Δt*a
  rhs = rhs - C*tv - K*td

  # rhs[[bdnode; bdnode.+(m+1)*(n+1)]] .= 0.0

  rhs = scatter_update(rhs, constant([bdnode; bdnode.+(m+1)*(n+1)]), constant(zeros(2*length(bdnode))))

  a = L\rhs # bottleneck
  d = td + β*Δt^2*a
  v = tv + γ*Δt*a
  U_new = d
  Varepsilon_new = eval_strain_on_gauss_pts(U_new, m, n, h)
  Sigma_new = Sigma*invG/Δt +  (Varepsilon_new-Varepsilon)*(invG*S)

  i+1, write(a_, i+1, a), write(v_, i+1, v), write(d_, i+1, d), write(U_, i+1, U_new),
        write(Sigma_, i+1, Sigma_new), write(Varepsilon_, i+1, Varepsilon_new)
end

i = constant(1, dtype=Int32)
# error()
_, _, _, _, u, sigma, varepsilon = while_loop(condition, body,
                  [i, a, v, d, U, Sigma, Varepsilon])

U = stack(u)
Sigma = stack(sigma)
Varepsilon = stack(varepsilon)

if mode!="data"
  global Uval = matread("U.mat")["U"]
  U.set_shape((NT+1, size(U, 2)))
  idx = 1:m+1
  global loss = sum((U[:, idx] - Uval[:, idx])^2)
end

# opt = AdamOptimizer(1.0).minimize(loss)
sess = Session(); init(sess)

if mode=="data"
  matwrite("U.mat", Dict("U"=>run(sess, U)))
  error()
end

@info run(sess, loss)
v_ = []
i_ = []
l_ = []
cb = (v, i, l)->begin
  println("[$i] loss = $l")
  println(v)
  push!(v_, [x for x in v])
  push!(i_, i)
  push!(l_, l)
end
loss_ = BFGS!(sess, loss, vars=[λ, μ, invη], callback=cb)
# matwrite("R.mat", Dict("V"=>v_, "L"=>l_))

# for i = 1:1000
#   _, l, invη_ = run(sess, [opt, loss, invη])
#   @show i, l, invη_
# end

# ηs = 11:0.1:13
# losses = []
# for η in ηs
#    push!(losses,run(sess, loss, invη=>η))
# end
# plot(ηs, losses)

# Uval, Sigmaval, Varepsilonval = run(sess, [U, Sigma, Varepsilon])
# Uval[idx]

# visualize_displacement(U, m, n, h; name = "_eta$η", xlim_=[-0.01,0.5], ylim_=[-0.05,0.22])
# # visualize_displacement(U, m, n, h;  name = "_viscoelasticity")
# # visualize_stress(H, U, m, n, h;  name = "_viscoelasticity")

# close("all")
# figure(figsize=(15,5))
# subplot(1,3,1)
# idx = div(n,2)*(m+1) + m+1
# plot((0:NT)*Δt, Uval[:, idx])
# xlabel("time")
# ylabel("\$u_x\$")

# subplot(1,3,2)
# idx = 4*(div(n,2)*m + m)
# plot((0:NT)*Δt, Sigmaval[:,idx,1])
# xlabel("time")
# ylabel("\$\\sigma_{xx}\$")

# subplot(1,3,3)
# idx = 4*(div(n,2)*m + m)
# plot((0:NT)*Δt, Varepsilonval[:,idx,1])
# xlabel("time")
# ylabel("\$\\varepsilon_{xx}\$")
# savefig("visco_eta$η.png")