text:
index(e) = e.index
index1(e) = e.index1
index2(e) = e.index2
inner(e) = e.inner
ref(e) = e.ref
ref1(e) = e.ref1
ref2(e) = e.ref2
# ref11(e) = e.ref11
# ref12(e) = e.ref12
# ref21(e) = e.ref21
# ref22(e) = e.ref22
refval(e) = e.ref.x
frefval(e) = e.fref.x
frefval1(e) = e.fref1.x
frefval2(e) = e.fref2.x
refval1(e) = e.ref1.x
refval2(e) = e.ref2.x
refval11(e) = e.ref11.x
refval12(e) = e.ref12.x
refval21(e) = e.ref21.x
refval22(e) = e.ref22.x
setrefval(e,val) = ref(e)[] = val
setrefval1(e,val) = ref1(e)[] = val
setrefval2(e,val) = ref2(e)[] = val
addrefval(e,val) = ref(e)[] += val
# addrefval1(e,val) = ref1(e)[] += val
# addrefval2(e,val) = ref2(e)[] += val

# non-caching variants
@inline non_caching_eval(f::Sink{Field},y,x,p=nothing) = inner(f)(y,x,p)
@inline function non_caching_eval(f::Field1{E,I},y,x,p=nothing) where {E,I}
    non_caching_eval(inner(f),y,x,p)
    @simd for i in eachindex(f.es)
        non_caching_eval(f.es[i],y,x,p)
    end
end
@inline function non_caching_eval(f::Field1{E,Nothing},y,x,p=nothing) where E
    @simd for i in eachindex(f.es)
        non_caching_eval(f.es[i],y,x,p)
    end
end

@inline non_caching_eval(e::Variable,x,p=nothing)  = @inbounds getindex(x,index(e))
@inline non_caching_eval(e::Parameter,x,p=nothing)  = @inbounds getindex(p,index(e))
@inline non_caching_eval(e::Constant,x,p=nothing)  = refval(e)
@inline function non_caching_eval(e::ExpressionSum{E,I},x,p=nothing) where {E,I}
    res = non_caching_eval(e.inner,x,p)
    @simd for i in eachindex(e.es)
        @inbounds res = add_sum(res,non_caching_eval(e.es[i],x,p))
    end
    return res
end
@inline function non_caching_eval(e::ExpressionSum{E,Nothing},x,p=nothing) where E
    res = non_caching_eval(e.es[1],x,p)
    @simd for i in 2:length(e.es)
        @inbounds res = add_sum(res,non_caching_eval(e.es[i],x,p))
    end
    return res
end

@inline non_caching_eval(::GradientNull,z,x,p=nothing,d0=1) = nothing
@inline non_caching_eval(d::Gradient0,y,x,p=nothing,d0=1) = (@inbounds y[d.offset] += d0; nothing)
@inline non_caching_eval(d::Gradient0,(j,y)::Tuple{Int,M},x,p=nothing,d0=1) where M <: AbstractMatrix = (@inbounds y[j,d.offset] += d0; nothing)
@inline non_caching_eval(d::Gradient0,(j,y)::Tuple{Int,M},x,p=nothing,d0=1) where M <: AbstractVector = (@inbounds y[d.offset] += d0; nothing)
@inline function non_caching_eval(d::GradientSum{D,I} where {D,I},y,x,p=nothing,d0=1)
    non_caching_eval(inner(d),y,x,p,d0)
    @simd for i in eachindex(d.ds)
        non_caching_eval(d.ds[i],y,x,p,d0)
    end
end
@inline function non_caching_eval(d::GradientSum{D,Nothing} where D,y,x,p=nothing,d0=1)
    @simd for i in eachindex(d.ds)
        non_caching_eval(d.ds[i],y,x,p,d0)
    end
end
@inline non_caching_eval(e::IndexedExpression{E},y,x,p=nothing) where {E} = (y[index(e)] = non_caching_eval(e.e,x,p))
@inline non_caching_eval(e::JacobianEntry{E},y,x,p=nothing) where {E} = non_caching_eval(e.e,(index(e),y),x,p)
@inline non_caching_eval(e::FieldNull,y,x,p=nothing) = nothing


text:
module Example2018

export Point
struct Point{T}
	x::T
	y::T
end

function Base. +(p::Point,q::Point)
	Point(p.x +q.x, p.y+q.y)
end

export MPoint
mutable struct MPoint{T}
	x::T
	y::T
end

end
# module

ext: jl


text:
using Wavelets

wf = waveletfilter(WT.db4)
N = 1024;
x0 = rand(N,N);
L = int(log2(N)-2)
tn = 10
#y=x0*0.0

println("iwt (N=",N,"x",N,"), ", L, " levels")
f(x0, L, wf) = for i = 1:tn; idwt(x0, wf, L); end
f(x0, L, wf);
@time f(x0, L, wf);


text:
"""
    s_per_degree(s_per_km, r_earth=$R_EARTH_KM_DEFAULT)

Convert horizontal slowness in s/km at the surface of the Earth to s/°.

`r_earth` is the radius of the spherical body in km.
"""
s_per_degree(s_km, r_earth=R_EARTH_KM_DEFAULT) = s_km*2π*r_earth/360

"""
    s_per_km(s_per_degree, r_earth=$R_EARTH_KM_DEFAULT)

Convert horizontal slowness in s/° at the surface of the Earth to s/km.

`r_earth` is the radius of the spherical body in km.
"""
s_per_km(s_degree, r_earth=R_EARTH_KM_DEFAULT) = s_degree*360/(2π*r_earth)

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
_tri(t, period) = (2 / π) * abs(asin(sin(π * (t - 1) / period)))
_sin(t, period) = abs(sin(π * (t - 1) / period))
_cycle(λ0, λ1, g) = abs(λ0 - λ1) * g + min(λ0, λ1)

"""
    Triangle{T, S<:Integer}(range0, range1, period)
    Triangle(;λ0, λ1, period)

A [triangle wave](https://en.wikipedia.org/wiki/Triangle_wave) schedule
with `period`.
The output conforms to
```text
abs(λ0 - λ1) * (2 / π) * abs(asin(sin(π * (t - 1) / period))) + min(λ0, λ1)
```

# Arguments
- `range0`/`λ0`: the first range endpoint
- `range1`/`λ1`: the second range endpoint
- `period::Integer`: the period
"""
struct Triangle{T, S<:Integer}
    range0::T
    range1::T
    period::S
end
Triangle(;λ0, λ1, period) = Triangle(λ0, λ1, period)

(schedule::Triangle)(t) = _cycle(schedule.range0, schedule.range1, _tri(t, schedule.period))

Base.eltype(::Type{<:Triangle{T}}) where T = T
Base.IteratorSize(::Type{<:Triangle}) = Base.IsInfinite()

Base.iterate(schedule::Triangle, t = 1) = schedule(t), t + 1

Base.axes(::Triangle) = (OneToInf(),)

"""
    TriangleDecay2{T, S<:Integer}(range0, range1, period)
    TriangleDecay2(;λ0, λ1, period)

A [triangle wave](https://en.wikipedia.org/wiki/Triangle_wave) schedule
with `period` and half the amplitude each cycle.
The output conforms to
```text
abs(λ0 - λ1) * Triangle(t) / (2^floor((t - 1) / period)) + min(λ0, λ1)
```
where `Triangle(t)` is `(2 / π) * abs(asin(sin(π * (t - 1) / schedule.period)))` (see [`Triangle`](#)).

# Arguments
- `range0`/`λ0`: the first range endpoint
- `range1`/`λ1`: the second range endpoint
- `period::Integer`: the period
"""
struct TriangleDecay2{T, S<:Integer}
    range0::T
    range1::T
    period::S
end
TriangleDecay2(;λ0, λ1, period) = TriangleDecay2(λ0, λ1, period)

(schedule::TriangleDecay2)(t) = _cycle(schedule.range0, schedule.range1,
                                       _tri(t, schedule.period) / (2^fld(t - 1, schedule.period)))

Base.eltype(::Type{<:TriangleDecay2{T}}) where T = T
Base.IteratorSize(::Type{<:TriangleDecay2}) = Base.IsInfinite()

Base.iterate(schedule::TriangleDecay2, t = 1) = schedule(t), t + 1

Base.axes(::TriangleDecay2) = (OneToInf(),)

"""
    TriangleExp{T, S<:Integer}(range0, range1, period, decay)
    TriangleExp(λ0, λ1, period, γ)
    TriangleExp(;λ0, λ1, period, γ)

A [triangle wave](https://en.wikipedia.org/wiki/Triangle_wave) schedule
with `period` and an exponentially decaying amplitude.
The output conforms to
```text
abs(λ0 - λ1) * Triangle(t) * γ^(t - 1) + min(λ0, λ1)
```
where `Triangle(t)` is `(2 / π) * abs(asin(sin(π * (t - 1) / schedule.period)))` (see [`Triangle`](#)).

# Arguments
- `range0`/`λ0`: the first range endpoint
- `range1`/`λ1`: the second range endpoint
- `period::Integer`: the period
- `decay`/`γ`: the decay rate
"""
struct TriangleExp{T, S<:Integer}
    range0::T
    range1::T
    period::S
    decay::T
end
TriangleExp(;λ0, λ1, period, γ) = TriangleExp(λ0, λ1, period, γ)

startvalue(schedule::TriangleExp) = schedule.tri.λ0
endvalue(schedule::TriangleExp) = schedule.tri.λ1
(schedule::TriangleExp)(t) = _cycle(schedule.range0, schedule.range1,
                                    _tri(t, schedule.period) * schedule.decay^(t - 1))

Base.eltype(::Type{<:TriangleExp{T}}) where T = T
Base.IteratorSize(::Type{<:TriangleExp}) = Base.IsInfinite()

Base.iterate(schedule::TriangleExp, t = 1) = schedule(t), t + 1

Base.axes(::TriangleExp) = (OneToInf(),)

"""
    Sin{T, S<:Integer}(range0, range1, period)
    Sin(;λ0, λ1, period)

A sine wave schedule with `period`.
The output conforms to
```text
abs(λ0 - λ1) * abs(sin(π * (t - 1) / period)) + min(λ0, λ1)
```

# Arguments
- `range0`/`λ0`: the first range endpoint
- `range1`/`λ1`: the second range endpoint
- `period::Integer`: the period
"""
struct Sin{T, S<:Integer}
    range0::T
    range1::T
    period::S
end
Sin(;λ0, λ1, period) = Sin(λ0, λ1, period)

(schedule::Sin)(t) = _cycle(schedule.range0, schedule.range1, _sin(t, schedule.period))

Base.eltype(::Type{<:Sin{T}}) where T = T
Base.IteratorSize(::Type{<:Sin}) = Base.IsInfinite()

Base.iterate(schedule::Sin, t = 1) = schedule(t), t + 1

Base.axes(::Sin) = (OneToInf(),)

"""
    SinDecay2{T, S<:Integer}(range0, range1, period)
    SinDecay2(;λ0, λ1, period)

A sine wave schedule with `period` and half the amplitude each cycle.
The output conforms to
```text
abs(λ0 - λ1) * Sin(t) / (2^floor((t - 1) / period)) + min(λ0, λ1)
```
where `Sin(t)` is `abs(sin(π * (t - 1) / period))` (see [`Sin`](#)).

# Arguments
- `range0`/`λ0`: the first range endpoint
- `range1`/`λ1`: the second range endpoint
- `period::Integer`: the period
"""
struct SinDecay2{T, S<:Integer}
    range0::T
    range1::T
    period::S
end
SinDecay2(;λ0, λ1, period) = SinDecay2(λ0, λ1, period)

(schedule::SinDecay2)(t) = _cycle(schedule.range0, schedule.range1,
                                  _sin(t, schedule.period) / (2^fld(t - 1, schedule.period)))

Base.eltype(::Type{<:SinDecay2{T}}) where T = T
Base.IteratorSize(::Type{<:SinDecay2}) = Base.IsInfinite()

Base.iterate(schedule::SinDecay2, t = 1) = schedule(t), t + 1

Base.axes(::SinDecay2) = (OneToInf(),)

"""
    SinExp{T, S<:Integer}(range0, range1, period, decay)
    SinDecay2(;λ0, λ1, period, γ)

A sine wave schedule with `period` and an exponentially decaying amplitude.
The output conforms to
```text
abs(λ0 - λ1) * Sin(t) * γ^(t - 1) + min(λ0, λ1)
```
where `Sin(t)` is `abs(sin(π * (t - 1) / period))` (see [`Sin`](#)).

# Arguments
- `range0`/`λ0`: the first range endpoint
- `range1`/`λ1`: the second range endpoint
- `period::Integer`: the period
- `decay`/`γ`: the decay rate
"""
struct SinExp{T, S<:Integer}
    range0::T
    range1::T
    period::S
    decay::T
end
SinExp(;λ0, λ1, period, γ) = SinExp(λ0, λ1, period, γ)

(schedule::SinExp)(t) = _cycle(schedule.range0, schedule.range1,
                               _sin(t, schedule.period) * schedule.decay^(t - 1))

Base.eltype(::Type{<:SinExp{T}}) where T = T
Base.IteratorSize(::Type{<:SinExp}) = Base.IsInfinite()

Base.iterate(schedule::SinExp, t = 1) = schedule(t), t + 1

Base.axes(::SinExp) = (OneToInf(),)

"""
    Cos{T, S<:Integer}(range0, range1, period)
    Cos(;λ0, λ1, period)

A cosine annealing schedule
(see ["SGDR: Stochastic Gradient Descent with Warm Restarts"](https://arxiv.org/abs/1608.03983v5))
The output conforms to
```text
abs(λ0 - λ1) * (1 + cos(π * mod(t - 1, period) / period)) / 2 + min(λ0, λ1)
```
This schedule is also referred to as "cosine annealing with warm restarts"
in machine learning literature.

# Arguments
- `range0`/`λ0`: the first range endpoint
- `range1`/`λ1`: the second range endpoint
- `period::Integer`: the period
"""
struct Cos{T, S<:Integer}
    range0::T
    range1::T
    period::S
end
Cos(;λ0, λ1, period) = Cos(λ0, λ1, period)

(schedule::Cos)(t) = _cycle(schedule.range0, schedule.range1,
                            (1 + cos(π * mod(t - 1, schedule.period) / schedule.period)) / 2)

Base.eltype(::Type{<:Cos{T}}) where T = T
Base.IteratorSize(::Type{<:Cos}) = Base.IsInfinite()

Base.iterate(schedule::Cos, t = 1) = schedule(t), t + 1

Base.axes(::Cos) = (OneToInf(),)


text:
"""
    agrees(h, u, p, t)

Determine whether history function evaluates to `u` at time point `t`
for parameters `p`.
"""
function agrees(h, u, p, t)
    # Obtain signatures of h
    sigs = [m.sig for m in methods(h)]

    # Compare evaluation of h at time point t with u
    if any(sig<:Tuple{Any, Any, Any} for sig in sigs)
        return h(p, t) == u
    elseif any(sig<:Tuple{Any, Any, Any, Any} for sig in sigs)
        val = recursivecopy(u)
        h(val, p, t)
        return val == u
    end

    return false
end

"""
    fsal_typeof(integrator::ODEIntegrator)

Return type of FSAL of `integrator`.
"""
function fsal_typeof(integrator::ODEIntegrator{<:OrdinaryDiffEq.OrdinaryDiffEqAlgorithm,IIP,
                                               uType,tType,P,eigenType,tTypeNoUnits,tdirType,ksEltype,
                                               SolType,F,CacheType,O,
                                               FSALType}) where {IIP,uType,
                                                                 tType,P,
                                                                 eigenType,
                                                                 tTypeNoUnits,
                                                                 tdirType,
                                                                 ksEltype,SolType,
                                                                 F,CacheType,O,
                                                                 FSALType}
    return FSALType
end

"""
    assign_expr(::Val{name}, ::Type{T}, ::Type{cache})

Create expression that extracts field `name` of type `T` from cache of type `cache`
to variable `name`.

Hereby u, uprev, uprev2, and function f are updated, if required.
"""
assign_expr(::Val{name}, ::Type, ::Type) where {name} =
    :($name = getfield(cache, $(Meta.quot(name))))

# update matrix exponential
assign_expr(::Val{:expA}, ::Type, ::Type) =
    :(A = f.f1; expA = exp(A*dt))
assign_expr(::Val{:phi1}, ::Type, ::Type{<:OrdinaryDiffEq.NorsettEulerCache}) =
    :(phi1 = ((expA-I)/A))

# update derivative wrappers
assign_expr(::Val{name}, ::Type{<:DiffEqDiffTools.TimeDerivativeWrapper}, ::Type) where name =
    :($name = DiffEqDiffTools.TimeDerivativeWrapper(f, u,p))
assign_expr(::Val{name}, ::Type{<:DiffEqDiffTools.UDerivativeWrapper}, ::Type) where name =
    :($name = DiffEqDiffTools.UDerivativeWrapper(f, t,p))
assign_expr(::Val{name}, ::Type{<:DiffEqDiffTools.TimeGradientWrapper}, ::Type) where name =
    :($name = DiffEqDiffTools.TimeGradientWrapper(
        f,uprev,p))
assign_expr(::Val{name}, ::Type{<:DiffEqDiffTools.UJacobianWrapper}, ::Type) where name =
    :($name = DiffEqDiffTools.UJacobianWrapper(
        f,t,p))

# create new config of Jacobian
assign_expr(::Val{name}, ::Type{<:ForwardDiff.JacobianConfig},
            ::Type) where {name} =
                :($name = OrdinaryDiffEq.build_jac_config(alg, f, uf, du1,
                                                            uprev, u, tmp, dz))
assign_expr(::Val{name}, ::Type{<:ForwardDiff.JacobianConfig},
            ::Type{<:OrdinaryDiffEq.RosenbrockMutableCache}) where {name} =
                :($name = OrdinaryDiffEq.build_jac_config(alg, f, uf, du1,
                                                            uprev, u, tmp, du2))
assign_expr(::Val{name}, ::Type{<:DiffEqDiffTools.JacobianCache},
            ::Type) where {name} =
                :($name = OrdinaryDiffEq.build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz))
assign_expr(::Val{name}, ::Type{<:DiffEqDiffTools.JacobianCache},
            ::Type{<:OrdinaryDiffEq.RosenbrockMutableCache}) where {name} =
                :($name = OrdinaryDiffEq.build_jac_config(alg, f,
                            uf, du1, uprev, u, tmp, du2))
# create new config of Gradient
assign_expr(::Val{name}, ::Type{<:ForwardDiff.DerivativeConfig},
           ::Type) where {name} =
               :($name = OrdinaryDiffEq.build_grad_config(alg, f, tf, du1, t))
assign_expr(::Val{name}, ::Type{<:DiffEqDiffTools.GradientCache},
          ::Type) where {name} =
              :($name = OrdinaryDiffEq.build_grad_config(alg, f, tf, du1, t))

# update implicit RHS
assign_expr(::Val{name}, ::Type{<:OrdinaryDiffEq.ImplicitRHS}, ::Type) where name =
    :($name = OrdinaryDiffEq.ImplicitRHS(f, cache.tmp, t, t, t, cache.dual_cache,p))
assign_expr(::Val{name}, ::Type{<:OrdinaryDiffEq.ImplicitRHS_Scalar}, ::Type) where name =
    :($name = OrdinaryDiffEq.ImplicitRHS_Scalar(f, zero(u), t, t, t,p))
assign_expr(::Val{name}, ::Type{<:OrdinaryDiffEq.RHS_IIF}, ::Type) where name =
    :($name = OrdinaryDiffEq.RHS_IIF(f, cache.tmp, t, t, cache.tmp, cache.dual_cache,p))
assign_expr(::Val{name}, ::Type{<:OrdinaryDiffEq.RHS_IIF_Scalar}, ::Type) where name =
    :($name = OrdinaryDiffEq.RHS_IIF_Scalar(f, zero(u), t, t,
                                            getfield(cache, $(Meta.quot(name))).a,p))

# create new NLsolve differentiable function
assign_expr(::Val{name}, ::Type{<:NLSolversBase.OnceDifferentiable},
            ::Type{<:OrdinaryDiffEq.OrdinaryDiffEqMutableCache}) where name =
                :($name = alg.nlsolve(Val{:init},rhs,u))
assign_expr(::Val{name}, ::Type{<:NLSolversBase.OnceDifferentiable},
            ::Type{<:OrdinaryDiffEq.OrdinaryDiffEqConstantCache}) where name =
                :($name = alg.nlsolve(Val{:init},rhs,uhold))

"""
    build_linked_cache(cache, alg, u, uprev, uprev2, f, t, dt)

Create cache for algorithm `alg` from existing cache `cache` with updated `u`, `uprev`,
`uprev2`, `f`, `t`, and `dt`.
"""
@generated function build_linked_cache(cache, alg, u, uprev, uprev2, f, t, dt,p)
    assignments = [assign_expr(Val{name}(), fieldtype(cache, name), cache)
                   for name in fieldnames(cache) if name ∉ [:u, :uprev, :uprev2, :t, :dt]]

    :($(assignments...); $(DiffEqBase.parameterless_type(cache))($(fieldnames(cache)...)))
end


text:
function AEMParser()
    return AEMParser(())
end

function AEMParser(arg0::DataContext)
    return AEMParser((DataContext,), arg0)
end

function parse(obj::AEMParser, arg0::BufferedReader, arg1::JString)
    return jcall(obj, "parse", AEMFile, (BufferedReader, JString), arg0, arg1)
end

function parse(obj::AEMParser, arg0::InputStream)
    return jcall(obj, "parse", AEMFile, (InputStream,), arg0)
end

function parse(obj::AEMParser, arg0::InputStream, arg1::JString)
    return jcall(obj, "parse", AEMFile, (InputStream, JString), arg0, arg1)
end

function parse(obj::AEMParser, arg0::JString)
    return jcall(obj, "parse", AEMFile, (JString,), arg0)
end

function with_conventions(obj::AEMParser, arg0::IERSConventions)
    return jcall(obj, "withConventions", AEMParser, (IERSConventions,), arg0)
end

function with_data_context(obj::AEMParser, arg0::DataContext)
    return jcall(obj, "withDataContext", AEMParser, (DataContext,), arg0)
end

function with_international_designator(obj::AEMParser, arg0::jint, arg1::jint, arg2::JString)
    return jcall(obj, "withInternationalDesignator", AEMParser, (jint, jint, JString), arg0, arg1, arg2)
end

function with_mission_reference_date(obj::AEMParser, arg0::AbsoluteDate)
    return jcall(obj, "withMissionReferenceDate", AEMParser, (AbsoluteDate,), arg0)
end

function with_mu(obj::AEMParser, arg0::jdouble)
    return jcall(obj, "withMu", AEMParser, (jdouble,), arg0)
end

function with_simple_eop(obj::AEMParser, arg0::jboolean)
    return jcall(obj, "withSimpleEOP", AEMParser, (jboolean,), arg0)
end


text:
module CORE

function uni(n) #QR decomposition. Q being the unitary matrix.
    H = randn(n,n) + 1im*randn(n,n)
    Q,R = qr(H) ### QR decomposition
    #return Q*Q.' #symmetric
    return Q #arbiturary
    #return diagm(Complex{Float64}[1,1,1,1])
end

#s = uni(4)

#call unitary matrix such that the random values are unanimous in all expression
# s here is 4 by 4 energy independent
# better run separately to for the sake of minor fixes in other codes while pinning s

function ground(E,x,y,z,s) # returns the ground energy

    e_phi(x,y) = exp.([x,y,0.5im*Float64(z),0]) # array list of phase-exponentials
    ephi = diagm(e_phi(1im*x/2.,1im*y/2.))
    ephiInv = inv(ephi)
    Asqr = 0.5*(1-E/sqrt(E^2+1))
    Asqrem = 0.5*(1+E/sqrt(E^2+1))
    Q =  Asqr*ephiInv*s*ephi  +  Asqrem*ephi*s.'*ephiInv
    return -0.5*log(det(Q))
end

function current1(E,x,y,z,s) # returns the current w.r.t the first phase change

    e_phi(x,y) = exp.([x,y,0.5im*Float64(z),0]) # array list of phase-exponentials
    ephi = diagm(e_phi(1im*x/2.,1im*y/2.))
    ephiInv = inv(ephi)
    Asqr = 0.5*(1-E/sqrt(E^2+1))
    Asqrem = 0.5*(1+E/sqrt(E^2+1))

    phiProj1 = diagm(Complex{Float64}[1im/2,0,0,0])

    Q =  Asqr*ephiInv*s*ephi  +  Asqrem*ephi*s.'*ephiInv
    Qinv = inv(Q)
    comm1 = ephiInv*s*ephi*phiProj1  -   phiProj1*ephiInv*s*ephi
    Qdphi1 =  Asqr*comm1 + Asqrem*comm1.'

    return -0.5*trace(Qinv*Qdphi1)
end

function current2(E,x,y,z,s) ## returns the current w.r.t the second phase change

    e_phi(x,y) = exp.([x,y,0.5im*Float64(z),0]) # array list of phase-exponentials
    ephi = diagm(e_phi(1im*x/2.,1im*y/2.))
    ephiInv = inv(ephi)
    Asqr = 0.5*(1-E/sqrt(E^2+1))
    Asqrem = 0.5*(1+E/sqrt(E^2+1))

    phiProj2 = diagm(Complex{Float64}[0,1im/2,0,0])

    Q =  Asqr*ephiInv*s*ephi  +  Asqrem*ephi*s.'*ephiInv
    Qinv = inv(Q)
    comm2 = ephiInv*s*ephi*phiProj2  -   phiProj2*ephiInv*s*ephi
    Qdphi2 =  Asqr*comm2 + Asqrem*comm2.'

    return -0.5*trace(Qinv*Qdphi2)
end

function invar(E,x,y,z,s) # invariant Q^{-1}.dQ/de.Q^{-1}.dQ/d(phi_1).Q^{-1}.dQ/d(phi_2), antisymmetrised upon 1,2
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)
    sphi = ephiInv*s*ephi
    Asqr = 0.5*(1.-E/sqrt(E^2+1.))
    Asqrem = 0.5*(1.+E/sqrt(E^2+1.))
    Q =  Asqr*sphi  +  Asqrem*sphi.'
    Qinv = inv(Q)
    Qde = 1/(2*(sqrt(E^2+1))^3)*(sphi.'-sphi)

    phiProj1 = diagm(Complex{Float64}[0.5*1im,0.,0.,0.])
    phiProj2 = diagm(Complex{Float64}[0.,0.5*1im,0.,0.])

    comm1 = sphi*phiProj1  -   phiProj1*sphi
    comm2 = sphi*phiProj2  -   phiProj2*sphi
    #comm(n) = ephiInv*s*ephi*phiProj(n)  -   phiProj(n)*ephiInv*s*ephi
    Qdphi1 =  Asqr*comm1 + Asqrem*comm1.'
    Qdphi2 =  Asqr*comm2 + Asqrem*comm2.'

    #Qdphi(n) =  Asqr*comm(n) + Asqrem*comm(n)'

    inv_ = Qinv*Qde*Qinv*Qdphi1*Qinv*Qdphi2  -  Qinv*Qde*Qinv*Qdphi2*Qinv*Qdphi1
    return -0.25*trace(inv_)
end

function SE(E_,x,y,z,scale_,s_)
    E = Float64(E_)
    scale = Float64(scale_)
    s = Array{Complex{Float64}}(s_)
    F::Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}} = eigfact(s)
    d::Array{Complex{Float64},1} = F[:values]
    v::Array{Complex{Float64},2} = F[:vectors]
    sgn::Array{Float64,2} = diagm(sign.(sin.(0.5*angle.(d))))
    e0 = scale*diagm(cos.(0.5*angle.(d)))*sgn
    gamma = 2*scale*diagm(sin.(0.5*angle.(d)))*sgn
    pole = e0 + 0.5*im*gamma
    S = v*((im*E*eye(4)-pole)/(im*E*eye(4)-conj(pole)))*v'
    return S
end

function LargeE(E_,x,y,z,scale_,s_) # invariant Q^{-1}.dQ/de.Q^{-1}.dQ/d(phi_1).Q^{-1}.dQ/d(phi_2), antisymmetrised upon 1,2
    E = Float64(E_)
    scale = Float64(scale_)
    s = Array{Complex{Float64}}(s_)
    F::Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}}= eigfact(s)
    d::Array{Complex{Float64},1} = F[:values]
    v::Array{Complex{Float64},2} = F[:vectors]
    sgn::Array{Float64,2} = diagm(sign.(sin.(0.5*angle.(d))))
    e0 = scale*diagm(cos.(0.5*angle.(d)))*sgn
    gamma = 2*scale*diagm(sin.(0.5*angle.(d)))*sgn
    pole = e0 + 0.5*im*gamma
    S = v*((im*E*eye(4)-pole)/(im*E*eye(4)-conj(pole)))*v'
    Sem = v*((-im*E*eye(4)-pole)/(-im*E*eye(4)-conj(pole)))*v'
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)
    sphi = ephiInv*S*ephi
    semphi = ephiInv*Sem*ephi
    Semde = ephiInv*v*(gamma/(-im*E*eye(4)-conj(pole))^2)*v'*ephi

    phiProj1 = diagm(Complex{Float64}[0.5*1im,0.,0.,0.])
    phiProj2 = diagm(Complex{Float64}[0.,0.5*1im,0.,0.])

    comm1 = sphi*phiProj1  -   phiProj1*sphi   #dS/dphi1
    comm2 = sphi*phiProj2  -   phiProj2*sphi
    commem1 = semphi*phiProj1  -   phiProj1*semphi #E -> -E
    commem2 = semphi*phiProj2  -   phiProj2*semphi

    inv_::Array{Complex{Float64},2} = Semde'*sphi*commem1'*comm2 - Semde'*sphi*commem2'*comm1
    return -0.5*trace(inv_)
end

function test(E_,x,y,z,scale_,s_)
    E = Float64(E_)
    scale = Float64(scale_)
    s = Array{Complex{Float64}}(s_)
    F = eigfact(s)
    d::Array{Complex{Float64},1} = F[:values]
    v::Array{Complex{Float64},2} = F[:vectors]
    sgn::Array{Float64,2} = diagm(sign.(sin.(0.5*angle.(d))))
    e0 = scale*diagm(cos.(0.5*angle.(d)))*sgn
    gamma = 2*scale*diagm(sin.(0.5*angle.(d)))*sgn
    pole = e0 + 0.5*im*gamma
    S = v*((im*E*eye(4)-pole)/(im*E*eye(4)-conj(pole)))*v'
    Semde = v*(-gamma/((-im*E*eye(4)-conj(pole))^2))*v'
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)
    phiProj1 = diagm(Complex{Float64}[0.5*1im,0.,0.,0.])
    phiProj2 = diagm(Complex{Float64}[0.,0.5*1im,0.,0.])
    sphi = ephiInv*S*ephi
    Semdephi = ephiInv*Semde*ephi
    Semdagphi = inv(sphi)
    ans = trace(Semdephi'*Semdagphi*(Semdagphi*phiProj1-phiProj1*Semdagphi)*(sphi*phiProj2-phiProj2*sphi) - Semdephi'*Semdagphi*(Semdagphi*phiProj2-phiProj2*Semdagphi)*(sphi*phiProj1-phiProj1*sphi))
    return -0.5*ans
end

function invarE(E::Float64, x::Float64, y::Float64, z::Float64, scale::Float64, s::Array{Complex{Float64},2}) # invariant Q^{-1}.dQ/de.Q^{-1}.dQ/d(phi_1).Q^{-1}.dQ/d(phi_2), antisymmetrised upon 1,2
    F::Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}} = eigfact(s)
    d::Array{Complex{Float64},1} = F[:values]::Array{Complex{Float64},1}
    v::Array{Complex{Float64},2} = F[:vectors]::Array{Complex{Float64},2}
    sgn::Array{Float64,2} = diagm(sign.(sin.(0.5*angle.(d))))
    e0 = scale*diagm(cos.(0.5*angle.(d)))*sgn
    gamma = 2.*scale*diagm(sin.(0.5*angle.(d)))*sgn
    pole = e0 + 0.5*im*gamma
    S = v*((im*E*eye(4)-pole)/(im*E*eye(4)-conj(pole)))*v'
    Sem = v*((-im*E*eye(4)-pole)/(-im*E*eye(4)-conj(pole)))*v'  #energy depend S transpose has to change E -> -E
    Sde = v*(-gamma/((im*E*eye(4)-conj(pole))^2))*v'
    Semde = v*(gamma/((-im*E*eye(4)-conj(pole))^2))*v'  #dS_{-E}/de
    ephi = diagm(exp.([0.5im*x,0.5im*y,0.5im*z,0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)
    sphi = ephiInv*S*ephi
    semphi = ephiInv*Sem*ephi  #sphi with -E
    sdephi = ephiInv*Sde*ephi
    semdephi = ephiInv*Semde*ephi   #sdephi with -E

    Asqr = 0.5*(1.-E/sqrt(E^2+1.))
    Asqrem = 0.5*(1.+E/sqrt(E^2+1.))
    Q =  Asqr*sphi  +  Asqrem*semphi.'
    Qinv = inv(Q)
    Qde = 1./(2.*(sqrt(E^2+1.))^3)*(semphi.'-sphi) + Asqr*sdephi + Asqrem*semdephi.'

    phiProj1 = diagm(Complex{Float64}[0.5*1im,0.,0.,0.])
    phiProj2 = diagm(Complex{Float64}[0.,0.5*1im,0.,0.])

    comm1 = sphi*phiProj1  -   phiProj1*sphi   #dS/dphi1
    comm2 = sphi*phiProj2  -   phiProj2*sphi
    commem1 = semphi*phiProj1  -   phiProj1*semphi #E -> -E
    commem2 = semphi*phiProj2  -   phiProj2*semphi

    Qdphi1 =  Asqr*comm1 + Asqrem*commem1.'
    Qdphi2 =  Asqr*comm2 + Asqrem*commem2.'

    #Qdphi(n) =  Asqr*comm(n) + Asqrem*comm(n)'

    inv_::Array{Complex{Float64},2} = Qinv*Qde*Qinv*Qdphi1*Qinv*Qdphi2  -  Qinv*Qde*Qinv*Qdphi2*Qinv*Qdphi1
    return -0.25*trace(inv_)
end

function boundary(E_,x,y,z,scale_,s_) # invariant Q^{-1}.dQ/de.Q^{-1}.dQ/d(phi_1).Q^{-1}.dQ/d(phi_2), antisymmetrised upon 1,2
    E = Float64(E_)
    scale = Float64(scale_)
    s = Array{Complex{Float64}}(s_)
    F::Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}} = eigfact(s)
    d::Array{Complex{Float64},1} = F[:values]
    v::Array{Complex{Float64},2} = F[:vectors]
    sgn::Array{Float64,2} = diagm(sign.(sin.(0.5*angle.(d))))
    e0 = scale*diagm(cos.(0.5*angle.(d)))*sgn
    gamma = 2*scale*diagm(sin.(0.5*angle.(d)))*sgn
    pole = e0 + 0.5*im*gamma
    S = v*((im*E*eye(4)-pole)/(im*E*eye(4)-conj(pole)))*v'
    Sem = v*((-im*E*eye(4)-pole)/(-im*E*eye(4)-conj(pole)))*v'  #energy depende transpose has to change E -> -E
    Sde = v*(-gamma/((im*E*eye(4)-conj(pole))^2))*v'
    Semde = v*(gamma/((-im*E*eye(4)-conj(pole))^2))*v'  #dS_{-E}/de
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)
    sphi = ephiInv*S*ephi
    semphi = ephiInv*Sem*ephi  #sphi with -E
    sdephi = ephiInv*Sde*ephi
    semdephi = ephiInv*Semde*ephi   #sdephi with -E
    Asqr = 0.5*(1.-E/sqrt(E^2+1.))
    Asqrem = 0.5*(1.+E/sqrt(E^2+1.))
    Q =  Asqr*sphi  +  Asqrem*semphi.'
    Qinv = inv(Q)
    Qde = 1/(2*(sqrt(E^2+1))^3)*(semphi.'-sphi) + Asqr*sdephi + Asqrem*semdephi.'

    phiProj1 = diagm(Complex{Float64}[0.5*1im,0.,0.,0.])
    phiProj2 = diagm(Complex{Float64}[0.,0.5*1im,0.,0.])

    comm1::Array{Complex{Float64},2} = sphi*phiProj1  -   phiProj1*sphi    #dS/dphi1
    comm2::Array{Complex{Float64},2} = sphi*phiProj2  -   phiProj2*sphi
    commem1 = semphi*phiProj1  -   phiProj1*semphi #E -> -E
    commem2 = semphi*phiProj2  -   phiProj2*semphi
    comme1::Array{Complex{Float64},2} = sdephi*phiProj1  -   phiProj1*sdephi  #d^2 S/ dedphi1
    comme2::Array{Complex{Float64},2} = sdephi*phiProj2  -   phiProj2*sdephi

    #comm(n) = ephiInv*s*ephi*phiProj(n)  -   phiProj(n)*ephiInv*s*ephi
    Qdphi1 =  Asqr*comm1 + Asqrem*commem1.'
    Qdphi2 =  Asqr*comm2 + Asqrem*commem2.'

    Bound::Array{Complex{Float64},2} = Asqr*(-Qinv*Qdphi2*Qinv*(sdephi*phiProj1+phiProj1*sdephi) + Qinv*comme2*phiProj1 + Qinv*phiProj1*comme2) - Asqr*(-Qinv*Qdphi1*Qinv*(sdephi*phiProj2+phiProj2*sdephi) + Qinv*comme1*phiProj2 + Qinv*phiProj2*comme1)

    #Qdphi(n) =  Asqr*comm(n) + Asqrem*comm(n)'

    inv_::Array{Complex{Float64},2} = 2*Bound
    return -0.25*trace(inv_)
end

function invarEcomp1(E_,x,y,z,scale_,s_) # invariant Q^{-1}.dQ/de.Q^{-1}.dQ/d(phi_1).Q^{-1}.dQ/d(phi_2), antisymmetrised upon 1,2
    E = Float64(E_)
    scale = Float64(scale_)
    s = Array{Complex{Float64}}(s_)
    F::Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}} = eigfact(s)
    d::Array{Complex{Float64},1} = F[:values]
    v::Array{Complex{Float64},2} = F[:vectors]
    sgn::Array{Float64,2} = diagm(sign.(sin.(0.5*angle.(d))))
    e0 = scale*diagm(cos.(0.5*angle.(d)))*sgn
    gamma = 2*scale*diagm(sin.(0.5*angle.(d)))*sgn
    pole = e0 + 0.5*im*gamma
    S = v*((im*0.5*(1-E^2)*eye(4)-pole*E)/(im*0.5*(1-E^2)*eye(4)-E*conj(pole)))*v'
    Sem = v*((-im*0.5*(1-E^2)*eye(4)-pole*E)/(-im*0.5*(1-E^2)*eye(4)-E*conj(pole)))*v'  #energy depend S transpose has to change E -> -E
    Sde = v*(-(gamma*E^2)/((im*0.5*(1-E^2)*eye(4)-E*conj(pole))^2))*v'
    Semde = v*((gamma*E^2)/((-im*0.5*(1-E^2)*eye(4)-E*conj(pole))^2))*v'  #dS_{-E}/de
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)
    sphi = ephiInv*S*ephi
    semphi = ephiInv*Sem*ephi  #sphi with -E
    sdephi = ephiInv*Sde*ephi
    semdephi = ephiInv*Semde*ephi   #sdephi with -E

    Asqr = E^2/(E^2+1)
    Asqrem = 1./(1+E^2)
    Q =  Asqr*sphi  +  Asqrem*semphi.'
    Qinv = inv(Q)
    Qde = 4*E^3/((1+E^2)^3)*(semphi.'-sphi) + Asqr*sdephi + Asqrem*semdephi.'

    phiProj1 = diagm(Complex{Float64}[0.5*1im,0.,0.,0.])
    phiProj2 = diagm(Complex{Float64}[0.,0.5*1im,0.,0.])

    comm1 = sphi*phiProj1  -   phiProj1*sphi   #dS/dphi1
    comm2 = sphi*phiProj2  -   phiProj2*sphi
    commem1 = semphi*phiProj1  -   phiProj1*semphi #E -> -E
    commem2 = semphi*phiProj2  -   phiProj2*semphi

    Qdphi1 =  Asqr*comm1 + Asqrem*commem1.'
    Qdphi2 =  Asqr*comm2 + Asqrem*commem2.'

    #Qdphi(n) =  Asqr*comm(n) + Asqrem*comm(n)'

    inv_::Array{Complex{Float64},2} = Qinv*Qde*Qinv*Qdphi1*Qinv*Qdphi2  -  Qinv*Qde*Qinv*Qdphi2*Qinv*Qdphi1
    return -0.25*trace(inv_)::Complex{Float64}
end

function invarEcomp2(E_,x,y,z,scale_,s_) # invariant Q^{-1}.dQ/de.Q^{-1}.dQ/d(phi_1).Q^{-1}.dQ/d(phi_2), antisymmetrised upon 1,2
    E = Float64(E_)
    scale = Float64(scale_)
    s = Array{Complex{Float64}}(s_)
    F::Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}} = eigfact(s)
    d::Array{Complex{Float64},1} = F[:values]
    v::Array{Complex{Float64},2} = F[:vectors]
    sgn::Array{Float64,2} = diagm(sign.(sin.(0.5*angle.(d))))
    e0 = scale*diagm(cos.(0.5*angle.(d)))*sgn
    gamma = 2*scale*diagm(sin.(0.5*angle.(d)))*sgn
    pole = e0 + 0.5*im*gamma
    S = v*((im*0.5*(E^2-1)*eye(4)-pole*E)/(im*0.5*(E^2-1)*eye(4)-E*conj(pole)))*v'
    Sem = v*((-im*0.5*(E^2-1)*eye(4)-pole*E)/(-im*0.5*(E^2-1)*eye(4)-E*conj(pole)))*v'  #energy depend S transpose has to change E -> -E
    Sde = v*(-(gamma*E^2)/((im*0.5*(E^2-1)*eye(4)-E*conj(pole))^2))*v'
    Semde = v*((gamma*E^2)/((-im*0.5*(E^2-1)*eye(4)-E*conj(pole))^2))*v'  #dS_{-E}/de
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)
    sphi = ephiInv*S*ephi
    semphi = ephiInv*Sem*ephi  #sphi with -E
    sdephi = ephiInv*Sde*ephi
    semdephi = ephiInv*Semde*ephi   #sdephi with -E

    Asqr = 1./(1+E^2)
    Asqrem = E^2/(1.+E^2)
    Q =  Asqr*sphi  +  Asqrem*semphi.'
    Qinv = inv(Q)
    Qde = 4*E^3/((1+E^2)^3)*(semphi.'-sphi) + Asqr*sdephi + Asqrem*semdephi.'

    phiProj1 = diagm(Complex{Float64}[0.5*1im,0.,0.,0.])
    phiProj2 = diagm(Complex{Float64}[0.,0.5*1im,0.,0.])

    comm1 = sphi*phiProj1  -   phiProj1*sphi   #dS/dphi1
    comm2 = sphi*phiProj2  -   phiProj2*sphi
    commem1 = semphi*phiProj1  -   phiProj1*semphi #E -> -E
    commem2 = semphi*phiProj2  -   phiProj2*semphi

    Qdphi1 =  Asqr*comm1 + Asqrem*commem1.'
    Qdphi2 =  Asqr*comm2 + Asqrem*commem2.'

    #Qdphi(n) =  Asqr*comm(n) + Asqrem*comm(n)'

    inv_::Array{Complex{Float64},2} = Qinv*Qde*Qinv*Qdphi1*Qinv*Qdphi2  -  Qinv*Qde*Qinv*Qdphi2*Qinv*Qdphi1
    return -0.25*trace(inv_)::Complex{Float64}
end

function boundarycomp1(E_,x,y,z,scale_,s_) # invariant Q^{-1}.dQ/de.Q^{-1}.dQ/d(phi_1).Q^{-1}.dQ/d(phi_2), antisymmetrised upon 1,2
    E = Float64(E_)
    scale = Float64(scale_)
    s = Array{Complex{Float64}}(s_)
    F::Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}} = eigfact(s)
    d::Array{Complex{Float64},1} = F[:values]
    v::Array{Complex{Float64},2} = F[:vectors]
    sgn::Array{Float64,2} = diagm(sign.(sin.(0.5*angle.(d))))
    e0 = scale*diagm(cos.(0.5*angle.(d)))*sgn
    gamma = 2*scale*diagm(sin.(0.5*angle.(d)))*sgn
    pole = e0 + 0.5*im*gamma
    S = v*((im*0.5*(1-E^2)*eye(4)-pole*E)/(im*0.5*(1-E^2)*eye(4)-E*conj(pole)))*v'
    Sem = v*((-im*0.5*(1-E^2)*eye(4)-pole*E)/(-im*0.5*(1-E^2)*eye(4)-E*conj(pole)))*v'  #energy depend S transpose has to change E -> -E
    Sde = v*(-(gamma*E^2)/((im*0.5*(1-E^2)*eye(4)-E*conj(pole))^2))*v'
    Semde = v*((gamma*E^2)/((-im*0.5*(1-E^2)*eye(4)-E*conj(pole))^2))*v'  #dS_{-E}/de
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)
    sphi = ephiInv*S*ephi
    semphi = ephiInv*Sem*ephi  #sphi with -E
    sdephi = ephiInv*Sde*ephi
    semdephi = ephiInv*Semde*ephi   #sdephi with -E

    Asqr = E^2/(E^2+1)
    Asqrem = 1./(1+E^2)
    Q =  Asqr*sphi  +  Asqrem*semphi.'
    Qinv = inv(Q)
    Qde = E^3/(2*(0.5*(1+E^2))^3)*(semphi.'-sphi) + Asqr*sdephi + Asqrem*semdephi.'

    phiProj1 = diagm(Complex{Float64}[0.5*1im,0.,0.,0.])
    phiProj2 = diagm(Complex{Float64}[0.,0.5*1im,0.,0.])

    comm1::Array{Complex{Float64},2} = sphi*phiProj1  -   phiProj1*sphi    #dS/dphi1
    comm2::Array{Complex{Float64},2} = sphi*phiProj2  -   phiProj2*sphi
    commem1 = semphi*phiProj1  -   phiProj1*semphi #E -> -E
    commem2 = semphi*phiProj2  -   phiProj2*semphi
    comme1::Array{Complex{Float64},2} = sdephi*phiProj1  -   phiProj1*sdephi  #d^2 S/ dedphi1
    comme2::Array{Complex{Float64},2} = sdephi*phiProj2  -   phiProj2*sdephi

    #comm(n) = ephiInv*s*ephi*phiProj(n)  -   phiProj(n)*ephiInv*s*ephi
    Qdphi1 =  Asqr*comm1 + Asqrem*commem1.'
    Qdphi2 =  Asqr*comm2 + Asqrem*commem2.'

    Bound::Array{Complex{Float64},2} = Asqr*(-Qinv*Qdphi2*Qinv*(sdephi*phiProj1+phiProj1*sdephi) + Qinv*comme2*phiProj1 + Qinv*phiProj1*comme2) - Asqr*(-Qinv*Qdphi1*Qinv*(sdephi*phiProj2+phiProj2*sdephi) + Qinv*comme1*phiProj2 + Qinv*phiProj2*comme1)

    #Qdphi(n) =  Asqr*comm(n) + Asqrem*comm(n)'

    inv_::Array{Complex{Float64},2} = 2*Bound
    return -0.25*trace(inv_)
end

function boundarycomp2(E_,x,y,z,scale_,s_) # invariant Q^{-1}.dQ/de.Q^{-1}.dQ/d(phi_1).Q^{-1}.dQ/d(phi_2), antisymmetrised upon 1,2
    E = Float64(E_)
    scale = Float64(scale_)
    s = Array{Complex{Float64}}(s_)
    F::Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}} = eigfact(s)
    d::Array{Complex{Float64},1} = F[:values]
    v::Array{Complex{Float64},2} = F[:vectors]
    sgn::Array{Float64,2} = diagm(sign.(sin.(0.5*angle.(d))))
    e0 = scale*diagm(cos.(0.5*angle.(d)))*sgn
    gamma = 2*scale*diagm(sin.(0.5*angle.(d)))*sgn
    pole = e0 + 0.5*im*gamma
    S = v*((im*0.5*(E^2-1)*eye(4)-pole*E)/(im*0.5*(E^2-1)*eye(4)-E*conj(pole)))*v'
    Sem = v*((-im*0.5*(E^2-1)*eye(4)-pole*E)/(-im*0.5*(E^2-1)*eye(4)-E*conj(pole)))*v'  #energy depend S transpose has to change E -> -E
    Sde = v*(-(gamma*E^2)/((im*0.5*(E^2-1)*eye(4)-E*conj(pole))^2))*v'
    Semde = v*((gamma*E^2)/((-im*0.5*(E^2-1)*eye(4)-E*conj(pole))^2))*v'  #dS_{-E}/de
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)
    sphi = ephiInv*S*ephi
    semphi = ephiInv*Sem*ephi  #sphi with -E
    sdephi = ephiInv*Sde*ephi
    semdephi = ephiInv*Semde*ephi   #sdephi with -E

    Asqr = 1./(1+E^2)
    Asqrem = E^2/(1.+E^2)
    Q =  Asqr*sphi  +  Asqrem*semphi.'
    Qinv = inv(Q)
    Qde = E^3/(2*(0.5*(1+E^2))^3)*(semphi.'-sphi) + Asqr*sdephi + Asqrem*semdephi.'

    phiProj1 = diagm(Complex{Float64}[0.5*1im,0.,0.,0.])
    phiProj2 = diagm(Complex{Float64}[0.,0.5*1im,0.,0.])

    comm1::Array{Complex{Float64},2} = sphi*phiProj1  -   phiProj1*sphi    #dS/dphi1
    comm2::Array{Complex{Float64},2} = sphi*phiProj2  -   phiProj2*sphi
    commem1 = semphi*phiProj1  -   phiProj1*semphi #E -> -E
    commem2 = semphi*phiProj2  -   phiProj2*semphi
    comme1::Array{Complex{Float64},2} = sdephi*phiProj1  -   phiProj1*sdephi  #d^2 S/ dedphi1
    comme2::Array{Complex{Float64},2} = sdephi*phiProj2  -   phiProj2*sdephi

    #comm(n) = ephiInv*s*ephi*phiProj(n)  -   phiProj(n)*ephiInv*s*ephi
    Qdphi1 =  Asqr*comm1 + Asqrem*commem1.'
    Qdphi2 =  Asqr*comm2 + Asqrem*commem2.'

    Bound::Array{Complex{Float64},2} = Asqr*(-Qinv*Qdphi2*Qinv*(sdephi*phiProj1+phiProj1*sdephi) + Qinv*comme2*phiProj1 + Qinv*phiProj1*comme2) - Asqr*(-Qinv*Qdphi1*Qinv*(sdephi*phiProj2+phiProj2*sdephi) + Qinv*comme1*phiProj2 + Qinv*phiProj2*comme1)

    #Qdphi(n) =  Asqr*comm(n) + Asqrem*comm(n)'

    inv_::Array{Complex{Float64},2} = 2*Bound
    return -0.25*trace(inv_)
end

function invarEindep(E::Float64, x::Float64, y::Float64, z::Float64, scale::Float64, s::Array{Complex{Float64},2}) # invariant Q^{-1}.dQ/de.Q^{-1}.dQ/d(phi_1).Q^{-1}.dQ/d(phi_2), antisymmetrised upon 1,2
    F::Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}} = eigfact(s)
    d::Array{Complex{Float64},1} = F[:values]::Array{Complex{Float64},1}
    v::Array{Complex{Float64},2} = F[:vectors]::Array{Complex{Float64},2}
    sgn::Array{Float64,2} = diagm(sign.(sin.(0.5*angle.(d))))
    e0 = diagm(cos.(0.5*angle.(d)))*sgn
    gamma = 2.*diagm(sin.(0.5*angle.(d)))*sgn
    pole = e0 + 0.5*im*gamma
    S = v*((im*E*eye(4)-pole)/(im*E*eye(4)-conj(pole)))*v'
    Sem = v*((-im*E*eye(4)-pole)/(-im*E*eye(4)-conj(pole)))*v'  #energy depend S transpose has to change E -> -E
    Sde = v*(-gamma/((im*E*eye(4)-conj(pole))^2))*v'
    Semde = v*(gamma/((-im*E*eye(4)-conj(pole))^2))*v'  #dS_{-E}/de
    ephi = diagm(exp.([0.5im*x,0.5im*y,0.5im*z,0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)
    sphi = ephiInv*S*ephi
    semphi = ephiInv*Sem*ephi  #sphi with -E
    sdephi = ephiInv*Sde*ephi
    semdephi = ephiInv*Semde*ephi   #sdephi with -E

    Asqr = 0.5*(1.-scale*E/sqrt(scale^2*E^2+1.))
    Asqrem = 0.5*(1.+scale*E/sqrt(scale^2*E^2+1.))
    Q =  Asqr*sphi  +  Asqrem*semphi.'
    Qinv = inv(Q)
    Qde = Asqr*sdephi + Asqrem*semdephi.' + 1./(2.*(sqrt(scale^2*E^2+1.))^3)*(semphi.'-sphi)*scale

    phiProj1 = diagm(Complex{Float64}[0.5*1im,0.,0.,0.])
    phiProj2 = diagm(Complex{Float64}[0.,0.5*1im,0.,0.])

    comm1 = sphi*phiProj1  -   phiProj1*sphi   #dS/dphi1
    comm2 = sphi*phiProj2  -   phiProj2*sphi
    commem1 = semphi*phiProj1  -   phiProj1*semphi #E -> -E
    commem2 = semphi*phiProj2  -   phiProj2*semphi

    Qdphi1 =  Asqr*comm1 + Asqrem*commem1.'
    Qdphi2 =  Asqr*comm2 + Asqrem*commem2.'

    #Qdphi(n) =  Asqr*comm(n) + Asqrem*comm(n)'

    inv_::Array{Complex{Float64},2} = Qinv*Qde*Qinv*Qdphi1*Qinv*Qdphi2  -  Qinv*Qde*Qinv*Qdphi2*Qinv*Qdphi1
    return -0.25*trace(inv_)
end

function invarcut(ie,x,y,z,s) # invariant Q^{-1}.dQ/de.Q^{-1}.dQ/d(phi_1).Q^{-1}.dQ/d(phi_2), antisymmetrised upon 1,2
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)
    sphi = ephiInv*s*ephi
    Asqr_p = 0.5*(1.-ie/sqrt(Float64(ie)^2-1.)) # A^2_e analy connti. e-> ie cut positive real, equivalent to A^2_-e negative real
    Asqr_m = 0.5*(1.+ie/sqrt(Float64(ie)^2-1.)) # A^2_e nagative real, equiv to A^2_-e posi. real
    Q_p =  Asqr_p*sphi  +  Asqr_m*sphi.' #Q analy connti. e-> ie cut positive real
    Q_m =  Asqr_m*sphi  +  Asqr_p*sphi.' #Q analy connti. e-> ie cut negative real
    Qp_inv = inv(Q_p)
    Qm_inv = inv(Q_m)
    Qp_de = 0.5*im/((Float64(ie)^2-1.)^1.5)*(sphi.'-sphi) #positive real dQ/de
    Qm_de = -0.5*im/((Float64(ie)^2-1.)^1.5)*(sphi.'-sphi)

    phiProj1 = diagm(Complex{Float64}[1im/2.,0.,0.,0.])
    phiProj2 = diagm(Complex{Float64}[0.,1im/2.,0.,0.])

    comm1 = sphi*phiProj1  -   phiProj1*sphi
    comm2 = sphi*phiProj2  -   phiProj2*sphi
    #comm(n) = ephiInv*s*ephi*phiProj(n)  -   phiProj(n)*ephiInv*s*ephi
    Qp_dphi1 =  Asqr_p*comm1 + Asqr_m*comm1.'
    Qm_dphi1 =  Asqr_m*comm1 + Asqr_p*comm1.'
    Qp_dphi2 =  Asqr_p*comm2 + Asqr_m*comm2.'
    Qm_dphi2 =  Asqr_m*comm2 + Asqr_p*comm2.'

    inv_p = Qp_inv*Qp_de*Qp_inv*Qp_dphi1*Qp_inv*Qp_dphi2  -  Qp_inv*Qp_de*Qp_inv*Qp_dphi2*Qp_inv*Qp_dphi1
    inv_m = Qm_inv*Qm_de*Qm_inv*Qm_dphi1*Qm_inv*Qm_dphi2  -  Qm_inv*Qm_de*Qm_inv*Qm_dphi2*Qm_inv*Qm_dphi1
    return -0.25*trace(inv_p-inv_m)
end

function avrginvar(x,y,z,s) #the average invariant part at the boundary: 0.5I_1 * S^\dag * 0.5I_2 S
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)
    sphi = ephiInv*s*ephi
    phiProj1 = diagm(Complex{Float64}[0.5*im,0.,0.,0.])
    phiProj2 = diagm(Complex{Float64}[0.,0.5*im,0.,0.])
    avrg_=(phiProj1*sphi'*phiProj2*sphi-phiProj2*sphi'*phiProj1*sphi)
    return -0.5*trace(avrg_)/(2.*pi)
end

function ssconj(x,y,z,s) # sum over all the positive ABS energies
    e_phi(x,y) = exp.([x,y,0.5im*Float64(z),0]) # array list of phase-exponentials
    ephi = diagm(e_phi(1im*x/2.,1im*y/2.))
    ephiInv = inv(ephi)

    Lambda = ephiInv*s*ephi*ephi*conj(s)*ephiInv
    chi = real(-0.5im*log.(eigvals(Lambda)))
    positive = chi[chi.>0]
    return sum(cos.(positive))+(2-ndims(positive))
end

function LowestABS(x,y,z,s) # Lowest ABS energies
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)

    Lambda = ephiInv*s*ephi*ephi*conj(s)*ephiInv
    chi = real(-0.5im*log.(eigvals(Lambda)))
    positive = chi[chi.>0]
    ABS = cos.(max(positive...))
    return ABS
end

function HighestABS(x,y,z,s) # Lowest ABS energies
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)

    Lambda = ephiInv*s*ephi*ephi*conj(s)*ephiInv
    chi = real(-0.5im*log.(eigvals(Lambda)))
    positive = chi[chi.>0]
    ABS = cos.(min(positive...))
    return ABS
end

function HighestIs1(x,y,z,s) # Lowest ABS energies
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)

    Lambda = ephiInv*s*ephi*ephi*conj(s)*ephiInv
    chi = real(-0.5im*log.(eigvals(Lambda)))
    positive = chi[chi.>0]
    ABS = cos.(min(positive...))
    if abs(ABS-1.)<0.00003
        return 1.
    else return 0.
    end
end

function LowestIs0(x,y,z,s) # Lowest ABS energies
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)

    Lambda = ephiInv*s*ephi*ephi*conj(s)*ephiInv
    chi = real(-0.5im*log.(eigvals(Lambda)))
    positive = chi[chi.>0]
    ABS = cos.(max(positive...))
    if abs(ABS-0)<0.02
        return 1.
    else return 0.
    end
end

function HisL(x,y,z,s) # Lowest ABS energies
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)

    Lambda = ephiInv*s*ephi*ephi*conj(s)*ephiInv
    chi = real(-0.5im*log.(eigvals(Lambda)))
    positive = chi[chi.>0]
    ABSL = cos.(max(positive...))
    ABSH = cos.(min(positive...))
    if abs(ABSL-ABSH)<0.02
        return 1.
    else return 0.
    end
end

function invarSum(x,y,z,s) #This should return the same value as function "intInvar" yet faster
    e_phi(x,y) = exp.([x,y,0.5im*Float64(z),0]) # array list of phase-exponentials
    ephi = diagm(e_phi(1im*x/2.,1im*y/2.))
    ephiInv = inv(ephi)
    sphi = ephiInv*s*ephi
    phiProj1 = diagm(Complex{Float64}[1im/2,0,0,0])
    phiProj2 = diagm(Complex{Float64}[0,1im/2,0,0])

    Lambda = sphi*conj(sphi)
    D,V = eig(Lambda)
    dLambda1 = 2*sphi*phiProj1*conj(sphi)-phiProj1*Lambda-Lambda*phiProj1
    dLambda2 = 2*sphi*phiProj2*conj(sphi)-phiProj2*Lambda-Lambda*phiProj2

    summation = 0         #This is summing over eigenvalues and eigenvectors' indeces
    for i = 1:4
        for j = 1:4
            if j == i
                continue
            end
            summation += 2*log.(D[i])/((D[i]-D[j])^2)*((V'[i,:].'*dLambda1*V[:,j])*(V'[j,:].'*dLambda2*V[:,i])
                        -(V'[i,:].'*dLambda2*V[:,j])*(V'[j,:].'*dLambda1*V[:,i]))
                        +1/(D[j]*(D[j]-D[i]))*((V'[i,:].'*dLambda1*V[:,j])*(V'[j,:].'*dLambda2*V[:,i])
                                    -(V'[i,:].'*dLambda2*V[:,j])*(V'[j,:].'*dLambda1*V[:,i]))
                        -(1-D[i]/D[j])*(1/(D[j]-D[i]))*((V'[i,:].'*(phiProj1-ephi*s.'*ephiInv*phiProj1*ephi*conj(s)*ephiInv)*V[:,j])*(V'[j,:].'*dLambda2*V[:,i])
                        -(V'[i,:].'*(phiProj2-ephi*s.'*ephiInv*phiProj2*ephi*conj(s)*ephiInv)*V[:,j])*(V'[j,:].'*dLambda1*V[:,i]))
        end
    end
    return -0.25*summation/(2.*pi)

end

function invarSum1(x,y,z,s) #This should return the same value as function "intSum" yet faster
    ephi = diagm(exp.([0.5im*Float64(x),0.5im*Float64(y),0.5im*Float64(z),0.])) # diag matrix of phase-exponentials
    ephiInv = inv(ephi)
    sphi = ephiInv*s*ephi
    phiProj1 = diagm(Complex{Float64}[1im/2,0,0,0])
    phiProj2 = diagm(Complex{Float64}[0,1im/2,0,0])

    Lambda = sphi*conj(sphi)
    D,V = eig(Lambda)
    dLambda1 = 2*sphi*phiProj1*conj(sphi)-phiProj1*Lambda-Lambda*phiProj1
    dLambda2 = 2*sphi*phiProj2*conj(sphi)-phiProj2*Lambda-Lambda*phiProj2

    P1 = V'*phiProj1*V
    P2 = V'*phiProj2*V
    Pbar1 = -conj(V')*sphi'*phiProj1*sphi*conj(V) #this Pbar[j,i] = <jbar| 0.5im*I |ibar> = -<i| s^T * 0.5im*I*s^* |j>
    Pbar2 = -conj(V')*sphi'*phiProj2*sphi*conj(V)

       #This is summing over eigenvalues and eigenvectors' indeces
    summation1 = 0
    for m = 1:4
        for n = 1:4
            if n == m
                continue
            end
            summation1 +=
             ((-2*log.(D[n]))/((D[m]-D[n])^2)*4*D[m]*D[n]*P1[m,n]*P2[n,m]+(2*log.(D[m]))/((D[m]-D[n])^2)*((D[m]+D[n])^2)*P1[m,n]*P2[n,m]+(2*log.(D[m])-2*log.(D[n]))/((D[m]-D[n])^2)*(2*D[m]*D[n]+2*D[m]^2)*Pbar1[n,m]*P2[n,m]) - ((-2*log.(D[n]))/((D[m]-D[n])^2)*4*D[m]*D[n]*P2[m,n]*P1[n,m]+(2*log.(D[m]))/((D[m]-D[n])^2)*((D[m]+D[n])^2)*P2[m,n]*P1[n,m]+(2*log.(D[m])-2*log.(D[n]))/((D[m]-D[n])^2)*(2*D[m]*D[n]+2*D[m]^2)*Pbar2[n,m]*P1[n,m]) + (-4*((D[m]+D[n])/(D[m]-D[n]))*P1[m,n]*P2[n,m]+(1/(D[n]*(D[n]-D[m])))*((D[m]+D[n])^2+4*D[m]*D[n])*P2[n,m]*Pbar1[n,m]) - (-4*((D[m]+D[n])/(D[m]-D[n]))*P2[m,n]*P1[n,m]+(1/(D[n]*(D[n]-D[m])))*((D[m]+D[n])^2+4*D[m]*D[n])*P1[n,m]*Pbar2[n,m])
        end
    end
    return -0.25*(summation1)/(2.*pi)

end

function invarSumP(x,y,z,s) # Invariant (Q^{-1}dQ)^3 summed over only Positive energies
    e_phi(x,y) = exp.([x,y,0.5im*Float64(z),0]) # array list of phase-exponentials
    ephi = diagm(e_phi(1im*x/2.,1im*y/2.))
    ephiInv = inv(ephi)
    phiProj1 = diagm(Complex{Float64}[1im/2,0,0,0])
    phiProj2 = diagm(Complex{Float64}[0,1im/2,0,0])

    Lambda = ephiInv*s*ephi*ephi*conj(s)*ephiInv
    D,V = eig(Lambda)
    chi = real(-0.5im*log.(D))
    dLambda1 = 2*ephiInv*s*ephi*phiProj1*ephi*conj(s)*ephiInv-phiProj1*Lambda-Lambda*phiProj1
    dLambda2 = 2*ephiInv*s*ephi*phiProj2*ephi*conj(s)*ephiInv-phiProj2*Lambda-Lambda*phiProj2

    summation = 0
    summation2 = 0         #This is summing over eigenvalues and eigenvectors' indeces
    for i = 1:4
        if chi[i]<0
            continue
        end
        for j = 1:4
            if j == i
                continue
            end
            summation += 4*log.(D[i])/((D[i]-D[j])^2)*((V'[i,:].'*dLambda1*V[:,j])*(V'[j,:].'*dLambda2*V[:,i])
                        -(V'[i,:].'*dLambda2*V[:,j])*(V'[j,:].'*dLambda1*V[:,i]))
        end
    end
    for m = 1:4
        if chi[m]>0
            continue
        end
        for n = 1:4
            if chi[n]<0
                continue
            end
            if n == m
                continue
            end
            summation2 += 1/(D[n]*(D[n]-D[m]))*((V'[m,:].'*dLambda1*V[:,n])*(V'[n,:].'*dLambda2*V[:,m])
                          -(V'[m,:].'*dLambda2*V[:,n])*(V'[n,:].'*dLambda1*V[:,m]))
                          -(1-D[m]/D[n])*(1/(D[n]-D[m]))*((V'[m,:].'*(phiProj1-ephi*s.'*ephiInv*phiProj1*ephi*conj(s)*ephiInv)*V[:,n])*(V'[n,:].'*dLambda2*V[:,m])
                          -(V'[m,:].'*(phiProj2-ephi*s.'*ephiInv*phiProj2*ephi*conj(s)*ephiInv)*V[:,n])*(V'[n,:].'*dLambda1*V[:,m]))
        end
    end
    return -0.25*(summation2)/(2.*pi)

end

function BerryP(x,y,z,s) # Berry curvature only positive ABS energies
    e_phi(x,y) = exp.([x,y,0.5im*Float64(z),0]) # array list of phase-exponentials
    ephi = diagm(e_phi(1im*x/2.,1im*y/2.))
    ephiInv = inv(ephi)
    phiProj1 = diagm(Complex{Float64}[1im/2,0,0,0])
    phiProj2 = diagm(Complex{Float64}[0,1im/2,0,0])

    Lambda = ephiInv*s*ephi*ephi*conj(s)*ephiInv
    D,V = eig(Lambda)
    chi = real(-0.5im*log.(D))
    dLambda1 = 2*ephiInv*s*ephi*phiProj1*ephi*conj(s)*ephiInv-phiProj1*Lambda-Lambda*phiProj1
    dLambda2 = 2*ephiInv*s*ephi*phiProj2*ephi*conj(s)*ephiInv-phiProj2*Lambda-Lambda*phiProj2

    summation = 0         #This is summing over eigenvalues and eigenvectors' indeces
    for i = 1:4
        if chi[i]<0
            continue
        end
        for j = 1:4
            if j == i
                continue
            end
            summation += 2im*pi/((D[i]-D[j])^2)*((V'[i,:].'*dLambda1*V[:,j])*(V'[j,:].'*dLambda2*V[:,i])
                        -(V'[i,:].'*dLambda2*V[:,j])*(V'[j,:].'*dLambda1*V[:,i]))
        end
    end
    return -0.25*summation/(2.*pi)

end

function BerryBranch(x,y,z,s) # Berry curvature only positive ABS energies
    e_phi(x,y) = exp.([x,y,0.5im*Float64(z),0]) # array list of phase-exponentials
    ephi = diagm(e_phi(1im*x/2.,1im*y/2.))
    ephiInv = inv(ephi)
    sphi = ephiInv*s*ephi
    phiProj1 = diagm(Complex{Float64}[1im/2,0,0,0])
    phiProj2 = diagm(Complex{Float64}[0,1im/2,0,0])

    Lambda = sphi*conj(sphi)
    D,V = eig(Lambda)
    dLambda1 = 2*sphi*phiProj1*conj(sphi)-phiProj1*Lambda-Lambda*phiProj1
    dLambda2 = 2*sphi*phiProj2*conj(sphi)-phiProj2*Lambda-Lambda*phiProj2

    P1 = V'*phiProj1*V
    P2 = V'*phiProj2*V
    Pbar1 = -conj(V')*sphi'*phiProj1*sphi*conj(V)
    Pbar2 = -conj(V')*sphi'*phiProj2*sphi*conj(V)

    summation = 0         #This is summing over eigenvalues and eigenvectors' indeces
    for i = 1:4
        if imag(D[i])<=0
            continue
        end
        for j = 1:4
            if j == i
                continue
            end # the 2 factor in "2* 2im* pi" comes from the factor in "2log<di | di>"
            summation += 2*2im*pi*((4*D[i]*D[j]*Pbar1[j,i]*Pbar2[i,j] + (D[i]+D[j])^2*P1[i,j]*P2[j,i] + 2*(D[i]*D[j]+D[j]^2)*P1[i,j]*Pbar2[i,j] + 2*(D[i]^2+D[i]*D[j])*P2[j,i]*Pbar1[j,i]) - (4*D[i]*D[j]*Pbar2[j,i]*Pbar1[i,j] + (D[i]+D[j])^2*P2[i,j]*P1[j,i] + 2*(D[i]*D[j]+D[j]^2)*P2[i,j]*Pbar1[i,j] + 2*(D[i]^2+D[i]*D[j])*P1[j,i]*Pbar2[j,i]))/((D[i]-D[j])^2)
        end
    end
    return -0.25*summation/(2.*pi)

end

function invarSumSec(x,y,z,s) #This should return the same value as function "intInvar" yet faster
    e_phi(x,y) = exp.([x,y,0.5im*Float64(z),0]) # array list of phase-exponentials
    ephi = diagm(e_phi(1im*x/2.,1im*y/2.))
    ephiInv = inv(ephi)
    sphi = ephiInv*s*ephi
    phiProj1 = diagm(Complex{Float64}[1im/2,0,0,0])
    phiProj2 = diagm(Complex{Float64}[0,1im/2,0,0])

    Lambda = sphi*conj(sphi)
    D,V = eig(Lambda)
    chi = real(-0.5im*log.(D))
    dLambda1 = 2*sphi*phiProj1*conj(sphi)-phiProj1*Lambda-Lambda*phiProj1
    dLambda2 = 2*sphi*phiProj2*conj(sphi)-phiProj2*Lambda-Lambda*phiProj2

    summation = 0         #This is summing over eigenvalues and eigenvectors' indeces
    for i = 1:4
        for j = 1:4
            if j == i
                continue
            end
            summation += 1/(D[j]*(D[j]-D[i]))*((V'[i,:].'*dLambda1*V[:,j])*(V'[j,:].'*dLambda2*V[:,i])
                        -(V'[i,:].'*dLambda2*V[:,j])*(V'[j,:].'*dLambda1*V[:,i]))   -(1-D[i]/D[j])*(1/(D[j]-D[i]))*((V'[i,:].'*(phiProj1-ephi*s.'*ephiInv*phiProj1*ephi*conj(s)*ephiInv)*V[:,j])*(V'[j,:].'*dLambda2*V[:,i])
                          -(V'[i,:].'*(phiProj2-ephi*s.'*ephiInv*phiProj2*ephi*conj(s)*ephiInv)*V[:,j])*(V'[j,:].'*dLambda1*V[:,i]))
        end
    end
    return -0.25*summation/(2.*pi)
end

function invarSumP1(x,y,z,s) # Invariant (Q^{-1}dQ)^3 summed over only Positive energies
    e_phi(x,y) = exp.([x,y,0.5im*Float64(z),0]) # array list of phase-exponentials
    ephi = diagm(e_phi(1im*x/2.,1im*y/2.))
    ephiInv = inv(ephi)
    phiProj1 = diagm(Complex{Float64}[1im/2,0,0,0])
    phiProj2 = diagm(Complex{Float64}[0,1im/2,0,0])

    Lambda = ephiInv*s*ephi*ephi*conj(s)*ephiInv
    D,V = eig(Lambda)
    chi = real(-0.5im*log.(D))
    dLambda1 = 2*ephiInv*s*ephi*phiProj1*ephi*conj(s)*ephiInv-phiProj1*Lambda-Lambda*phiProj1
    dLambda2 = 2*ephiInv*s*ephi*phiProj2*ephi*conj(s)*ephiInv-phiProj2*Lambda-Lambda*phiProj2

    summation = 0       #This is summing over eigenvalues and eigenvectors' indeces
    for i = 1:4
        if chi[i]>0
            continue
        end
        for j = 1:4
            if j == i
                continue
            end
            summation += 2*(2*log.(D[i])/((D[i]-D[j])^2)*((V'[i,:].'*dLambda1*V[:,j])*(V'[j,:].'*dLambda2*V[:,i])
                        -(V'[i,:].'*dLambda2*V[:,j])*(V'[j,:].'*dLambda1*V[:,i]))
                        +1/(D[j]*(D[j]-D[i]))*((V'[i,:].'*dLambda1*V[:,j])*(V'[j,:].'*dLambda2*V[:,i])
                          -(V'[i,:].'*dLambda2*V[:,j])*(V'[j,:].'*dLambda1*V[:,i]))
                          -(1-D[i]/D[j])*(1/(D[j]-D[i]))*((V'[i,:].'*(phiProj1-ephi*s.'*ephiInv*phiProj1*ephi*conj(s)*ephiInv)*V[:,j])*(V'[j,:].'*dLambda2*V[:,i])
                          -(V'[i,:].'*(phiProj2-ephi*s.'*ephiInv*phiProj2*ephi*conj(s)*ephiInv)*V[:,j])*(V'[j,:].'*dLambda1*V[:,i])))
        end
    end
    return -0.25*summation/(2.*pi)

end

using QuadGK
using Cuba

function ChernS(z::Float64,s::Array{Complex{Float64},2})
    result::Array{Float64,1}, = cuhre((x::Array{Float64,1},f::Array{Float64,1})->f[1] = ((2.*x[1]^2-2.*x[1]+1.)/((x[1]-x[1]^2)^2))*2*pi*real(invar((2.*x[1]-1.)/(x[1]-x[1]^2),2.*pi*x[2],2.*pi*x[3],z,s)::Complex{Float64}),3,abstol = 1e-12, reltol = 1e-10)

    return result[1] #the result should /2pi for energy integral and *(2pi)^2 for [0,1] phase compactify. So, overall 2pi
end

function ChernL(z::Float64,s::Array{Complex{Float64},2})
    result::Array{Float64,1}, = cuhre((x::Array{Float64,1},f::Array{Float64,1})->f[1] = 4*pi^2*real(avrginvar(2.*pi*x[1],2.*pi*x[2],z,s)::Complex{Float64}),2,abstol = 1e-8, reltol = 1e-7)

    return result[1] #the result should /2pi for energy integral and *(2pi)^2 for [0,1] phase compactify. So, overall 2pi
end

function ChernSum(z::Float64,s::Array{Complex{Float64},2})
    result::Array{Float64,1}, = cuhre((x::Array{Float64,1},f::Array{Float64,1})->f[1] = 4*pi^2*real(invarSum1(2.*pi*x[1],2.*pi*x[2],z,s)::Complex{Float64}),2,abstol = 1e-8, reltol = 1e-7)

    return result[1] #the result should /2pi for energy integral and *(2pi)^2 for [0,1] phase compactify. So, overall 2pi
end

function ChernBerry(z::Float64,s::Array{Complex{Float64},2})
    result::Array{Float64,1}, = cuhre((x::Array{Float64,1},f::Array{Float64,1})->f[1] = 4*pi^2*real(BerryBranch(2.*pi*x[1],2.*pi*x[2],z,s)::Complex{Float64}),2,abstol = 1e-8, reltol = 1e-7)

    return result[1] #the result should /2pi for energy integral and *(2pi)^2 for [0,1] phase compactify. So, overall 2pi
end

function ChernCut(z::Float64,s::Array{Complex{Float64},2})
    result::Array{Float64,1}, = cuhre((x::Array{Float64,1},f::Array{Float64,1})->f[1] = 1./((1.-x[1])^2)*2*pi*real(invarcut(1./(1.-x[1]),2.*pi*x[2],2.*pi*x[3],z,s)::Complex{Float64}) - 1./((1.-x[1])^2)*2*pi*real(invarcut(1+1e-10,2.*pi*x[2],2.*pi*x[3],z,s))*1e-10/(1./(1.-x[1])-1),3,abstol = 1e-8, reltol = 1e-7)

    return result[1] #the result should /2pi for energy integral and *(2pi)^2 for [0,1] phase compactify. So, overall 2pi
end

function ChernE(z::Float64,scale::Float64,s::Array{Complex{Float64},2})
    result1::Array{Float64,1}, = cuhre((x::Array{Float64,1},f::Array{Float64,1})->f[1] = 2*pi*real(invarEcomp2(x[1],2.*pi*x[2],2.*pi*x[3],z,scale,s)::Complex{Float64}*(x[1]^2+1)/(2.*x[1]^2)),3,atol = 1e-8, rtol = 1e-7)

    result2::Array{Float64,1}, = cuhre((x::Array{Float64,1},f::Array{Float64,1})->f[1] = 2*pi*real(invarEcomp1(x[1],2.*pi*x[2],2.*pi*x[3],z,scale,s)::Complex{Float64}*(x[1]^2+1)/(2.*x[1]^2)),3,abstol = 1e-8, reltol = 1e-7)

    return result1[1] + result2[1] #the result should /2pi for energy integral and *(2pi)^2 for [0,1] phase compactify. So, overall 2pi
end

function ChernEindep(z::Float64,scale::Float64,s::Array{Complex{Float64},2})
    result::Array{Float64,1}, = cuhre((x::Array{Float64,1},f::Array{Float64,1})->f[1] = 20*2*pi*real(invarEindep(20*x[1]-10,2.*pi*x[2],2.*pi*x[3],z,scale,s)::Complex{Float64}),3,abstol = 1e-8, reltol = 1e-7)

    return result[1] #the result should /2pi for energy integral and *(2pi)^2 for [0,1] phase compactify. So, overall 2pi
end

function ChernBoundE(z::Float64,scale::Float64,s::Array{Complex{Float64},2})
    result1::Array{Float64,1}, = cuhre((x::Array{Float64,1},f::Array{Float64,1})->f[1] = 2*pi*real(boundarycomp2(x[1],2.*pi*x[2],2.*pi*x[3],z,scale,s)::Complex{Float64}*(x[1]^2+1)/(2.*x[1]^2)),3,atol = 1e-8, rtol = 1e-7)

    result2::Array{Float64,1}, = cuhre((x::Array{Float64,1},f::Array{Float64,1})->f[1] = 2*pi*real(boundarycomp1(x[1],2.*pi*x[2],2.*pi*x[3],z,scale,s)::Complex{Float64}*(x[1]^2+1)/(2.*x[1]^2)),3,abstol = 1e-8, reltol = 1e-7)

    return result1[1] + result2[1] #the result should /2pi for energy integral and *(2pi)^2 for [0,1] phase compactify. So, overall 2pi
end

function intGrnd(x,y,z,s)
    intg, err = quadgk(e->ground(e,x,y,z,s), -50, 50, atol=0.0000001)
    return intg/(2.*pi)
end

function intCurrent1(x,y,z,s)
    intg, err = quadgk(e->current1(e,x,y,z,s), -50, 50, atol=0.0000001)
    return intg/(2.*pi)
end

function intCurrent2(x,y,z,s)
    intg, err = quadgk(e->current2(e,x,y,z,s), -50, 50, atol=0.0000001)
    return intg/(2.*pi)
end

function intInvar(x,y,z,s) # integrate invariant (Q^{-1}dQ)^3 w.r.t energy for phases (x,y), constant S matrix
    intgS::Complex{Float64}, err = quadgk(E->invar(E,x,y,z,s), -2., 2., atol=0.00000001)::Tuple{Complex{Float64},Float64}
    intgL1::Complex{Float64}, err = quadgk(E->invar(E,x,y,z,s), -50., -2., atol=0.00000001)::Tuple{Complex{Float64},Float64}
    intgL2::Complex{Float64}, err = quadgk(E->invar(E,x,y,z,s), 2., 50., atol=0.00000001)::Tuple{Complex{Float64},Float64}
    return (intgS+intgL1+intgL2)/(2.*pi)
end

function intInvarES(x::Float64, y::Float64, z::Float64, scale::Float64, s::Array{Complex{Float64},2}) # integrate invariant of the response function at all scale energy dependent scale<1
    intg1::Float64, = quadgk(E::Float64->real(invarE(E,x,y,z,scale,s)), -2.*scale, 2.*scale)::Tuple{Float64,Float64}
    intg2::Float64, = quadgk(E->real(invarE(E,x,y,z,scale,s)), -2., -2.*scale)::Tuple{Float64,Float64}
    intg3::Float64, = quadgk(E->real(invarE(E,x,y,z,scale,s)), -100., -2.)::Tuple{Float64,Float64}
    intg4::Float64, = quadgk(E->real(invarE(E,x,y,z,scale,s)), 2.*scale, 2.)::Tuple{Float64,Float64}
    intg5::Float64, = quadgk(E->real(invarE(E,x,y,z,scale,s)), 2., 100.)::Tuple{Float64,Float64}
    return (intg1+intg2+intg3+intg4+intg5)/(2.*pi)
end

function intInvarEB(x::Float64, y::Float64, z::Float64, scale::Float64, s::Array{Complex{Float64},2}) # integrate invariant of the response function at all scale energy dependent large scale>1
    intg1::Float64, = quadgk(E::Float64->real(invarE(E,x,y,z,scale,s)), -2., 2., atol=0.0000001)::Tuple{Float64,Float64}
    intg2::Float64, = quadgk(E->real(invarE(E,x,y,z,scale,s)), -2.*scale, -2., atol=0.0000001)::Tuple{Float64,Float64}
    intg3::Float64, = quadgk(E->real(invarE(E,x,y,z,scale,s)), 2.,2.*scale, atol=0.0000001)::Tuple{Float64,Float64}
    intg4::Float64, = quadgk(E->real(invarE(E,x,y,z,scale,s)), -100.*scale,-2.*scale, atol=0.0000001)::Tuple{Float64,Float64}
    intg5::Float64, = quadgk(E->real(invarE(E,x,y,z,scale,s)), 2.*scale, 100.*scale, atol=0.0000001)::Tuple{Float64,Float64}
    return (intg1+intg2+intg3+intg4+intg5)/(2.*pi)

end

function intBoundE(x,y,z,scale,s) # integrate boundary at all scale
    if scale < 1.
        intg1, err = quadgk(E->boundary(E,x,y,z,scale,s), -2.*scale, 2.*scale)::Tuple{Complex{Float64},Float64}
        intg2, err = quadgk(E->boundary(E,x,y,z,scale,s), -2., -2.*scale)::Tuple{Complex{Float64},Float64}
        intg3, err = quadgk(E->boundary(E,x,y,z,scale,s), -100., -2.)::Tuple{Complex{Float64},Float64}
        intg4, err = quadgk(E->boundary(E,x,y,z,scale,s), 2.*scale, 2.)::Tuple{Complex{Float64},Float64}
        intg5, err = quadgk(E->boundary(E,x,y,z,scale,s), 2., 100.)::Tuple{Complex{Float64},Float64}
        sum = (intg1+intg2+intg3+intg4+intg5)/(2.*pi)
    else
        intg1, err = quadgk(E->boundary(E,x,y,z,scale,s), -2., 2.)::Tuple{Complex{Float64},Float64}
        intg2, err = quadgk(E->boundary(E,x,y,z,scale,s), -2.*scale, -2.)::Tuple{Complex{Float64},Float64}
        intg3, err = quadgk(E->boundary(E,x,y,z,scale,s), 2.,2.*scale)::Tuple{Complex{Float64},Float64}
        intg4, err = quadgk(E->boundary(E,x,y,z,scale,s), -100.*scale,-2.*scale)::Tuple{Complex{Float64},Float64}
        intg5, err = quadgk(E->boundary(E,x,y,z,scale,s), 2.*scale, 100.*scale)::Tuple{Complex{Float64},Float64}
        sum = (intg1+intg2+intg3+intg4+intg5)/(2.*pi)
    end
    return sum
end

function intInvarEcomp(x::Float64,y::Float64,z::Float64,scale::Float64,s::Array{Complex{Float64},2}) # integrate invarE at all scale
    intg1::Complex{Float64}, = quadgk(E::Float64->invarEcomp1(E,x,y,z,scale,s)*(E^2+1)/(-2.*E^2), 0., 1.)::Tuple{Complex{Float64},Float64}
    intg2::Complex{Float64}, = quadgk(E::Float64->invarEcomp2(E,x,y,z,scale,s)*(E^2+1)/(2.*E^2), 0., 1.)::Tuple{Complex{Float64},Float64}
    return (-intg1+intg2)/(2.*pi)
end

function intBoundEcomp(x,y,z,scale,s) # integrate invarE at all scale
    intg1::Complex{Float64}, err = quadgk(E->boundarycomp1(E,x,y,z,scale,s)*(E^2+1)/(-2.*E^2), 0., 1.)::Tuple{Complex{Float64},Float64}
    intg2::Complex{Float64}, err = quadgk(E->boundarycomp2(E,x,y,z,scale,s)*(E^2+1)/(2.*E^2), 0., 1.)::Tuple{Complex{Float64},Float64}

    return (-intg1+intg2)/(2.*pi)
end

function intInvarEL(x,y,z,scale,s) # integrate invarE at large scale
    intg1::Complex{Float64}, err = quadgk(E->invarE(E,x,y,z,scale,s), 50., 100000., atol=0.00000001)::Tuple{Complex{Float64},Float64}
    intg2::Complex{Float64}, err = quadgk(E->invarE(E,x,y,z,scale,s), -100000., -50., atol=0.00000001)::Tuple{Complex{Float64},Float64}
    return (intg1+intg2)/(2.*pi)
end

function intBoundEL(x,y,z,scale,s) # integrate boundary at large scale
    intg1::Complex{Float64}, err = quadgk(E->boundary(E,x,y,z,scale,s), 50., 100000., atol=0.00000001)::Tuple{Complex{Float64},Float64}
    intg2::Complex{Float64}, err = quadgk(E->boundary(E,x,y,z,scale,s), -100000., -50., atol=0.00000001)::Tuple{Complex{Float64},Float64}

    return (intg1+intg2)/(2.*pi)
end

function intLargeE(x,y,z,scale,s) # integrate largeE
    intg::Complex{Float64}, err = quadgk(E->LargeE(E,x,y,z,scale,s), 20., 10000., atol=0.00000001)::Tuple{Complex{Float64},Float64}
    return -intg/(2.*pi)
end

function intInvarCut(x,y,z,s)
    intg, err = quadgk(ie->invarcut(ie,x,y,z,s), 1.00000001, 50.00000001, atol=0.00000001)
    return intg/(2.*pi)
end

function intInvarerr(x,y,z,s) # integrate invariant (Q^{-1}dQ)^3 w.r.t energy for phases (x,y), limit(-50, 50)
    intg, err = quadgk(e->invar(e,x,y,z,s), -50, 50, atol=0.0000001)
    return err/(2.*pi)
end

function invarSumB(x,y,z,s) #This should return the same value as function "intInvar" yet faster
    e_phi(x,y) = exp.([x,y,0.5im*Float64(z),0]) # array list of phase-exponentials
    ephi = diagm(e_phi(1im*x/2.,1im*y/2.))
    ephiInv = inv(ephi)
    phiProj1 = diagm(Complex{Float64}[1im/2,0,0,0])
    phiProj2 = diagm(Complex{Float64}[0,1im/2,0,0])

    Lambda = ephiInv*s*ephi*ephi*conj(s)*ephiInv
    D,V = eig(Lambda)
    dLambda1 = 2*ephiInv*s*ephi*phiProj1*ephi*conj(s)*ephiInv-phiProj1*Lambda-Lambda*phiProj1
    dLambda2 = 2*ephiInv*s*ephi*phiProj2*ephi*conj(s)*ephiInv-phiProj2*Lambda-Lambda*phiProj2

    dL1 = V'*dLambda1*V
    dL2 = V'*dLambda2*V

    summation = 0         #This is summing over eigenvalues and eigenvectors' indeces
    for i = 1:4
        for j = 1:4
            if j == i
                continue
            end
            summation += 2*log.(D[i])/((D[i]-D[j])^2)*(dL1[i,j]*dL2[j,i]
                        -dL2[i,j]*dL1[j,i])
        end
    end
    return -0.25*summation/(2.*pi)

end

#intInvar(1.2,1.3,s)

GrndGrid(xx,yy,z,s) = map(x->intGrnd(x...,z,s),Iterators.product(yy,xx))

Current1Grid(xx,yy,z,s) = map(x->intCurrent1(x...,z,s),Iterators.product(yy,xx))

Current2Grid(xx,yy,z,s) = map(x->intCurrent2(x...,z,s),Iterators.product(yy,xx))

InvarGrid(xx,yy,z,s) = map(x->intInvar(x...,z,s),Iterators.product(yy,xx))
#@time InvarGrid = pmap(x->intInvar(x...,s),Iterators.product(1:3,1:4))
InvarCutGrid(xx,yy,z,s) = map(x->intInvarCut(x...,z,s),Iterators.product(yy,xx))

InvarGriderr(xx,yy,z,s) = map(x->intInvarerr(x...,z,s),Iterators.product(yy,xx))

#intInvar(1.2,1.3)
avrgInvarGrid(xx,yy,z,s) = map(x->avrginvar(x...,z,s),Iterators.product(yy,xx))

LambdaGrid(xx,yy,z,s) = map(x->ssconj(x...,z,s),Iterators.product(yy,xx))
LowestABSGrid(xx,yy,z,s) = map(x->LowestABS(x...,z,s),Iterators.product(yy,xx))
HighestABSGrid(xx,yy,z,s) = map(x->HighestABS(x...,z,s),Iterators.product(yy,xx))
HighestIs1Grid(xx,yy,z,s) = map(x->HighestIs1(x...,z,s),Iterators.product(yy,xx))
LowestIs0Grid(xx,yy,z,s) = map(x->LowestIs0(x...,z,s),Iterators.product(yy,xx))
HisLGrid(xx,yy,z,s) = map(x->HisL(x...,z,s),Iterators.product(yy,xx))

InvarSumGrid(xx,yy,z,s) = map(x->invarSum(x...,z,s),Iterators.product(yy,xx))
InvarSumPGrid(xx,yy,z,s) = map(x->invarSumP(x...,z,s),Iterators.product(yy,xx))
InvarSum1Grid(xx,yy,z,s) = map(x->invarSum1(x...,z,s),Iterators.product(yy,xx))
InvarSumBGrid(xx,yy,z,s) = map(x->invarSumB(x...,z,s),Iterators.product(yy,xx))

invarSumP1Grid(xx,yy,z,s) = map(x->invarSumP1(x...,z,s),Iterators.product(yy,xx))

BerryPGrid(xx,yy,z,s) = map(x->BerryP(x...,z,s),Iterators.product(yy,xx))
BerryBranchGrid(xx,yy,z,s) = map(x->BerryBranch(x...,z,s),Iterators.product(yy,xx))
InvarSumSecGrid(xx,yy,z,s) = map(x->invarSumSec(x...,z,s),Iterators.product(yy,xx))

IntInvarESGrid(xx,yy,z,scale,s) = map(x->intInvarES(x...,z,scale,s),Iterators.product(yy,xx))
IntInvarEBGrid(xx,yy,z,scale,s) = map(x->intInvarEB(x...,z,scale,s),Iterators.product(yy,xx))
IntBoundEGrid(xx,yy,z,scale,s) = map(x->intBoundE(x...,z,scale,s),Iterators.product(yy,xx))

IntInvarEcompGrid(xx,yy,z,scale,s) = map(x->intInvarEcomp(x...,z,scale,s),Iterators.product(yy,xx))

export avrginvar, avrgInvarGrid, BerryBranchGrid, BerryPGrid, Current1Grid, Current2Grid, GrndGrid, HighestABSGrid, HighestIs1Grid, intBoundE, intBoundEcomp, IntBoundEGrid, intBoundEL, intInvar, intInvarE, intInvarE1, intInvarEcomp, intInvarEL, IntInvarEBGrid, IntInvarESGrid, IntInvarEcompGrid, intLargeE, invar, invarcut, InvarCutGrid, invarE, invarEcomp1, invarEcomp2, invarEGrid, InvarGrid, InvarGriderr, invarSum1, InvarSum1Grid, InvarSumBGrid, InvarSumGrid, InvarSumPGrid, invarSumP1Grid, InvarSumSecGrid, LambdaGrid, LowestABSGrid, uni, s, SE, HisLGrid, LowestIs0Grid, invarEindep, ChernE, ChernEindep, ChernBoundE

end

#using QDPHI

max_stars_repo_path: all in one/core.jl

max_issues_repo_path: all in one/core.jl

max_forks_repo_path: all in one/core.jl


text:
include("common.jl")

example = """
0,9 -> 5,9
8,0 -> 0,8
9,4 -> 3,4
2,2 -> 2,1
7,0 -> 7,4
6,4 -> 2,0
0,9 -> 2,9
3,4 -> 1,4
0,0 -> 8,8
5,5 -> 8,2
"""

struct Line
    start_x::Int64
    start_y::Int64
    end_x::Int64
    end_y::Int64
end

function Base.parse(::Type{Line}, line::AbstractString)
    l_start, l_end = split(line, " -> ")
    start_x, start_y = split(l_start, ",")
    end_x, end_y = split(l_end, ",")
    return Line(map(x -> parse(Int64, x), (start_x, start_y, end_x, end_y))...)
end

function Base.parse(::Type{Vector{Line}}, input::AbstractString)
    lines = split(input, "\n")
    result = Line[]
    for line in lines
        isempty(line) && continue
        push!(result, parse(Line, line))
    end
    return result
end

horizontal(line::Line) = line.start_y == line.end_y
vertical(line::Line) = line.start_x == line.end_x

rectilinear(line::Line) = horizontal(line) || vertical(line)

function range_index_x(line::Line)
    ix1 = line.start_x + 1
    ix2 = line.end_x + 1
    return ix1 < ix2 ? (ix1:ix2) : (ix1:-1:ix2)
end

function range_index_y(line::Line)
    iy1 = line.start_y + 1
    iy2 = line.end_y + 1
    return iy1 < iy2 ? (iy1:iy2) : (iy1:-1:iy2)
end

function make_grid(lines::AbstractVector{Line})
    max_x = 0
    max_y = 0
    for line in lines
        max_x = max(max_x, line.start_x)
        max_x = max(max_x, line.end_x)
        max_y = max(max_y, line.start_y)
        max_y = max(max_y, line.end_y)
    end
    # Note that we have a zero-based indexing in terms of Line, so add one.
    return zeros(Int64, (max_x + 1, max_y + 1))
end

function populate_rectilinear!(grid::Matrix{Int64}, lines::AbstractVector{Line})
    for line in lines
        if horizontal(line)
            i_y = line.start_y + 1
            @inbounds for i_x in range_index_x(line)
                grid[i_x, i_y] += 1
            end
        elseif vertical(line)
            i_x = line.start_x + 1
            @inbounds for i_y in range_index_y(line)
                grid[i_x, i_y] += 1
            end
        end
        # Line gets skipped if it isn't rectilinear.
    end
    return grid
end

function populate_rectilinear(lines::AbstractVector{Line})
    return populate_rectilinear!(make_grid(lines), lines)
end

function get_ans(input::AbstractString)
    lines = parse(Vector{Line}, input)
    grid = populate_rectilinear(lines)
    return sum(grid .> 1)
end

get_ans(example)
@assert get_ans(example) == 5
get_ans(load_input(5))

# Part 2

function populate_all!(grid::Matrix{Int64}, lines::AbstractVector{Line})
    for line in lines
        if horizontal(line)
            i_y = line.start_y + 1
            @inbounds for i_x in range_index_x(line)
                grid[i_x, i_y] += 1
            end
        elseif vertical(line)
            i_x = line.start_x + 1
            @inbounds for i_y in range_index_y(line)
                grid[i_x, i_y] += 1
            end
        else
            # It should be the case that the line as at 45 degrees...
            Ix = range_index_x(line)
            Iy = range_index_y(line)
            @assert length(Ix) == length(Iy)
            for (i_x, i_y) in zip(Ix, Iy)
                grid[i_x, i_y] += 1
            end
        end
    end
    return grid
end

populate_all(lines) = populate_all!(make_grid(lines), lines)

function get_ans2(input::AbstractString)
    lines = parse(Vector{Line}, input)
    grid = populate_all(lines)
    return sum(grid .> 1)
end

@assert get_ans2(example) == 12
get_ans2(load_input(5))


text:
using NPZ
using CSV
using DataFrames
using Random
using FileIO
using HDF5

# Parameter setup
csv_path = "C:/Users/Al-Najar/PycharmProjects/bathymetry_estimation/models/v8/data/v8_raw_40x41.csv"
dataset_size = 10
random_seed = 1
Random.seed!(random_seed)

# Dataset preparation
dataset_csv = CSV.File(csv_path; header=["col_path", "col_depth", "col_hmo", "col_freqpeak", "col_thetapeak"])
df = DataFrame(dataset_csv)
# df = df[shuffle(1:size(df, 1)),:]  # Shuffling
# df = df[1:dataset_size, [:col_path, :col_depth]]  # Reducing

for i in 1:size(df, 1)
    path = df.col_path[i]
    println(string(i, " -- ", size(df)))
    npyfile = npzread(string(path, ".npy"))
    h5open(string(path, ".h5"), "w") do file
        write(file, "data", npyfile)  # alternatively, say "@write file A"
    end
end


text:
"""
	MCprobability!(state::QState; kwargs...)
Get probability distribution of qubit configurations by Monte Carlo sampling.

# Argument
- `state`: Any QState with **measure!** defined
- `samplesize`: Keyword argument to specify sample size. By default `samplesize = 1024`.
- `percentage`: Keyword argument. By default `percentage = false`. Return counts of configurations if set to `false`; return percentage frequency if set to `true`.

# Return Value
Returns a dinctionary data structure. The key of each pair is the qubit configuration; the value of the pair is frequency of appearance in sampling.

# Example
``` julia
julia> state = MPSState(2)
2-qubit MPSState

julia> apply!(state, H, ActPosition(1))
2-qubit MPSState

julia> d =MCprobability!(state; samplesize=1000)
Dict{Int64,Real} with 2 entries:
  0 => 509
  1 => 491

julia> d =MCprobability!(state; samplesize=1000, percentage = true)
Dict{Int64,Real} with 2 entries:
  0 => 0.47168
  1 => 0.504883
```
"""
function MCprobability!(state::QState; kwargs...)
	samplesize = get(kwargs, :samplesize, 1024)
	result = measure!(state, [1:length(state);], samplesize; binary = false)
	return frequency(result; kwargs...)
end

function frequency(v::Vector; kwargs...)
	percentage = get(kwargs, :percentage, false)
 	samplesize = get(kwargs, :samplesize, length(v))
 	cutoff = get(kwargs, :cutoff, 0)
 	abscutoff = get(kwargs, :abscutoff, 0)
 	if abscutoff !=0 && cutoff !=0
 		error("cannot specify both cutoff and abscutoff")
 	end
 	abscutoff ==0 && (abscutoff = cutoff*samplesize)
	dict = Dict{Int, Real}()
	for key ∈ v
		push!(dict, key => get!(dict, key, 0)+1)
	end
	for (k,s) ∈ dict
		if s <= abscutoff
			delete!(dict, k)
		elseif percentage
			push!(dict, k => s/samplesize)
		end
	end

	return dict
end


text:
@testset "The 2-Sphere Sphere(2)" begin
x = SnPoint([1.,0.,0.])
y = SnPoint([0.,1.,0.])
z = SnPoint([0.,0.,1.])
M = Sphere(2)
ξ = log(M,x,y)
y2 = exp(M,x,ξ)
ν = log(M,x,z)
# Test unary operator
@test getValue( (-ν) ) == getValue(-ν)
@test opposite(M,x) == SnPoint(-getValue(x))
@test norm( getValue(y) - getValue(y2) ) ≈ 0 atol=10.0^(-16)
@test norm(M,x,ξ) ≈ norm( getValue(ξ) ) atol=10.0^(-16)
@test distance(M,x,y) ≈ norm(M,x,ξ) atol=10.0^(-16)
@test dot(M,x,ξ,ν) ≈ 0 atol = 10.0^(-16)
@test norm( getValue( mean(M,[x,y,z]) ) - 1/sqrt(3)*ones(3)) ≈ 0 atol=10.0^(-7)
@test norm( getValue( mean(M,[x,y]) ) - [1/sqrt(2),1/sqrt(2),0] ) ≈ 0 atol=10.0^(-15)
@test manifoldDimension(M) == manifoldDimension(x)
# Test randoms
@test validateMPoint(M,randomMPoint(M))
@test validateTVector(M,x,randomTVector(M,x))
# Test extended
xT = MPointE(x); yT = MPointE(y); zT = MPointE(z);
@test_throws DomainError dot(M,xT,log(M,xT,zT),log(M,yT,zT) )
@test dot(M,x,log(M,x,z),log(M,x,y) ) ≈ 0 atol=10.0^(-15)
# Tst ONB
@test tangentONB(M,x,y) == ( [SnTVector([0.,1.,0.]), SnTVector([0., 0., 1.])], [0.,1.] )
#
@test_throws ErrorException validateMPoint(M,SnPoint([1.,0.]))
@test_throws ErrorException validateMPoint(M,SnPoint([2.,0.,0.]))
@test_throws ErrorException validateTVector(M,x,SnTVector([1., 0.]))
@test_throws ErrorException validateTVector(M,x,log(M,z,x))

#check that PT(q->p, log_qp) = -log_pq (unitary minus already checked)
@test parallelTransport(M,y,x,log(M,y,x)) == -ξ
@test typeofTVector(typeof(x)) == SnTVector{Float64}
@test typeofMPoint(typeof(ξ)) == SnPoint{Float64}
@test typeofTVector(x) == SnTVector{Float64}
@test typeofMPoint(ξ) == SnPoint{Float64}
@test typicalDistance(M) == π
@test "$M" == "The 2-Sphere"
@test "$x" == "Sn([1.0, 0.0, 0.0])"
@test "$ξ" == "SnT([0.0, 1.5708, 0.0])"
#
#
#
@test_throws ErrorException x⊗y
@test_throws ErrorException x⊗y

end