text:
using Test
using StaticArrays, IntervalSets, LinearAlgebra

import ClimaCore:
    ClimaCore,
    slab,
    Spaces,
    Domains,
    Meshes,
    Geometry,
    Topologies,
    Spaces,
    Fields,
    Operators,
    level
import ClimaCore.Domains.Geometry: ⊗
import ClimaCore.Utilities: half

function hvspace_2D(;
    xlim = (-π, π),
    zlim = (0, 4π),
    helem = 10,
    velem = 64,
    npoly = 7,
)
    FT = Float64
    vertdomain = Domains.IntervalDomain(
        Geometry.ZPoint{FT}(zlim[1]),
        Geometry.ZPoint{FT}(zlim[2]);
        boundary_tags = (:bottom, :top),
    )
    vertmesh = Meshes.IntervalMesh(vertdomain, nelems = velem)
    vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)

    horzdomain = Domains.IntervalDomain(
        Geometry.XPoint{FT}(xlim[1]) .. Geometry.XPoint{FT}(xlim[2]),
        periodic = true,
    )
    horzmesh = Meshes.IntervalMesh(horzdomain, nelems = helem)
    horztopology = Topologies.IntervalTopology(horzmesh)

    quad = Spaces.Quadratures.GLL{npoly + 1}()
    horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)

    hv_center_space =
        Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)
    hv_face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(hv_center_space)
    return (hv_center_space, hv_face_space)
end

# V_face = Covariant3Vector
# V_center = Covariant12Vector
# V = V_center + V_face

# divergence(V) = horz_div(V) + vert_div(V)
# divergence(V) = horz_div(V_center) + horz_div(V_face) + vert_div(V_center) + vert_div(V_face)

# 1) horz_div(V_center): project to Contravariant1 + Contravariant2, take spectral derivative
# 2) horz_div(V_face): project to Contravariant1 + Contravariant2, interpolate to center, take spectral derivative
#   - will be zero if orthogional geom
# 3) vert_div(V_center): project to Contravariant3, interpolate to face, take FD deriv
#   - will be zero if orthogional geom
# 4) vert_div(V_face): project to Contravariant3, take FD deriv

@testset "1D SE, 1D FD Extruded Domain level extraction" begin
    hv_center_space, hv_face_space =
        hvspace_2D(npoly = 4, velem = 10, helem = 10)

    fcoord = Fields.coordinate_field(hv_face_space)
    ccoord = Fields.coordinate_field(hv_center_space)
    @test parent(Fields.field_values(level(fcoord.x, half))) == parent(
        Fields.field_values(
            Fields.coordinate_field(hv_face_space.horizontal_space).x,
        ),
    )
    @test parent(Fields.field_values(level(ccoord.x, 1))) == parent(
        Fields.field_values(
            Fields.coordinate_field(hv_center_space.horizontal_space).x,
        ),
    )
    @test parent(Fields.field_values(level(fcoord.z, half))) ==
          parent(
        Fields.field_values(
            Fields.coordinate_field(hv_face_space.horizontal_space).x,
        ),
    ) .* 0
end

@testset "1D SE, 1D FD Extruded Domain matrix interpolation" begin
    hv_center_space, hv_face_space =
        hvspace_2D(npoly = 4, velem = 10, helem = 10)

    center_field = sin.(Fields.coordinate_field(hv_center_space).z)
    face_field = sin.(Fields.coordinate_field(hv_face_space).z)

    for npoints in (3, 8)
        M_center = Operators.matrix_interpolate(center_field, npoints)
        M_face = Operators.matrix_interpolate(face_field, npoints)
        @test size(M_center) == (10, 10 * npoints)
        @test size(M_face) == (10 + 1, 10 * npoints)
    end
end

@testset "1D SE, 1D FD Extruded Domain vertical advection operator" begin

    # Advection Operator
    # c ∂_z f
    # for this test, we use f(z) = sin(z) and c = 1, a WVector
    # => c ∂_z f = cos(z)
    hv_center_space, hv_face_space = hvspace_2D()

    function advection(c, f)
        adv = zeros(eltype(f), hv_center_space)
        A = Operators.AdvectionC2C(
            bottom = Operators.SetValue(0.0),
            top = Operators.Extrapolate(),
        )
        return @. adv = A(c, f)
    end

    # advective velocity
    c = Geometry.WVector.(ones(Float64, hv_face_space),)
    # scalar-valued field to be advected
    f = sin.(Fields.coordinate_field(hv_center_space).z)

    # Call the advection operator
    adv = advection(c, f)

    @test norm(adv .- cos.(Fields.coordinate_field(hv_center_space).z)) ≤ 5e-2
end

@testset "1D SE, 1D FD Extruded Domain horizontal divergence operator" begin

    # Divergence Operator
    # ∂_x (c f)
    # for this test, we use f(x) = sin(x) and c = 2, a UVector
    # => ∂_x (c f) = 2 * cos(x)

    function divergence(c, f)
        divf = zeros(eltype(f), hv_center_space)
        # horizontal divergence operator applied to all levels
        hdiv = Operators.Divergence()
        divf .= hdiv.(f .* Ref(c)) # Ref is needed to treat c as a scalar for the broadcasting operator
        Spaces.weighted_dss!(divf)
        return divf
    end

    hv_center_space, _ = hvspace_2D()
    f = sin.(Fields.coordinate_field(hv_center_space).x)

    # Set up constant velocity field
    c = Geometry.UVector(2.0)

    # Call the divergence operator
    divf = divergence(c, f)

    @test norm(
        divf .- 2.0 .* cos.(Fields.coordinate_field(hv_center_space).x),
    ) ≤ 5e-5
end

@testset "1D SE, 1D FD Extruded Domain horz & vert divergence operator" begin

    # Divergence operator in 2D Cartesian domain with
    # Cₕ = (cₕ, 0), Cᵥ = (0, cᵥ)
    # ∇ₕ⋅(Cₕ * f) + ∇ᵥ⋅(Cᵥ * f)
    # here cₕ == cᵥ == 1

    # NOTE: the equation setup is only correct for Cartesian domains!

    hv_center_space, hv_face_space = hvspace_2D(xlim = (-1, 1), zlim = (-1, 1))

    function divergence(f)

        divuf = Fields.FieldVector(h = zeros(eltype(f), hv_center_space))
        h = f.h
        dh = divuf.h

        # vertical advection no inflow at bottom
        # and outflow at top
        Ic2f = Operators.InterpolateC2F(
            top = Operators.Extrapolate(),
            bottom = Operators.Extrapolate(),
        )
        divf2c = Operators.DivergenceF2C()
        # only upward component of divergence
        @. dh = divf2c(Ic2f(h) * Geometry.WVector(1.0))

        # only horizontal component of divergence
        hdiv = Operators.Divergence()
        @. dh += hdiv(h * Geometry.UVector(1.0)) # add the two components for full divergence +=
        Spaces.weighted_dss!(dh)

        return divuf
    end

    # A horizontal Gaussian blob, centered at (-x_0, -z_0),
    # with `a` the height of the blob and σ the standard deviation
    function h_blob(x_0, z_0, a = 1.0, σ = 0.2)
        coords = Fields.coordinate_field(hv_center_space)
        f = map(coords) do coord
            a * exp(-((coord.x + x_0)^2 + (coord.z + z_0)^2) / (2 * σ^2))
        end

        return f
    end

    # Spatial derivative of a horizontal Gaussian blob, centered at (-x_0, -z_0),
    # with `a` the height of the blob and σ the standard deviation
    function div_h_blob(x_0, z_0, a = 1.0, σ = 0.2)
        coords = Fields.coordinate_field(hv_center_space)
        div_uf = map(coords) do coord
            -a * (exp(-((coord.x + x_0)^2 + (coord.z + z_0)^2) / (2 * σ^2))) /
            (σ^2) * ((coord.x + x_0) + (coord.z + z_0))
        end

        return div_uf
    end

    f = Fields.FieldVector(h = h_blob(0.5, 0.5))

    divuf = divergence(f)
    exact_divuf = div_h_blob(0.5, 0.5)

    @test norm(exact_divuf .- divuf.h) ≤ 5e-2
end

@testset "1D SE, 1D FD Extruded Domain vertical vector advection" begin

    # Vector Advection Operator
    # C ∂_z F
    # for this test, we use F(z) = sin(z),
    # with F a UVector and C = ones, a WVector constant field
    # => C ∂_z F = cos(z)
    hv_center_space, hv_face_space = hvspace_2D()

    function advect(f)
        advf = zeros(eltype(f), hv_center_space)
        A = Operators.AdvectionC2C(
            bottom = Operators.SetValue(Geometry.UVector(0.0)), # value of f at the boundary (UVector field)
            top = Operators.Extrapolate(),
        )
        @. advf = A(vC, f)
    end

    # Vertical advective velocity
    vC = Geometry.WVector.(ones(Float64, hv_face_space),)
    # vector-valued field to be advected (one component only, a UVector)
    f = Geometry.UVector.(sin.(Fields.coordinate_field(hv_center_space).z),)

    advf = advect(f)

    function div!(F)
        vecdivf = zeros(eltype(F), hv_center_space)
        Ic2f = Operators.InterpolateC2F()
        divf2c = Operators.DivergenceF2C(
            bottom = Operators.SetValue(
                Geometry.WVector(1.0) ⊗ Geometry.UVector(0.0),
            ),
            top = Operators.Extrapolate(),
        )
        # only upward advection
        @. vecdivf = divf2c(vC ⊗ Ic2f(F))

        hdiv = Operators.Divergence()
        @. vecdivf += hdiv(Geometry.UVector(1.0) ⊗ F)
        Spaces.weighted_dss!(vecdivf)
        return vecdivf
    end

    F = Geometry.UVector.(sin.(Fields.coordinate_field(hv_center_space).z),)

    # Call the divergence operator
    vecdivf = div!(F)

    @test norm(
        vecdivf .-
        Geometry.UVector.(cos.(Fields.coordinate_field(hv_center_space).z),),
    ) ≤ 6.5e-2

end

@testset "1D SE, 1D FD Extruded Domain scalar diffusion" begin

    # Scalar diffusion operator in 2D
    # ∂_xx u + ∂_zz u

    hv_center_space, hv_face_space = hvspace_2D(xlim = (-5, 5), zlim = (-5, 5))

    K = 1.0

    function diff(K, u)

        diffu = zeros(eltype(u), hv_center_space)
        gradc2f = Operators.GradientC2F(
            top = Operators.SetValue(0.0),
            bottom = Operators.SetValue(0.0),
        )
        divf2c = Operators.DivergenceF2C()
        @. diffu = divf2c(K * gradc2f(u))

        hgrad = Operators.Gradient()
        hdiv = Operators.Divergence()

        @. diffu += hdiv(K * hgrad(u))
        Spaces.weighted_dss!(diffu)
        return diffu
    end

    # scalar-valued field to be diffused, u = exp(-(coord.x^2 + coord.z^2) / 2)
    u = map(Fields.coordinate_field(hv_center_space)) do coord
        exp(-(coord.x^2 + coord.z^2) / 2)
    end

    # Laplacian of u = exp(-(coord.x^2 + coord.z^2) / 2)
    function laplacian_u()
        coords = Fields.coordinate_field(hv_center_space)
        laplacian_u = map(coords) do coord
            ((coord.x^2 - 1) + (coord.z^2 - 1)) *
            (exp(-((coord.x^2) / 2 + (coord.z^2) / 2)))
        end

        return laplacian_u
    end

    # Call the diffusion operator
    diffu = diff(K, u)

    exact_laplacian = laplacian_u()

    @test norm(exact_laplacian .- diffu) ≤ 1e-2

end

@testset "1D SE, 1D FD Extruded Domain vector diffusion" begin

    # Vector diffusion operator in 2D, of vector-valued field U, with U == u
    # (one component only)
    # ∂_xx u + ∂_zz u

    hv_center_space, hv_face_space = hvspace_2D(xlim = (-5, 5), zlim = (-5, 5))

    K = 1.0

    function vec_diff(K, U)

        vec_diff = zeros(eltype(U), hv_center_space)
        gradc2f = Operators.GradientC2F(
            top = Operators.SetValue(Geometry.UVector(0.0)),
            bottom = Operators.SetValue(Geometry.UVector(0.0)),
        )
        divf2c = Operators.DivergenceF2C()
        @. vec_diff = divf2c(K * gradc2f(U))

        hgrad = Operators.Gradient()
        hdiv = Operators.Divergence()

        @. vec_diff += hdiv(K * hgrad(U))
        Spaces.weighted_dss!(vec_diff)

        return vec_diff
    end

    # vector-valued field to be diffused (one component only) u = exp(-(coord.x^2 + coord.z^2) / 2)
    U = map(Fields.coordinate_field(hv_center_space)) do coord
        Geometry.UVector(exp(-(coord.x^2 + coord.z^2) / 2))
    end

    # Laplacian of u = exp(-(coord.x^2 + coord.z^2) / 2)
    function hessian_u()
        coords = Fields.coordinate_field(hv_center_space)
        hessian_u = map(coords) do coord
            Geometry.UVector(
                ((coord.x^2 - 1) + (coord.z^2 - 1)) *
                (exp(-((coord.x^2) / 2 + (coord.z^2) / 2))),
            )
        end

        return hessian_u
    end

    # Call the vector diffusion operator
    vec_diff = vec_diff(K, U)

    exact_hessian = hessian_u()

    @test norm(exact_hessian .- vec_diff) ≤ 1e-2

end

@testset "curl" begin
    hv_center_space, hv_face_space =
        hvspace_2D(xlim = (-pi, pi), zlim = (-pi, pi))

    ccoords = Fields.coordinate_field(hv_center_space)
    fcoords = Fields.coordinate_field(hv_face_space)

    u =
        Geometry.transform.(
            Ref(Geometry.Covariant1Axis()),
            Geometry.UVector.(sin.(ccoords.x .+ 2 .* ccoords.z)),
        )
    w =
        Geometry.transform.(
            Ref(Geometry.Covariant3Axis()),
            Geometry.WVector.(cos.(3 .* fcoords.x .+ 4 .* fcoords.z)),
        )

    curl = Operators.Curl()
    curlC2F = Operators.CurlC2F(
        bottom = Operators.SetValue(Geometry.Covariant1Vector(0.0)),
        top = Operators.SetValue(Geometry.Covariant1Vector(0.0)),
    )

    curlu = curlC2F.(u)
    curlw = curl.(w)
    curluw = curlu .+ curlw
    Spaces.weighted_dss!(curluw)

    curlw_ref =
        Geometry.Contravariant2Vector.(
            3 .* sin.(3 .* fcoords.x .+ 4 .* fcoords.z),
        )
    curlu_ref =
        Geometry.Contravariant2Vector.(2 .* cos.(fcoords.x .+ 2 .* fcoords.z))

    curluw_ref = curlu_ref .+ curlw_ref

    # TODO: make SetValue do reasonable things on Extruded spaces (#79)
    zeroboundary = Operators.SetBoundaryOperator(
        bottom = Operators.SetValue(Geometry.Contravariant2Vector(0.0)),
        top = Operators.SetValue(Geometry.Contravariant2Vector(0.0)),
    )
    curluw_diff = zeroboundary.(curluw .- curluw_ref)

    @test norm(curluw_diff) < norm(curluw_ref) * 1e-2
end

@testset "curl-cross" begin
    hv_center_space, hv_face_space =
        hvspace_2D(xlim = (-pi, pi), zlim = (-pi, pi))

    ccoords = Fields.coordinate_field(hv_center_space)
    fcoords = Fields.coordinate_field(hv_face_space)

    a = 1.0
    b = 2.0
    c = 3.0
    d = 4.0

    u =
        Geometry.transform.(
            Ref(Geometry.Covariant1Axis()),
            Geometry.UVector.(a .* ccoords.x .+ b .* ccoords.z),
        )
    w =
        Geometry.transform.(
            Ref(Geometry.Covariant3Axis()),
            Geometry.WVector.(c .* fcoords.x .+ d .* fcoords.z),
        )

    curl = Operators.Curl()
    curlC2F = Operators.CurlC2F(
        bottom = Operators.SetCurl(Geometry.Contravariant2Vector(b)),
        top = Operators.SetCurl(Geometry.Contravariant2Vector(b)),
    )

    curlu = curlC2F.(u)
    curlw = curl.(w)
    curluw = curlu .+ curlw
    @test norm(curluw .- Ref(Geometry.Contravariant2Vector(b - c))) < 1e-10

    Ic2f = Operators.InterpolateC2F(
        bottom = Operators.Extrapolate(),
        top = Operators.Extrapolate(),
    )

    fu¹² = Geometry.Contravariant1Vector.(Geometry.Covariant13Vector.(Ic2f.(u))) # Contravariant12Vector in 3D
    fu³ = Geometry.Contravariant3Vector.(Geometry.Covariant13Vector.(w))

    Geometry.WVector.(Geometry.Covariant13Vector.(curluw .× fu¹²))

    curl_cross =
        Geometry.UWVector.(
            Geometry.Covariant13Vector.(curluw .× fu¹²) .+
            Geometry.Covariant13Vector.(curluw .× fu³),
        )

    curl_cross_ref =
        Geometry.UWVector.(
            (b - c) .* (c .* fcoords.x .+ d .* fcoords.z),
            .-(b - c) .* (a .* fcoords.x .+ b .* fcoords.z),
        )

    @test curl_cross ≈ curl_cross_ref rtol = 1e-2

end

@testset "2D hybrid hyperdiffusion" begin
    hv_center_space, hv_face_space =
        hvspace_2D(xlim = (-pi, pi), zlim = (-pi, pi))

    coords = Fields.coordinate_field(hv_center_space)
    k = 2
    y = @. sin(k * coords.x)
    ∇⁴y_ref = @. k^4 * sin(k * coords.x)

    wdiv = Operators.WeakDivergence()
    grad = Operators.Gradient()
    χ = Spaces.weighted_dss!(@. wdiv(grad(y)))
    ∇⁴y = Spaces.weighted_dss!(@. wdiv(grad(χ)))

    @test ∇⁴y_ref ≈ ∇⁴y rtol = 2e-2
end


text:
module EquivariantOperators

include("operators.jl")
export Grid, get,put!, fieldnorm
export FieldProd, fieldprod, FieldConv, fieldconv,⊗,⋅,⨉
export Op,Radfunc,remake!,nae
#
end # module

ext: jl


text: println("hello there")


text:
"""

Calculate the parent probabilites of each event.

    parent_probs(events, bg, kappa, kern)

# Arguments

* `events` Array of event times.
* `bg` Background intensity
* `kappa` Child excitation parameter.
* `kern` Kernel decay function

# Examples

```julia
kern_f(x) = pdf.(Distributions.Exponential(1/0.5), x)
parent_probs([1,2,3], 0.5, 0.5, )
```

"""
function parent_probs(events::Array{<:Number, 1}, bg, kappa::Float64, kern::Function)
    eventdifferencelist = event_difference_list(events)
    parent_probs(events, bg, kappa, kern, eventdifferencelist)
end

function parent_probs(events::Array{<:Number, 1}, bg, kappa::Float64, kern::Function, eventdifferencelist)
    parentprobs = Array{Array{Float64, 1}}(undef, length(events))
    for i in 1:length(events)
        parentprobs[i] = event_parent_prob(events[i], bg, kappa, kern, eventdifferencelist[i])
    end
    return parentprobs
end

function event_parent_prob(event::Number, bg::Number, kappa::Float64, kern::Function, eventdifference::Array{<:Number, 1})
#    probs = Array{Float64, 1}(undef, length(eventdifference) + 1)
#    probs[1] = bg
#    probs[2:end] .= kappa * kern.(eventdifference)
#    probs
    n = length(eventdifference)
    probs = Array{Float64, 1}(undef, n + 1)
    probs[1] = bg
    for i = 2:(n+1)
        probs[i] = kappa * kern(eventdifference[i-1])
    end
    probs
end

function event_parent_prob(event::Number, bg::Function, kappa::Float64, kern::Function, eventdifference::Array{<:Number, 1})
    event_parent_prob(event, bg(event), kappa, kern, eventdifference)
end


text:
mldatasets = ["MNIST", "FashionMNIST", "CIFAR10", "SVHN2"]

"""
	load_mldatasets_data(dataset::String[; anomaly_class_ind=1])

Loads MNIST, FMNIST, SVHN2 and CIFAR10 datasets.
"""
function load_mldatasets_data(dataset::String; anomaly_class_ind::Int=1)
	(dataset in mldatasets) ? nothing : error("$dataset not available in MLDatasets.jl")
	# since the functions for MLDatasets.MNIST, MLDatasets.CIFAR10 are the same
	sublib = getfield(MLDatasets, Symbol(dataset))

	# do we need to download it?
	isdir(joinpath(first(DataDeps.standard_loadpath), dataset)) ?
		nothing : sublib.download()

	# now get the data
	tr_x, tst_x = sublib.traintensor(Float32), sublib.testtensor(Float32)
	if ndims(tr_x) == 3 # some datasets are 3D tensors :(
		tr_x = reshape(tr_x, size(tr_x,1), size(tr_x,2), 1, :)
		tst_x = reshape(tst_x, size(tst_x,1), size(tst_x,2), 1, :)
	end
	data = cat(tr_x, tst_x, dims=4)
	labels = cat(sublib.trainlabels(), sublib.testlabels(), dims=1)

	# return the normal and anomalous data
	label_list = unique(labels)
	aclass_inds = labels .== label_list[anomaly_class_ind] # indices of anomalous data
	data[:,:,:,.!aclass_inds], data[:,:,:,aclass_inds]
end


text:
include("LinearSolver.jl")

abstract type AbstractBorderedLinearSolver <: AbstractLinearSolver end

# call for using BorderedArray input, specific to Arclength Continuation
(lbs::AbstractBorderedLinearSolver)(J, dR::vec1, dz::vec2, R, n::T, theta::T; shift::Ts = nothing) where {vec1, vec2, T, Ts} = (lbs)(J, dR, dz.u, dz.p, R, n, theta / length(dz.u), one(T) - theta; shift = shift)

####################################################################################################
@with_kw struct BorderingBLS{S <: AbstractLinearSolver, Ttol} <: AbstractBorderedLinearSolver
	solver::S
	tol::Ttol = 1e-12
	checkPrecision::Bool = false
end

# dummy constructor to simplify user passing options to continuation
BorderingBLS() = BorderingBLS(solver = DefaultLS())
BorderingBLS(ls::AbstractLinearSolver) = BorderingBLS(solver = ls)

# solve in dX, dl
#          J  * dX +       dR   * dl = R
# xiu * dz.u' * dX + xip * dz.p * dl = n
function (lbs::BorderingBLS{S, Ttol})(  J, dR,
								dzu, dzp::T, R, n::T,
								xiu::T = T(1), xip::T = T(1); shift::Ts = nothing)  where {T, S, Ts, Ttol}
	# the following parameters are used for the pseudo arc length continuation
	# xiu = theta / length(dz.u)
	# xip = one(T) - theta

	# we make this branching to avoid applying a zero shift
	if isnothing(shift)
		x1, x2, _, (it1, it2) = lbs.solver(J, R, dR)
	else
		x1, x2, _, (it1, it2) = lbs.solver(J, R, dR; a₀ = shift)
	end

	dl = (n - dot(dzu, x1) * xiu) / (dzp * xip - dot(dzu, x2) * xiu)

	# dX = x1 .- dl .* x2
	axpy!(-dl, x2, x1)

	# we check the precision of the solution by the bordering algorithm
	# mainly for debugging purposes
	if lbs.checkPrecision
		# at this point, x2 is not used anymore, we can use it for computing the residual
		# hence x2 = J*x1 + dl*dR - R
		x2 = apply(J, x1)
		axpy!(dl, dR, x2)
		axpy!(-1, R, x2)

		printstyled(color=:red,"--> res = ", ( norm(x2), abs(n - xip*dzp*dl -xiu* dot(dzu, x1))), "\n")
	end

	return x1, dl, true, (it1, it2)
end
####################################################################################################
# this interface should work for Sparse Matrices as well as for Matrices
struct MatrixBLS{S <: AbstractLinearSolver} <: AbstractBorderedLinearSolver
	solver::S
end

# dummy constructor to simplify user passing options to continuation
MatrixBLS() = MatrixBLS(DefaultLS())

# case of a scalar additional linear equation
function (lbs::MatrixBLS)(J, dR,
						dzu, dzp::T, R::vectype, n::T,
						xiu::T = T(1), xip::T = T(1); shift::Ts = nothing)  where {T <: Number, vectype <: AbstractVector, S, Ts}
	N = length(dzu)

	rhs = vcat(R, n)
	# A = similar(J, N + 1, N + 1)
	# if shift == nothing
	# 	A[1:N, 1:N] .= J
	# else
	# 	A[1:N, 1:N] .= J + shift * I
	# end
	#
	# A[1:N, end] .= dR
	# A[end, 1:N] .= vec(dzu) .* xiu
	# A[end, end]  = dzp * xip

	# A[1:N, end] .= dR
	# A[end, 1:N] .= vec(dzu) .* xiu
	# A[end, end]  = dzp * xip

	if isnothing(shift)
		A = J
	else
		A = J + shift * I
	end
	A = hcat(A, dR)
	A = vcat(A, vcat(vec(dzu) .* xiu, dzp * xip)')

	res = A \ rhs
	return res[1:end-1], res[end], true, 1
end

# case of a scalar additional linear equation
function (lbs::MatrixBLS)(J, dR,
						dzu, dzp, R::vectype, n::Tv,
						xiu::T = T(1), xip::T = T(1); shift::Ts = nothing)  where {Tv, T <: Number, vectype <: AbstractVector, S, Ts}
	@warn "Experimental!! Solving bordered linear system with non scalar bordered equation"
	N = size(J, 1)
	n1, n2 = size(dzp)

	rhs = vcat(R, n)
	A = similar(J, N + n1, N + n2)

	if shift == nothing
		A[1:N, 1:N] .= J
	else
		A[1:N, 1:N] .= J + shift * I
	end

	A[1:N, end-n1+1:end] .= dR
	A[end-n1+1:end, 1:N] .= (dzu) .* xiu
	A[end-n1+1:end, end-n1+1:end]  .= dzp .* xip
	res = A \ rhs
	return res[1:end-n1], res[end-n1+1:end], true, 1
end
####################################################################################################
# composite type to save the bordered linear system with expression
# [ J	a]
# [b'	c]
# It then solved using Matrix Free algorithm applied to the full operator and not just J as for MatrixFreeBLS
#
struct MatrixFreeBLSmap{Tj, Ta, Tb, Tc}
	J::Tj
	a::Ta
	b::Tb
	c::Tc
end

function (lbmap::MatrixFreeBLSmap{Tj, Ta, Tb, Tc})(x::BorderedArray{Ta, Tc}) where {Tj, Ta, Tb, Tc <: Number}
	out = similar(x)
	copyto!(out.u, apply(lbmap.J, x.u) .+ lbmap.a .* x.p)
	out.p  =   dot(lbmap.b, x.u)  + lbmap.c  * x.p
	return out
end

struct MatrixFreeBLS{S} <: AbstractBorderedLinearSolver
	solver::S
end

MatrixFreeBLS() = MatrixFreeBLS(DefaultLS())

# dummy constructor to simplify user passing options to continuation
# We restrict to bordered systems where the added component is scalar
# For now, we restrict to KrylovKit iterative solver because we did not make BorderedArray a subtype of AbstractVector
function (lbs::MatrixFreeBLS{S})(J, 		dR,
								dzu, 	dzp::T, R, n::T,
								xiu::T = T(1), xip::T = T(1); shift::Ts = nothing) where {T <: Number, S, Ts}
	linearmap = MatrixFreeBLSmap(J, dR, dzu * xiu, dzp * xip)
	rhs = BorderedArray(copy(R), n)
	sol, cv, it = lbs.solver(linearmap, rhs)
	return sol.u, sol.p, cv, it
end
####################################################################################################
# Nested algorithm for solving the bordered linear system


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Mortar3D.jl/blob/master/LICENSE

using Base.Test
using Mortar3D: calculate_normals

@testset "calculate normals" begin
    X = Dict(1 => [0.0, 0.0, 0.0], 2 => [1.0, 0.0, 0.0], 3 => [0.0, 1.0, 0.0])
    elements = Dict(1 => [1, 2, 3])
    element_types = Dict(1 => :Tri3)
    normals = calculate_normals(elements, element_types, X)
    @test isapprox(normals[1], [0.0, 0.0, 1.0])
end


text:
# TODO: Reviewers: Is this a correct list of keys to verify?
POWER_MODELS_KEYS = [
    "baseMVA",
    "branch",
    "bus",
    "dcline",
    "gen",
    "load",
    "name",
    "per_unit",
    "shunt",
    "source_type",
    "source_version",
    "storage",
]

badfiles = Dict("case30.m" => PSY.InvalidValue)

@testset "Parse Matpower data files" begin
    files = [x for x in readdir(joinpath(MATPOWER_DIR)) if splitext(x)[2] == ".m"]
    if length(files) == 0
        @error "No test files in the folder"
    end

    for f in files
        @info "Parsing $f..."
        path = joinpath(MATPOWER_DIR, f)

        pm_dict = PowerSystems.parse_file(path)
        for key in POWER_MODELS_KEYS
            @test haskey(pm_dict, key)
        end
        @info "Successfully parsed $path to PowerModels dict"

        if f in keys(badfiles)
            @test_logs(
                (:error, r"cannot create Line"),
                match_mode = :any,
                @test_throws(badfiles[f], System(PowerSystems.PowerModelsData(pm_dict)))
            )
        else
            sys = System(PowerSystems.PowerModelsData(pm_dict))
            @info "Successfully parsed $path to System struct"
        end
    end
end

@testset "Parse PowerModel Matpower data files" begin
    files = [
        x for x in readdir(joinpath(DATA_DIR, "pm_data", "matpower")) if
        splitext(x)[2] == ".m"
    ]
    if length(files) == 0
        @error "No test files in the folder"
    end

    for f in files
        @info "Parsing $f..."
        path = joinpath(joinpath(DATA_DIR, "pm_data", "matpower"), f)

        pm_dict = PowerSystems.parse_file(path)
        for key in POWER_MODELS_KEYS
            @test haskey(pm_dict, key)
        end
        @info "Successfully parsed $path to PowerModels dict"

        if f in keys(badfiles)
            @test_logs(
                (:error, r"cannot create Line"),
                match_mode = :any,
                @test_throws(badfiles[f], System(PowerSystems.PowerModelsData(pm_dict)))
            )
        else
            sys = System(PowerSystems.PowerModelsData(pm_dict))
            @info "Successfully parsed $path to System struct"
        end
    end
end


text:
# # Anyons

# We solve the almost-bosonic anyon model of https://arxiv.org/pdf/1901.10739.pdf

using DFTK
using StaticArrays
using Plots

# Unit cell. Having one of the lattice vectors as zero means a 2D system
a = 14
lattice = a .* [[1 0 0.]; [0 1 0]; [0 0 0]];

# Confining scalar potential
pot(x, y, z) = ((x - a/2)^2 + (y - a/2)^2)

# Parameters
Ecut = 50
n_electrons = 1
β = 5

terms = [Kinetic(2),
         ExternalFromReal(X -> pot(X...)),
         Anyonic(1, β)
]
model = Model(lattice; n_electrons=n_electrons,
              terms=terms, spin_polarization=:spinless)  # "spinless electrons"
basis = PlaneWaveBasis(model; Ecut, kgrid=(1, 1, 1), supersampling=1.2)
scfres = direct_minimization(basis, tol=1e-14)  # Reduce tol for production
E = scfres.energies.total
s = 2
E11 = π/2 * (2(s+1)/s)^((s+2)/s) * (s/(s+2))^(2(s+1)/s) * E^((s+2)/s) / β
println("e(1,1) / (2π)= ", E11 / (2π))
display(heatmap(scfres.ρ[:, :, 1, 1], c=:blues))


text:
## Definition of a phase-type distribution, and related functions.
using Distributions
using Statistics: cor, cov, median, std, quantile

import Base: minimum, maximum
import Distributions: cdf, insupport, logpdf, pdf
import Statistics: mean

macro check_args(D, cond)
    quote
        if !($(esc(cond)))
            throw(ArgumentError(string(
                $(string(D)), ": the condition ", $(string(cond)),
                " is not satisfied.")))
        end
    end
end

struct PhaseType <: ContinuousUnivariateDistribution
    # Defining properties
    π # entry probabilities
    T # transition probabilities

    # Derived properties
    t # exit probabilities
    p # number of phases

    function PhaseType(π, T, t, p)
        @check_args(PhaseType, all(π .>= zero(π[1])))
        @check_args(PhaseType, isapprox(sum(π), 1.0, atol=1e-4))
        @check_args(PhaseType, p == length(π) && p == length(t) &&
            all(p .== size(T)))
        @check_args(PhaseType, all(t .>= zero(t[1])))
        @check_args(PhaseType, all(isapprox(t, -T*ones(p))))
        new(π, T, t, p)
    end
end

PhaseType(π, T) = PhaseType(π, T, -T*ones(length(π)), length(π))
PhaseType(π, T, t) = PhaseType(π, T, t, length(π))

minimum(d::PhaseType) = 0
maximum(d::PhaseType) = Inf
insupport(d::PhaseType, x::Real) = x > 0 && x < Inf

pdf(d::PhaseType, x::Real) = transpose(d.π) * exp(d.T * x) * d.t
logpdf(d::PhaseType, x::Real) = log(pdf(d, x))
cdf(d::PhaseType, x::Real) = 1 - transpose(d.π) * exp(d.T * x) * ones(d.p)

mean(d::PhaseType) = -transpose(d.π) * inv(d.T) * ones(d.p)

iscoxian(d::PhaseType) = all(isapprox.(d.T[diagm(0 => ones(d.p),
    1 => ones(d.p-1)) .< 1], 0))

function reverse_coxian(π::AbstractArray{Float64}, T::AbstractArray{Float64},
        t::AbstractArray{Float64})
    # Get the time-reversed distribution
    Δ_v = Diagonal((-π' * inv(T))')
    πStar = (t' * Δ_v)'
    TStar = inv(Δ_v) * T' * Δ_v

    # Relabel the states from 1,..,p as p,...,1
    πCanon = reverse(πStar)
    λ = -reverse(diag(TStar))
    TCanon = diagm(0 => -λ, 1 => λ[1:end-1])
    (πCanon, TCanon, -TCanon*ones(length(π)))
end

function reverse_coxian(dist::PhaseType)
    PhaseType(reverse_coxian(dist.π, dist.T, dist.t) ...)
end

function sort_into_canonical_form(π::AbstractArray{Float64},
        T::AbstractArray{Float64}, t::AbstractArray{Float64})
    p = length(π)
    λ = -diag(T)
    π = π[:]

    while true
        swapped = false
        for i = 1:(p-1)
            if λ[i] > λ[i+1]
                w = λ[i+1]/λ[i]
                π[i] += (1-w)*π[i+1]
                π[i+1] *= w
                λ[i:i+1] = λ[i+1:-1:i]
                swapped = true
            end
        end

        if swapped == false
            break
        end
    end

    TSort = diagm(0 => -λ, 1 => λ[1:end-1])
    (π, TSort, -TSort*ones(p))
end

function sort_into_canonical_form(dist::PhaseType)
    PhaseType(sort_into_canonical_form(dist.π, dist.T, dist.t) ...)
end


text:
using ACCESS
using Test

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


text:
"""
    Ellipse{T} <: Surface{T}

Elliptical surface, not a valid CSG object.
The rotation of the rectangle around its normal is defined by `rotationvec`.
`rotationvec×surfacenormal` is taken as the vector along the u axis.

**Can be used as a detector in [`OpticalSystem`](@ref)s.**

```julia
Ellipse(halfsizeu::T, halfsizev::T, [surfacenormal::SVector{3,T}, centrepoint::SVector{3,T}]; interface::NullOrFresnel{T} = nullinterface(T))
```

The minimal case returns a ellipse centered at the origin with `surfacenormal = [0, 0, 1]`.
"""
struct Ellipse{T} <: Surface{T}
    plane::Plane{T,3}
    halfsizeu::T
    halfsizev::T
    uvec::SVector{3,T}
    vvec::SVector{3,T}

    function Ellipse(halfsizeu::T, halfsizev::T; interface::NullOrFresnel{T} = NullInterface(T)) where {T<:Real}
        @assert halfsizeu > zero(T) && halfsizev > zero(T)
        new{T}(Plane(SVector{3,T}(0, 0, 1), SVector{3,T}(0, 0, 0), interface = interface), halfsizeu, halfsizev, SVector{3,T}(1.0, 0.0, 0.0), SVector{3,T}(0.0, 1.0, 0.0))
    end

    function Ellipse(halfsizeu::T, halfsizev::T, surfacenormal::SVector{3,T}, centrepoint::SVector{3,T}; interface::NullOrFresnel{T} = NullInterface(T), rotationvec::SVector{3,T} = SVector{3,T}(0.0, 1.0, 0.0)) where {T<:Real}
        @assert halfsizeu > zero(T) && halfsizev > zero(T)
        n̂ = normalize(surfacenormal)
        if abs(dot(rotationvec, n̂)) == one(T)
            rotationvec = SVector{3,T}(1.0, 0.0, 0.0)
        end
        uvec = normalize(cross(normalize(rotationvec), n̂))
        vvec = normalize(cross(n̂, uvec))
        new{T}(Plane(n̂, centrepoint, interface = interface), halfsizeu, halfsizev, uvec, vvec)
    end

    function Ellipse(plane::Plane{T,3}, halfsizeu::T, halfsizev::T, uvec::SVector{3,T}, vvec::SVector{3,T}) where {T<:Real}
        new{T}(plane, halfsizeu, halfsizev, uvec, vvec)
    end
end
export Ellipse

Base.show(io::IO, a::Ellipse{T}) where {T<:Real} = print(io, "Ellipse{$T}($(centroid(a)), $(normal(a)), $(a.halfsizeu), $(a.halfsizev), $(interface(a)))")

interface(a::Ellipse{T}) where {T<:Real} = interface(a.plane)

uvrange(::Type{Ellipse{T}}) where {T<:Real} = ((-T(π), T(π)), (zero(T), one(T))) # θ and ρ

normal(r::Ellipse{T}) where {T<:Real} = normal(r.plane)
normal(r::Ellipse{T}, ::T, ::T) where {T<:Real} = normal(r)

point(r::Ellipse{T}, θ::T, ρ::T) where {T<:Real} = centroid(r) + ρ * (r.halfsizeu * cos(θ) * r.uvec + r.halfsizev * sin(θ) * r.vvec)
partials(r::Ellipse{T}, θ::T, ρ::T) where {T<:Real} = (ρ * (r.halfsizeu * -sin(θ) * r.uvec + r.halfsizev * cos(θ) * r.vvec), (r.halfsizeu * cos(θ) * r.uvec + r.halfsizev * sin(θ) * r.vvec))

uv(r::Ellipse{T}, x::T, y::T, z::T) where {T<:Real} = uv(r, SVector{3,T}(x, y, z))
function uv(r::Ellipse{T}, p::SVector{3,T}) where {T<:Real}
    v = dot(p - centroid(r), r.vvec)
    u = dot(p - centroid(r), r.uvec)
    θ = NaNsafeatan(v, u)
    rad = norm(r.halfsizeu * cos(θ) * r.uvec + r.halfsizev * sin(θ) * r.vvec)
    return SVector{2,T}(θ, norm(p - centroid(r)) / rad)
end

onsurface(a::Ellipse{T}, point::SVector{3,T}) where {T<:Real} = onsurface(a.plane, point) && zero(T) <= uv(a, point)[2] <= one(T)

function uvtopix(::Ellipse{T}, uv::SVector{2,T}, imsize::Tuple{Int,Int}) where {T<:Real}
    θ, ρ = uv
    h, w = imsize
    u = (cos(θ) * ρ + one(T)) / 2
    v = (sin(θ) * ρ + one(T)) / 2
    pixu = Int(floor((w - 1) * u)) + 1
    pixv = h - Int(floor((h - 1) * v))
    return pixu, pixv
end

centroid(r::Ellipse{T}) where {T<:Real} = r.plane.pointonplane

function surfaceintersection(ell::Ellipse{T}, r::AbstractRay{T,3}) where {T<:Real}
    interval = surfaceintersection(ell.plane, r)
    if interval isa EmptyInterval{T} || isinfiniteinterval(interval)
        return EmptyInterval(T) # no ray plane intersection or inside plane but no hit
    else
        intersect = halfspaceintersection(interval)
        p = point(intersect)
        θ, ρ = uv(ell, p)
        if ρ > one(T)
            return EmptyInterval(T) # no ray plane intersection
        else
            intuv = Intersection(α(intersect), p, normal(ell), θ, ρ, interface(ell))
            if dot(normal(ell), direction(r)) < zero(T)
                return positivehalfspace(intuv)
            else
                return rayorigininterval(intuv)
            end
        end
    end
end

function makemesh(c::Ellipse{T}, subdivisions::Int = 30) where {T<:Real}
    dθ = T(2π) / subdivisions
    centre = point(c, zero(T), zero(T))
    tris = Vector{Triangle{T}}(undef, subdivisions)
    for i in 0:(subdivisions - 1)
        θ1 = i * dθ - π
        θ2 = (i + 1) * dθ - π
        p1 = point(c, θ1, one(T))
        p2 = point(c, θ2, one(T))
        tris[i + 1] = Triangle(centre, p1, p2)
    end
    return TriangleMesh(tris)
end

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

function Circle(radius::T; interface::NullOrFresnel{T} = NullInterface(T)) where {T<:Real}
    return Ellipse(radius, radius; interface = interface)
end

"""
    Circle(radius, [surfacenormal, centrepoint]; interface = nullinterface(T))

Shortcut method to create a circle. The minimal case returns a circle centred at the origin with `normal = [0, 0, 1]`.
"""
function Circle(radius::T, surfacenormal::SVector{3,T}, centrepoint::SVector{3,T}; interface::NullOrFresnel{T} = NullInterface(T)) where {T<:Real}
    return Ellipse(radius, radius, surfacenormal, centrepoint; interface = interface, rotationvec = SVector{3,T}(0.0, 1.0, 0.0))
end

export Circle


text:
got_coords = [947 1023;
939 1017;
1017 1017;
916 1016;
948 1016;
973 1016;
931 1015;
975 1015;
1022 1015;
872 1014;
1023 1014;
875 1013;
869 1012;
864 1011;
903 1011;
926 1010;
969 1010;
973 1010;
1000 1010;
877 1009;
883 1009;
914 1009;
962 1009;
1005 1009;
976 1008;
1000 1008;
1001 1008;
891 1007;
901 1007;
902 1007;
999 1007;
1001 1007;
1014 1007;
923 1006;
975 1006;
863 1005;
982 1005;
872 1004;
914 1004;
985 1003;
926 1002;
932 1002;
904 1001;
919 1001;
989 1001;
879 1000;
986 1000;
919 999;
925 999;
901 998;
1009 998;
869 997;
1016 996;
1004 995;
999 994;
1018 993;
995 992;
1005 991;
829 990;
908 990;
939 990;
1029 990;
940 989;
984 988;
1000 986;
829 985;
908 985;
961 983;
900 982;
928 982;
1011 980;
1003 978;
934 976;
986 976;
1002 976;
1026 976;
983 975;
993 975;
1036 975;
1094 975;
750 974;
846 974;
929 974;
1020 974;
1031 973;
1037 973;
1041 973;
835 972;
868 972;
939 971;
967 971;
1124 971;
1126 971;
1127 971;
957 970;
1040 970;
855 969;
945 969;
957 969;
745 968;
1077 968;
1078 968;
1092 968;
909 967;
931 967;
946 967;
1130 967;
925 966;
1096 966;
1129 966;
741 965;
939 965;
1026 964;
929 963;
933 963;
1010 963;
905 962;
930 962;
970 962;
1041 962;
897 961;
934 961;
1018 961;
1035 961;
1038 961;
1112 961;
1118 961;
758 960;
766 960;
874 960;
877 960;
1034 960;
741 959;
1000 959;
1082 959;
734 958;
914 958;
1036 958;
1071 958;
1140 958;
1012 957;
1073 957;
1161 957;
935 956;
1104 956;
1107 956;
897 955;
1091 954;
1102 954;
718 953;
893 953;
911 953;
1015 953;
1019 953;
1132 953;
722 952;
886 952;
906 952;
1096 952;
1098 952;
1134 952;
757 951;
901 951;
1023 951;
1026 951;
771 950;
1022 950;
875 949;
902 949;
903 949;
917 949;
1063 949;
1169 949;
766 948;
759 947;
1096 947;
872 946;
898 946;
1122 946;
713 945;
1063 945;
784 944;
815 944;
924 944;
927 944;
1147 944;
1150 944;
1062 943;
706 942;
1011 942;
1022 942;
1115 942;
1136 942;
1116 941;
1016 940;
1029 940;
929 939;
1059 939;
761 938;
703 937;
743 937;
768 937;
880 937;
923 937;
940 937;
1144 937;
697 936;
760 936;
772 936;
780 936;
798 936;
936 936;
1049 936;
699 935;
702 935;
796 935;
1052 935;
770 934;
1029 933;
768 932;
807 932;
721 931;
768 931;
1038 931;
677 930;
722 930;
674 929;
687 929;
689 929;
695 929;
794 929;
815 929;
1016 929;
1052 929;
672 928;
707 928;
747 928;
799 928;
814 928;
1187 928;
1046 927;
1156 927;
691 926;
758 926;
873 926;
700 925;
795 925;
802 925;
815 924;
1038 924;
1141 924;
761 922;
794 922;
822 922;
1112 922;
1135 922;
663 921;
739 921;
1183 921;
1202 921;
1205 921;
695 920;
721 920;
750 920;
1104 920;
1179 920;
1200 920;
666 919;
686 919;
695 919;
1134 919;
1199 919;
748 918;
1161 918;
1193 918;
1199 918;
1051 917;
1060 917;
1112 917;
749 916;
1154 916;
1221 916;
671 915;
789 915;
830 915;
1053 915;
1064 915;
1100 915;
1133 915;
1160 915;
670 914;
738 914;
780 914;
775 913;
782 913;
799 913;
1110 913;
1138 913;
693 912;
826 912;
1055 912;
1183 912;
793 911;
799 911;
688 910;
1008 909;
658 908;
667 908;
1042 908;
1150 908;
1161 907;
1221 907;
649 906;
782 906;
831 906;
1125 906;
1159 906;
1181 906;
767 905;
783 905;
1212 905;
682 904;
785 904;
1123 904;
1174 904;
728 903;
1026 903;
1036 903;
1191 902;
779 901;
815 901;
1153 901;
1232 901;
716 900;
818 900;
843 900;
1043 900;
703 899;
804 899;
823 899;
1127 899;
733 898;
737 898;
792 898;
1191 898;
826 897;
1037 897;
1239 897;
696 896;
833 896;
1061 896;
665 895;
667 895;
672 895;
718 895;
1016 895;
1044 895;
1114 895;
1191 895;
1207 895;
1215 895;
646 894;
697 894;
827 894;
846 894;
1042 894;
1204 894;
1219 894;
680 893;
729 893;
810 893;
821 893;
1018 893;
1037 892;
1041 892;
1044 892;
1108 892;
1204 892;
660 891;
799 891;
808 891;
833 891;
1058 891;
1111 891;
1196 891;
1209 891;
1236 891;
667 889;
670 889;
843 889;
1009 889;
1118 889;
1217 889;
1223 889;
1251 889;
842 888;
636 887;
654 887;
849 887;
1086 887;
1244 887;
812 886;
1047 886;
713 885;
828 885;
839 885;
1122 885;
1257 885;
642 884;
687 884;
1026 884;
1244 884;
997 883;
1026 883;
1036 883;
1208 883;
990 882;
1193 882;
1031 881;
1223 881;
692 880;
826 880;
642 879;
1080 879;
1106 879;
1259 879;
1114 878;
1186 878;
1244 878;
1271 878;
1098 877;
1100 877;
1255 877;
635 876;
1089 876;
1242 876;
646 875;
997 875;
1018 875;
1200 875;
1259 875;
648 874;
1098 874;
1191 874;
819 873;
981 873;
996 873;
1257 873;
626 872;
723 872;
654 871;
724 871;
804 871;
988 871;
1246 871;
684 870;
813 870;
844 870;
1032 870;
632 869;
800 869;
830 869;
1022 869;
1031 869;
1276 869;
662 868;
823 868;
1081 868;
624 867;
1093 867;
1205 867;
1222 867;
991 866;
697 864;
699 864;
857 864;
1028 864;
1118 864;
1258 864;
977 863;
1195 863;
1264 863;
1277 863;
826 860;
1249 860;
1269 860;
823 859;
963 859;
986 859;
1217 859;
1302 859;
836 858;
844 858;
1009 858;
1093 858;
624 857;
683 857;
816 857;
1004 857;
1110 857;
1218 857;
1258 857;
1261 857;
1215 856;
1280 856;
1185 855;
1197 855;
668 854;
859 854;
954 854;
984 854;
1186 854;
1215 854;
602 853;
608 853;
626 853;
652 853;
983 853;
685 852;
828 851;
831 851;
846 851;
1299 851;
843 850;
849 850;
974 850;
1111 850;
1117 850;
1284 850;
1291 850;
614 849;
1191 849;
845 848;
968 848;
973 848;
1091 848;
1096 848;
1203 848;
701 847;
1270 847;
845 846;
945 846;
947 846;
603 845;
684 845;
1261 845;
1262 845;
849 844;
930 844;
818 843;
939 843;
971 843;
1182 843;
1245 843;
922 842;
984 842;
983 841;
828 840;
941 840;
1262 840;
1001 839;
1007 839;
1134 839;
685 838;
805 838;
1258 838;
1305 838;
617 837;
670 837;
805 837;
848 837;
1135 837;
1168 837;
1177 837;
1261 837;
821 836;
922 836;
618 835;
658 835;
1129 835;
1205 835;
1286 835;
1296 835;
1305 835;
825 834;
927 834;
995 834;
1183 834;
660 833;
817 833;
919 833;
986 833;
1112 833;
1153 833;
1177 833;
676 832;
976 832;
914 831;
920 831;
1096 831;
1099 831;
1174 831;
1194 831;
828 830;
905 830;
942 830;
962 830;
1161 830;
1280 830;
1303 830;
817 829;
963 829;
1188 829;
1281 828;
959 827;
1154 827;
1300 827;
957 826;
1284 826;
825 825;
855 825;
1110 825;
1113 824;
914 823;
1186 823;
854 822;
667 821;
1117 821;
599 820;
660 820;
831 820;
953 820;
924 819;
1169 819;
1318 819;
921 818;
1304 818;
1311 818;
902 817;
1325 817;
845 816;
824 815;
829 815;
821 814;
838 814;
895 814;
957 814;
985 814;
1113 814;
1141 814;
817 813;
916 813;
983 813;
842 812;
857 812;
889 812;
1282 812;
1158 811;
1163 811;
1310 811;
840 810;
845 810;
1172 810;
914 809;
1295 809;
839 808;
907 808;
942 808;
1317 808;
893 807;
922 807;
941 807;
1158 807;
1174 807;
822 806;
833 806;
893 806;
914 806;
938 806;
1148 806;
1183 806;
845 804;
889 804;
914 804;
944 804;
1299 804;
807 803;
828 803;
876 803;
892 803;
950 803;
1132 803;
1170 803;
803 802;
882 802;
1138 802;
1346 802;
819 801;
954 801;
1338 801;
1340 801;
827 800;
881 800;
933 800;
1335 800;
1333 799;
828 798;
918 798;
830 796;
1343 796;
805 795;
832 794;
885 794;
924 794;
1085 794;
1093 793;
1109 793;
816 792;
818 792;
888 792;
1038 792;
1106 792;
1345 792;
911 791;
889 790;
905 790;
1065 789;
1120 789;
1311 789;
850 788;
873 788;
912 788;
1129 788;
826 787;
908 787;
1014 787;
1048 787;
1092 787;
1108 787;
1342 787;
866 786;
874 786;
880 786;
891 786;
903 786;
1056 786;
1090 786;
1112 786;
1305 786;
1330 786;
832 785;
854 784;
912 784;
1023 784;
1028 784;
1043 784;
1317 784;
1330 784;
921 783;
1061 783;
1065 783;
1079 783;
1334 783;
860 782;
1101 782;
1118 782;
1135 782;
1138 782;
1151 782;
1158 782;
1336 782;
839 781;
897 781;
1041 781;
1047 781;
1316 781;
548 780;
844 780;
1040 780;
1061 780;
1082 780;
1095 780;
1140 780;
813 779;
990 779;
1042 779;
1058 779;
818 778;
843 778;
988 778;
1051 778;
1127 778;
1177 778;
1329 778;
878 777;
1117 777;
1120 777;
1310 777;
1338 777;
1085 776;
1100 776;
1176 776;
1306 776;
1320 776;
1010 775;
1072 775;
1094 775;
1104 775;
855 774;
915 774;
1003 774;
1344 774;
666 773;
886 773;
903 773;
1167 773;
1189 773;
913 772;
1055 772;
1073 772;
1096 772;
1118 772;
1138 772;
1178 772;
1185 772;
1338 772;
1350 772;
807 771;
883 771;
889 771;
989 771;
1016 771;
1063 771;
1089 771;
1143 771;
1145 771;
841 770;
1124 770;
914 769;
1106 769;
1166 769;
1193 769;
653 768;
807 768;
812 768;
990 768;
1144 768;
1164 768;
1172 768;
1198 768;
1200 768;
1359 768;
586 767;
1085 767;
1326 767;
621 766;
831 766;
914 766;
1033 766;
1202 766;
1208 766;
1344 766;
543 765;
548 765;
805 765;
975 765;
1037 765;
1091 765;
1102 765;
1162 765;
581 764;
875 764;
985 764;
1129 764;
1176 764;
579 763;
611 763;
623 763;
639 763;
659 763;
678 763;
821 763;
876 763;
1063 763;
1066 763;
690 762;
871 762;
1027 762;
1101 762;
1211 762;
610 761;
1095 761;
809 760;
887 760;
1007 760;
1148 760;
1175 760;
1176 760;
1233 760;
691 759;
699 759;
811 759;
814 759;
1119 759;
1215 759;
1349 759;
814 758;
1024 758;
1028 758;
1116 758;
1122 758;
1181 758;
1198 758;
1339 758;
626 757;
634 757;
641 757;
808 757;
825 757;
994 757;
1107 757;
1175 757;
1190 757;
1349 757;
615 756;
817 756;
1068 756;
1102 756;
1187 756;
1204 756;
662 755;
993 755;
1343 755;
601 754;
636 754;
656 754;
663 754;
1011 754;
1137 754;
1186 754;
1204 754;
566 753;
630 753;
679 753;
826 753;
829 753;
832 753;
1024 753;
1086 753;
1191 753;
609 752;
869 752;
871 752;
874 752;
887 752;
1032 752;
1155 752;
1236 752;
1372 752;
724 751;
733 751;
869 751;
1062 751;
1198 751;
1325 751;
1348 751;
1354 751;
533 750;
573 750;
625 750;
832 750;
842 750;
1077 750;
1156 750;
1224 750;
1329 750;
882 749;
954 749;
1048 749;
1062 749;
1233 749;
1338 749;
565 748;
969 748;
1012 748;
1108 748;
1331 748;
1332 748;
533 747;
535 747;
633 747;
721 747;
743 747;
1180 747;
1192 747;
563 746;
611 746;
660 746;
704 746;
979 746;
999 746;
1112 746;
1115 746;
1121 746;
1327 746;
639 745;
669 745;
729 745;
747 745;
819 745;
978 745;
1025 745;
1026 745;
1067 745;
1068 745;
1370 745;
535 744;
592 744;
593 744;
712 744;
721 744;
727 744;
825 744;
852 744;
1133 744;
1224 744;
1261 744;
836 743;
845 743;
1035 743;
1067 743;
1113 743;
1114 743;
1116 743;
1131 743;
1181 743;
566 742;
617 742;
668 742;
712 742;
995 742;
1028 742;
1256 742;
1342 742;
536 741;
698 741;
838 741;
1071 741;
1102 741;
1167 741;
1268 741;
1370 741;
707 740;
978 740;
990 740;
1029 740;
1069 740;
1081 740;
1121 740;
1132 740;
672 739;
678 739;
827 739;
841 739;
1079 739;
1102 739;
1206 739;
579 738;
605 738;
650 738;
681 738;
986 738;
1033 738;
1168 738;
535 737;
674 737;
679 737;
824 737;
885 737;
1015 737;
1112 737;
1139 737;
695 736;
714 736;
866 736;
1032 736;
1055 736;
1192 736;
1362 736;
528 735;
610 735;
640 735;
728 735;
826 735;
980 735;
1032 735;
1066 735;
1120 735;
1166 735;
1255 735;
528 734;
530 734;
808 734;
870 734;
978 734;
1078 734;
1082 734;
1156 734;
1159 734;
857 733;
975 733;
1079 733;
1340 733;
632 732;
662 732;
804 732;
814 732;
1028 732;
1256 732;
1353 732;
1373 732;
557 731;
562 731;
643 731;
690 731;
992 731;
1012 731;
1015 731;
1103 731;
1126 731;
1129 731;
1221 731;
1392 731;
680 730;
1176 730;
1177 730;
1192 730;
1196 730;
1232 730;
1237 730;
1271 730;
610 729;
984 729;
994 729;
1111 729;
1368 729;
555 728;
587 728;
806 728;
811 728;
841 728;
1003 728;
1128 728;
1369 728;
1378 728;
611 727;
667 727;
809 727;
1002 727;
1205 727;
1391 727;
650 726;
1011 726;
1150 726;
1368 726;
634 725;
733 725;
806 725;
1049 725;
525 724;
742 724;
761 724;
812 724;
1161 724;
1188 724;
1205 724;
538 723;
627 723;
733 723;
751 723;
994 723;
1119 723;
1383 723;
651 722;
683 722;
717 722;
1175 722;
1189 722;
1207 722;
1236 722;
1247 722;
1249 722;
642 721;
644 721;
815 721;
833 721;
1104 721;
585 720;
611 720;
837 720;
1044 720;
1074 720;
1384 720;
530 719;
533 719;
562 719;
568 719;
730 719;
1025 719;
1031 719;
1126 719;
1139 719;
1180 719;
1253 719;
1365 719;
1373 719;
737 718;
758 718;
813 718;
1042 718;
1133 718;
525 717;
736 717;
828 717;
830 717;
863 717;
1114 717;
521 716;
691 716;
751 716;
851 716;
1017 716;
1044 716;
1055 716;
1139 716;
1243 716;
1251 716;
1272 716;
1305 716;
542 715;
1172 715;
1207 715;
1299 715;
1370 715;
1396 715;
690 714;
750 714;
824 714;
1008 714;
1093 714;
1115 714;
1143 714;
547 713;
605 713;
1032 713;
1054 713;
1110 713;
1143 713;
1161 713;
1211 713;
1243 713;
724 712;
799 712;
848 712;
864 712;
1089 712;
1118 712;
1128 712;
1250 712;
1301 712;
596 711;
760 711;
1050 711;
1103 711;
1266 711;
1275 711;
543 710;
568 710;
727 710;
1090 710;
1135 710;
1243 710;
1257 710;
1358 710;
734 709;
1011 709;
1122 709;
1168 709;
712 708;
1043 708;
1129 708;
874 707;
1177 707;
1211 707;
1213 707;
1275 707;
1285 707;
1310 707;
1353 707;
1369 707;
523 706;
828 706;
838 706;
853 706;
1163 706;
1210 706;
514 705;
795 705;
861 705;
1014 705;
1055 705;
1241 705;
1362 705;
1369 705;
1373 705;
1393 705;
703 704;
855 704;
875 704;
1033 704;
1116 704;
1184 704;
1260 704;
1283 704;
533 703;
537 703;
1128 703;
1173 703;
1196 703;
1208 703;
1183 702;
525 701;
833 701;
856 701;
857 701;
1004 701;
1029 701;
1086 701;
1092 701;
1377 701;
523 700;
549 700;
739 700;
741 700;
743 700;
844 700;
853 700;
994 700;
1215 700;
1250 700;
1313 700;
1359 700;
1391 700;
527 699;
579 699;
1048 699;
1186 699;
1228 699;
739 698;
759 698;
787 698;
791 698;
861 698;
1063 698;
1082 698;
1124 698;
1244 698;
1256 698;
1274 698;
1278 698;
1374 698;
581 697;
729 697;
738 697;
779 697;
1129 697;
1187 697;
1231 697;
1240 697;
1309 697;
1310 697;
1364 697;
529 696;
534 696;
716 696;
824 696;
849 696;
1119 696;
1192 696;
1411 696;
535 695;
571 695;
850 695;
860 695;
1038 695;
1060 695;
1088 695;
1105 695;
562 694;
815 694;
1059 694;
1114 694;
1115 694;
1228 694;
1293 694;
562 693;
1048 693;
1079 693;
1265 693;
1272 693;
1373 693;
827 692;
1012 692;
1191 692;
1292 692;
1280 691;
501 690;
534 690;
557 690;
578 690;
580 690;
785 690;
861 690;
1200 690;
1389 690;
746 689;
763 689;
1048 689;
1248 689;
1390 689;
500 688;
570 688;
762 688;
783 688;
833 688;
1112 688;
1268 688;
504 687;
543 687;
789 687;
832 687;
1037 687;
1108 687;
1374 687;
1402 687;
836 686;
842 686;
1064 686;
1208 686;
1231 686;
1331 686;
1406 686;
526 685;
565 685;
586 685;
794 685;
828 685;
1009 685;
1020 685;
1212 685;
1219 685;
1300 685;
551 684;
561 684;
801 684;
815 684;
1035 684;
1116 684;
1296 684;
1413 684;
571 683;
773 683;
780 683;
784 683;
831 683;
1000 683;
1013 683;
1034 683;
1328 683;
507 682;
587 682;
744 682;
778 682;
787 682;
821 682;
830 682;
1003 682;
1037 682;
1133 682;
1201 682;
1213 682;
1257 682;
1274 682;
1335 682;
527 681;
743 681;
1024 681;
1065 681;
1221 681;
1302 681;
1387 681;
790 680;
1335 680;
560 679;
768 679;
850 679;
1033 679;
1050 679;
1340 679;
503 678;
793 678;
795 678;
868 678;
1044 678;
1285 678;
558 677;
733 677;
834 677;
1005 677;
1048 677;
1247 677;
1333 677;
508 676;
814 676;
831 676;
1000 676;
1057 676;
1065 676;
1066 676;
1237 676;
1240 676;
1383 676;
497 675;
737 675;
753 675;
1239 675;
1319 675;
508 674;
524 674;
549 674;
1006 674;
1218 674;
1243 674;
514 673;
566 673;
803 673;
816 673;
1258 673;
1304 673;
1395 673;
499 672;
775 672;
1031 672;
1040 672;
1255 672;
1408 672;
774 671;
835 671;
1018 671;
1019 671;
1312 671;
1395 671;
1407 671;
991 670;
1272 670;
1294 670;
503 669;
1213 669;
1250 669;
1287 669;
552 668;
801 668;
743 667;
846 667;
851 667;
1003 667;
1006 667;
1019 667;
1396 667;
755 666;
770 666;
845 666;
1293 666;
1323 666;
1347 666;
1380 666;
506 665;
554 665;
757 665;
785 665;
800 665;
849 665;
872 665;
1037 665;
1238 665;
805 664;
807 664;
1047 664;
1391 664;
505 663;
768 663;
812 663;
874 663;
1030 663;
758 662;
1298 662;
1316 662;
542 661;
1221 661;
1225 661;
1233 661;
1227 660;
1236 660;
1237 660;
1256 660;
1260 660;
1287 660;
1402 660;
506 659;
996 659;
1023 659;
1310 659;
1333 659;
1391 659;
761 658;
765 658;
990 658;
1023 658;
1030 658;
1285 658;
787 657;
834 657;
1000 657;
1011 657;
1273 657;
744 656;
767 656;
1034 656;
842 655;
1002 655;
1049 655;
1317 655;
1339 655;
1408 655;
1022 654;
541 653;
760 653;
761 653;
771 653;
807 653;
810 653;
823 653;
1277 653;
836 652;
1042 652;
1228 652;
1306 652;
747 651;
1014 651;
1298 651;
1393 651;
510 650;
838 650;
848 650;
1016 650;
1041 650;
1065 650;
501 649;
774 649;
782 649;
826 649;
848 649;
1011 649;
1042 649;
1229 649;
1299 649;
1428 649;
790 648;
1020 648;
1165 648;
1288 648;
1303 648;
512 647;
848 647;
1002 647;
514 646;
774 646;
788 646;
528 644;
1159 644;
1347 644;
501 643;
883 643;
1158 643;
1428 643;
505 642;
764 642;
798 642;
843 642;
1186 642;
1295 642;
784 641;
1168 641;
1424 641;
510 640;
846 640;
880 640;
885 640;
1013 640;
1028 640;
1289 640;
1307 640;
1329 640;
1331 640;
529 639;
826 639;
850 639;
892 639;
1327 639;
806 638;
997 638;
998 638;
1007 638;
774 637;
824 637;
857 637;
872 637;
894 637;
1031 637;
1046 637;
1182 637;
1360 637;
1416 637;
756 636;
827 636;
887 636;
1173 636;
497 635;
525 635;
547 635;
1431 635;
546 634;
841 634;
900 634;
1026 634;
1356 634;
774 633;
870 633;
1018 633;
546 632;
549 632;
809 632;
854 632;
873 632;
1170 632;
1318 632;
1347 632;
556 631;
769 631;
770 631;
823 631;
859 631;
933 631;
984 631;
1015 631;
1158 631;
1297 631;
1308 631;
549 630;
819 630;
1165 630;
1355 630;
499 629;
888 629;
1016 629;
1031 629;
1169 629;
1184 629;
1326 629;
1346 629;
1410 629;
799 628;
829 628;
1159 628;
1433 628;
766 627;
963 627;
983 627;
1189 627;
872 626;
900 626;
979 626;
1173 626;
1353 626;
497 625;
788 625;
797 625;
826 625;
894 625;
967 625;
1200 625;
1209 625;
1320 625;
556 624;
838 624;
1161 624;
1174 624;
1190 624;
1207 624;
1299 624;
927 623;
1409 623;
555 622;
777 622;
817 622;
875 622;
885 622;
1199 622;
1326 622;
1370 622;
1399 622;
710 621;
867 621;
868 621;
874 621;
972 621;
989 621;
1212 621;
1351 621;
549 620;
564 620;
831 620;
995 620;
1146 620;
1374 620;
555 619;
877 619;
909 619;
912 619;
930 619;
964 619;
1366 619;
497 618;
869 618;
874 618;
895 618;
946 618;
1024 618;
1151 618;
1309 618;
1347 618;
815 617;
904 617;
914 617;
968 617;
1179 617;
1315 617;
1368 617;
1420 617;
767 616;
1036 616;
776 615;
855 615;
928 615;
961 615;
1175 615;
1191 615;
1410 615;
1412 615;
810 614;
857 614;
1001 614;
1006 614;
1316 614;
491 613;
541 613;
965 613;
1323 613;
1341 613;
666 612;
779 612;
816 612;
677 611;
812 611;
843 611;
903 611;
1201 611;
501 610;
666 610;
678 610;
722 610;
998 610;
1172 610;
663 609;
773 609;
801 609;
920 609;
933 609;
992 609;
1004 609;
1420 609;
662 608;
714 608;
879 608;
884 608;
944 608;
1409 608;
557 607;
657 607;
699 607;
774 607;
819 607;
1201 607;
1329 607;
555 606;
819 606;
847 606;
852 606;
932 606;
1344 606;
541 605;
665 605;
789 605;
877 605;
907 605;
957 605;
1189 605;
1212 605;
505 604;
920 604;
1004 604;
1011 604;
1026 604;
1074 604;
1430 604;
491 603;
493 603;
693 603;
843 603;
859 603;
864 603;
949 603;
955 603;
1132 603;
1150 603;
1210 603;
1320 603;
968 602;
984 602;
1138 602;
673 601;
725 601;
791 601;
1311 601;
1416 601;
487 600;
1223 600;
1324 600;
1410 600;
826 599;
833 599;
857 599;
945 599;
963 599;
1037 599;
1055 599;
1195 599;
524 598;
695 598;
993 598;
1010 598;
1021 598;
1048 598;
1067 598;
1077 598;
1148 598;
1229 598;
1416 598;
524 597;
714 597;
734 597;
861 597;
914 597;
983 597;
1359 597;
874 596;
1351 596;
1357 596;
1436 596;
498 595;
707 595;
923 595;
972 595;
1010 595;
1055 595;
1142 595;
1310 595;
800 594;
810 594;
1045 594;
1118 594;
1192 594;
1229 594;
559 593;
928 593;
931 593;
991 593;
1158 593;
550 592;
985 592;
994 592;
1082 592;
1217 592;
1332 592;
729 591;
818 591;
830 591;
875 591;
915 591;
942 591;
1059 591;
1307 591;
920 590;
1057 590;
1122 590;
1133 590;
1136 590;
1175 590;
1188 590;
1419 590;
717 589;
788 589;
832 589;
886 589;
896 589;
911 589;
944 589;
971 589;
1315 589;
1366 589;
1416 589;
834 588;
913 588;
1138 588;
1425 588;
557 587;
564 587;
803 587;
871 587;
900 587;
948 587;
1016 587;
1057 587;
1146 587;
1344 587;
1378 587;
558 586;
670 586;
853 586;
985 586;
1002 586;
680 585;
795 585;
823 585;
958 585;
1034 585;
1168 585;
1335 585;
486 584;
656 584;
947 584;
1030 584;
1038 584;
1123 584;
1167 584;
1225 584;
1360 584;
911 583;
919 583;
1024 583;
1147 583;
1202 583;
1240 583;
673 582;
688 582;
856 582;
561 581;
956 581;
1020 581;
1186 581;
1237 581;
1351 581;
1364 581;
1430 581;
806 580;
810 580;
929 580;
1009 580;
1055 580;
1316 580;
1366 580;
1379 580;
725 579;
733 579;
792 579;
917 579;
962 579;
1062 579;
1078 579;
1146 579;
1149 579;
1174 579;
634 578;
700 578;
856 578;
946 578;
1182 578;
1336 578;
649 577;
674 577;
996 577;
1083 577;
1125 577;
1213 577;
1361 577;
1371 577;
787 576;
814 576;
835 576;
850 576;
927 576;
959 576;
1035 576;
1118 576;
1125 576;
1226 576;
1313 576;
819 575;
885 575;
917 575;
918 575;
1188 575;
1207 575;
1316 575;
1355 575;
1393 575;
638 574;
682 574;
720 574;
725 574;
821 574;
966 574;
987 574;
1016 574;
1060 574;
1098 574;
1111 574;
1136 574;
1204 574;
1342 574;
919 573;
925 573;
991 573;
1036 573;
1126 573;
1163 573;
1192 573;
1205 573;
1238 573;
1426 573;
693 572;
694 572;
733 572;
993 572;
1084 572;
1087 572;
1123 572;
1187 572;
1318 572;
1345 572;
636 571;
704 571;
722 571;
985 571;
1014 571;
1029 571;
1032 571;
1151 571;
1376 571;
810 570;
837 570;
894 570;
963 570;
1000 570;
1127 570;
1250 570;
648 569;
796 569;
817 569;
837 569;
888 569;
895 569;
914 569;
971 569;
1013 569;
1106 569;
1134 569;
1233 569;
1247 569;
734 568;
833 568;
843 568;
863 568;
908 568;
910 568;
930 568;
937 568;
960 568;
1048 568;
1061 568;
1342 568;
708 567;
720 567;
809 567;
882 567;
1050 567;
1085 567;
1164 567;
1428 567;
705 566;
816 566;
869 566;
942 566;
964 566;
1059 566;
1063 566;
1146 566;
1313 566;
1343 566;
673 565;
866 565;
881 565;
892 565;
907 565;
939 565;
943 565;
947 565;
985 565;
1068 565;
1109 565;
1157 565;
1168 565;
1174 565;
1313 565;
1346 565;
1387 565;
666 564;
792 564;
853 564;
869 564;
891 564;
1011 564;
1013 564;
1023 564;
1101 564;
1138 564;
1173 564;
1211 564;
1348 564;
1389 564;
856 563;
1031 563;
1071 563;
1077 563;
1145 563;
1154 563;
1199 563;
1383 563;
712 562;
1028 562;
1067 562;
1119 562;
1120 562;
1195 562;
1258 562;
1345 562;
975 561;
1026 561;
1128 561;
1132 561;
1325 561;
691 560;
829 560;
904 560;
995 560;
1040 560;
1063 560;
1147 560;
1197 560;
516 559;
983 559;
1114 559;
1136 559;
1223 559;
1328 559;
892 558;
914 558;
1162 558;
1175 558;
670 557;
674 557;
714 557;
732 557;
846 557;
867 557;
1146 557;
1383 557;
670 556;
697 556;
787 556;
860 556;
878 556;
1049 556;
1082 556;
1121 556;
1227 556;
1253 556;
1426 556;
657 555;
713 555;
874 555;
876 555;
884 555;
1024 555;
1048 555;
1086 555;
1125 555;
1258 555;
683 554;
853 554;
875 554;
911 554;
1004 554;
1067 554;
1138 554;
1174 554;
1237 554;
1369 554;
797 553;
886 553;
994 553;
1054 553;
1055 553;
1094 553;
1255 553;
1379 553;
678 552;
830 552;
832 552;
918 552;
1056 552;
1060 552;
1089 552;
1129 552;
1134 552;
1193 552;
1341 552;
1366 552;
1432 552;
504 551;
667 551;
738 551;
1018 551;
1107 551;
1137 551;
1144 551;
1377 551;
691 550;
728 550;
870 550;
1019 550;
1131 550;
1148 550;
1161 550;
1342 550;
1385 550;
1431 550;
503 549;
927 549;
963 549;
996 549;
1198 549;
681 548;
689 548;
873 548;
899 548;
1011 548;
1117 548;
1118 548;
1156 548;
1167 548;
1177 548;
1351 548;
1392 548;
825 547;
830 547;
891 547;
965 547;
1023 547;
1104 547;
1117 547;
1136 547;
1266 547;
1314 547;
1357 547;
500 546;
938 546;
1022 546;
1069 546;
1086 546;
1208 546;
1369 546;
1388 546;
496 545;
719 545;
854 545;
893 545;
966 545;
1434 545;
680 544;
837 544;
893 544;
1084 544;
1351 544;
1390 544;
703 543;
918 543;
948 543;
961 543;
974 543;
1190 543;
1194 543;
1213 543;
1267 543;
923 542;
967 542;
1021 542;
1081 542;
1146 542;
1209 542;
1210 542;
1220 542;
1354 542;
1400 542;
687 541;
703 541;
989 541;
1127 541;
1151 541;
1350 541;
1369 541;
736 540;
838 540;
870 540;
1048 540;
1119 540;
1121 540;
1144 540;
1148 540;
1162 540;
1249 540;
1346 540;
1358 540;
1377 540;
1384 540;
884 539;
915 539;
943 539;
1027 539;
1198 539;
1258 539;
1435 539;
737 538;
865 538;
900 538;
1021 538;
1111 538;
1393 538;
690 537;
707 537;
729 537;
839 537;
860 537;
862 537;
939 537;
1055 537;
1095 537;
1100 537;
1230 537;
904 536;
944 536;
1010 536;
1071 536;
1390 536;
863 535;
930 535;
1040 535;
1222 535;
1360 535;
858 534;
983 534;
992 534;
1129 534;
1133 534;
1149 534;
868 533;
1074 533;
1097 533;
1163 533;
1399 533;
740 532;
819 532;
876 532;
934 532;
1124 532;
1189 532;
1358 532;
706 531;
713 531;
862 531;
918 531;
926 531;
1038 531;
1178 531;
1194 531;
1219 531;
843 530;
940 530;
967 530;
1047 530;
1083 530;
1100 530;
1355 530;
840 529;
1194 529;
1252 529;
1366 529;
1404 529;
711 528;
712 528;
995 528;
1130 528;
716 527;
908 527;
1046 527;
1354 527;
705 526;
832 526;
872 526;
1351 526;
1360 526;
1048 525;
1082 525;
1272 525;
1359 525;
1368 525;
737 524;
824 524;
1041 524;
1049 524;
1149 524;
1194 524;
1222 524;
1377 524;
841 523;
863 523;
922 523;
949 523;
978 523;
1203 523;
1352 523;
1372 523;
1393 523;
1408 523;
740 522;
1031 522;
1062 522;
1107 522;
1236 522;
1271 522;
490 521;
826 521;
828 521;
893 521;
977 521;
1084 521;
1149 521;
1399 521;
937 520;
984 520;
1004 520;
1192 520;
1391 520;
880 519;
961 519;
1107 519;
1145 519;
947 518;
1021 518;
1173 518;
1227 518;
1253 518;
1265 518;
895 517;
1147 517;
1196 517;
1205 517;
1229 517;
1248 517;
1269 517;
1358 517;
1395 517;
1400 517;
506 516;
725 516;
853 516;
878 516;
885 516;
944 516;
1108 516;
1197 516;
1207 516;
1240 516;
1355 516;
603 515;
734 515;
1013 515;
1056 515;
1158 515;
1204 515;
1256 515;
750 514;
850 514;
885 514;
919 514;
970 514;
1014 514;
1065 514;
1098 514;
491 513;
590 513;
717 513;
848 513;
885 513;
959 513;
1027 513;
1094 513;
710 512;
838 512;
889 512;
955 512;
1060 512;
1065 512;
1075 512;
1083 512;
1103 512;
1112 512;
1166 512;
1240 512;
1252 512;
1366 512;
1396 512;
741 511;
848 511;
850 511;
1057 511;
1065 511;
1145 511;
1249 511;
591 510;
936 510;
966 510;
1000 510;
1053 510;
1168 510;
1172 510;
1403 510;
502 509;
870 509;
1038 509;
1068 509;
1179 509;
1376 509;
1059 508;
1079 508;
1120 508;
1217 508;
1230 508;
483 507;
486 507;
862 507;
891 507;
1200 507;
1409 507;
596 506;
731 506;
745 506;
1009 506;
1033 506;
1035 506;
1041 506;
1365 506;
1379 506;
1410 506;
497 505;
551 505;
723 505;
862 505;
875 505;
971 505;
1124 505;
887 504;
934 504;
969 504;
1011 504;
1062 504;
1064 504;
1177 504;
1209 504;
1271 504;
1273 504;
1379 504;
908 503;
998 503;
1085 503;
1249 503;
1384 503;
485 502;
871 502;
880 502;
903 502;
1067 502;
1164 502;
1175 502;
1230 502;
1402 502;
636 501;
742 501;
906 501;
990 501;
1093 501;
1203 501;
1258 501;
1373 501;
507 500;
558 500;
728 500;
752 500;
856 500;
972 500;
1001 500;
1120 500;
1125 500;
1203 500;
1204 500;
1396 500;
558 499;
623 499;
736 499;
909 499;
949 499;
1080 499;
1109 499;
1174 499;
1204 499;
1210 499;
1267 499;
490 498;
573 498;
740 498;
901 498;
914 498;
1018 498;
1066 498;
1192 498;
575 497;
720 497;
1104 497;
1404 497;
602 496;
752 496;
888 496;
951 496;
1030 496;
1072 496;
1079 496;
1092 496;
1124 496;
1214 496;
1215 496;
1226 496;
1235 496;
1377 496;
1378 496;
503 495;
640 495;
714 495;
1092 495;
1102 495;
506 494;
545 494;
641 494;
642 494;
731 494;
941 494;
960 494;
1011 494;
1156 494;
1404 494;
556 493;
650 493;
741 493;
895 493;
1097 493;
1157 493;
1264 493;
1364 493;
1365 493;
1371 493;
566 492;
603 492;
639 492;
874 492;
892 492;
895 492;
989 492;
997 492;
1047 492;
1048 492;
1053 492;
1058 492;
1118 492;
1219 492;
1229 492;
582 491;
882 491;
890 491;
966 491;
1065 491;
1170 491;
1198 491;
1267 491;
909 490;
1071 490;
1274 490;
1294 490;
539 489;
741 489;
985 489;
996 489;
1007 489;
1051 489;
1226 489;
1383 489;
1390 489;
537 488;
717 488;
898 488;
1044 488;
1048 488;
1198 488;
1248 488;
1373 488;
1380 488;
1410 488;
924 487;
1026 487;
1104 487;
1114 487;
1137 487;
1255 487;
1398 487;
609 486;
1019 486;
1111 486;
1412 486;
612 485;
718 485;
892 485;
925 485;
1056 485;
1063 485;
1106 485;
1191 485;
1250 485;
916 484;
1074 484;
1377 484;
558 483;
746 483;
1032 483;
1103 483;
1172 483;
1196 483;
897 482;
1066 482;
1120 482;
584 481;
623 481;
646 481;
912 481;
923 481;
937 481;
1011 481;
1077 481;
1091 481;
1214 481;
1266 481;
554 480;
563 480;
624 480;
739 480;
943 480;
981 480;
1114 480;
1411 480;
497 479;
594 479;
627 479;
740 479;
1015 479;
1184 479;
495 478;
560 478;
625 478;
974 478;
1021 478;
1135 478;
1243 478;
573 477;
734 477;
870 477;
940 477;
961 477;
1002 477;
1097 477;
622 476;
650 476;
893 476;
901 476;
1105 476;
1152 476;
1172 476;
1217 476;
1270 476;
1284 476;
1385 476;
559 475;
1020 475;
1044 475;
1169 475;
540 474;
578 474;
871 474;
896 474;
952 474;
1058 474;
1368 474;
1395 474;
501 473;
552 473;
630 473;
777 473;
908 473;
1146 473;
1267 473;
1298 473;
1407 473;
973 472;
975 472;
1183 472;
1207 472;
1231 472;
1244 472;
1278 472;
510 471;
728 471;
1207 471;
1257 471;
561 470;
717 470;
885 470;
936 470;
1029 470;
1185 470;
1389 470;
732 469;
777 469;
970 469;
1168 469;
1224 469;
627 468;
633 468;
719 468;
783 468;
917 468;
1019 468;
1026 468;
1043 468;
1082 468;
1169 468;
1236 468;
947 467;
1133 467;
1374 467;
774 466;
1057 466;
1112 466;
659 465;
673 465;
961 465;
1046 465;
1056 465;
1178 465;
1194 465;
565 464;
713 464;
726 464;
767 464;
961 464;
1081 464;
1128 464;
1407 464;
508 463;
758 463;
776 463;
871 463;
944 463;
954 463;
964 463;
998 463;
1035 463;
1252 463;
633 462;
643 462;
707 462;
1011 462;
1071 462;
1098 462;
1126 462;
1215 462;
1230 462;
1235 462;
546 461;
892 461;
1105 461;
1144 461;
1246 461;
1399 461;
642 460;
764 460;
884 460;
897 460;
1115 460;
1166 460;
510 459;
549 459;
878 459;
913 459;
928 459;
1035 459;
1043 459;
1152 459;
1223 459;
1373 459;
1375 459;
491 458;
640 458;
722 458;
981 458;
999 458;
1077 458;
1139 458;
1173 458;
1179 458;
1191 458;
1275 458;
1392 458;
1413 458;
539 457;
555 457;
715 457;
736 457;
972 457;
1021 457;
1053 457;
1156 457;
1276 457;
1390 457;
512 456;
871 456;
942 456;
995 456;
1039 456;
1167 456;
1211 456;
495 455;
546 455;
716 455;
734 455;
795 455;
1009 455;
1146 455;
1215 455;
1228 455;
1253 455;
1378 455;
507 454;
518 454;
790 454;
891 454;
966 454;
1053 454;
1114 454;
1117 454;
1149 454;
1407 454;
515 453;
538 453;
718 453;
950 453;
1064 453;
1173 453;
1255 453;
491 452;
889 452;
940 452;
1190 452;
1228 452;
1274 452;
926 451;
1003 451;
1026 451;
1073 451;
1077 451;
1081 451;
1084 451;
1130 451;
1164 451;
1225 451;
1306 451;
1366 451;
544 450;
656 450;
778 450;
971 450;
1037 450;
1066 450;
1145 450;
1204 450;
1228 450;
963 449;
985 449;
990 449;
1225 449;
1364 449;
491 448;
721 448;
801 448;
908 448;
987 448;
994 448;
1253 448;
1297 448;
501 447;
872 447;
989 447;
1101 447;
1169 447;
1257 447;
1366 447;
536 446;
794 446;
951 446;
1050 446;
510 445;
661 445;
926 445;
1053 445;
1171 445;
1396 445;
516 444;
535 444;
810 444;
991 444;
1263 444;
1369 444;
916 443;
920 443;
936 443;
1035 443;
1385 443;
532 442;
668 442;
784 442;
902 442;
1117 442;
655 441;
809 441;
906 441;
1000 441;
1156 441;
1225 441;
1400 441;
511 440;
513 440;
909 440;
1232 440;
1412 440;
507 439;
787 439;
967 439;
1116 439;
1129 439;
1156 439;
1241 439;
1258 439;
1308 439;
655 438;
715 438;
717 438;
790 438;
806 438;
982 438;
1167 438;
1279 438;
648 437;
779 437;
1153 437;
1253 437;
1263 437;
667 436;
822 436;
1018 436;
1175 436;
512 435;
795 435;
874 435;
891 435;
944 435;
1120 435;
1131 435;
1142 435;
1241 435;
1374 435;
1409 435;
532 434;
937 434;
1000 434;
1019 434;
1286 434;
497 433;
507 433;
649 433;
665 433;
899 433;
901 433;
971 433;
979 433;
1044 433;
1114 433;
1142 433;
1174 433;
1233 433;
1243 433;
1402 433;
716 432;
877 432;
1002 432;
1098 432;
519 431;
523 431;
669 431;
922 431;
950 431;
1032 431;
1201 431;
1246 431;
529 430;
531 430;
834 430;
1116 430;
1136 430;
1254 430;
1265 430;
1375 430;
1390 430;
718 429;
1191 429;
1236 429;
882 428;
948 428;
1030 428;
1048 428;
1055 428;
1059 428;
1391 428;
653 427;
828 427;
996 427;
1006 427;
1012 427;
1049 427;
1084 427;
1122 427;
1168 427;
528 426;
661 426;
999 426;
1015 426;
1076 426;
1269 426;
1405 426;
655 425;
796 425;
915 425;
917 425;
993 425;
1203 425;
1403 425;
509 424;
667 424;
1129 424;
794 423;
889 423;
919 423;
990 423;
1127 423;
1176 423;
1195 423;
1203 423;
1305 423;
1363 423;
1380 423;
679 422;
881 422;
896 422;
1106 422;
1114 422;
1159 422;
1205 422;
1277 422;
1286 422;
537 421;
796 421;
838 421;
1018 421;
1150 421;
1169 421;
1250 421;
1259 421;
655 420;
658 420;
816 420;
877 420;
900 420;
916 420;
919 420;
920 420;
1001 420;
1013 420;
1107 420;
1129 420;
1248 420;
1260 420;
1367 420;
1386 420;
1398 420;
835 419;
921 419;
982 419;
1018 419;
1158 419;
1167 419;
1257 419;
1258 419;
654 418;
662 418;
684 418;
892 418;
1052 418;
1171 418;
1244 418;
1270 418;
1384 418;
1410 418;
992 417;
1093 417;
1114 417;
1355 417;
520 416;
682 416;
854 416;
1037 416;
1123 416;
1205 416;
888 415;
905 415;
924 415;
937 415;
1028 415;
1087 415;
1200 415;
1383 415;
681 414;
878 414;
965 414;
1029 414;
1134 414;
1142 414;
1185 414;
1309 414;
1366 414;
1381 414;
516 413;
537 413;
838 413;
929 413;
968 413;
998 413;
1047 413;
1064 413;
1072 413;
1077 413;
1081 413;
1102 413;
1145 413;
1295 413;
1353 413;
1354 413;
1406 413;
804 412;
815 412;
841 412;
859 412;
917 412;
1041 412;
1051 412;
1128 412;
501 411;
540 411;
667 411;
1005 411;
1046 411;
1132 411;
1163 411;
1396 411;
658 410;
853 410;
864 410;
901 410;
1023 410;
1266 410;
501 409;
535 409;
648 409;
908 409;
1306 409;
1388 409;
1409 409;
512 408;
540 408;
976 408;
992 408;
1006 408;
1024 408;
1058 408;
1148 408;
1206 408;
1367 408;
869 407;
916 407;
1042 407;
1105 407;
1152 407;
1162 407;
1253 407;
1305 407;
1405 407;
674 406;
873 406;
921 406;
963 406;
1004 406;
1047 406;
1138 406;
1158 406;
1164 406;
1245 406;
1249 406;
1254 406;
1294 406;
1379 406;
988 405;
993 405;
1056 405;
1067 405;
1081 405;
1083 405;
1162 405;
1182 405;
679 404;
889 404;
980 404;
1052 404;
1155 404;
511 403;
985 403;
1270 403;
1362 403;
1366 403;
647 402;
660 402;
850 402;
852 402;
855 402;
913 402;
962 402;
1043 402;
1288 402;
1358 402;
847 401;
908 401;
1073 401;
1090 401;
527 400;
653 400;
904 400;
1029 400;
1163 400;
1182 400;
1189 400;
1270 400;
516 399;
889 399;
1000 399;
1008 399;
1038 399;
1060 399;
1277 399;
508 398;
824 398;
837 398;
970 398;
1127 398;
1145 398;
1263 398;
840 397;
870 397;
887 397;
922 397;
1026 397;
1060 397;
1168 397;
1268 397;
1348 397;
856 396;
967 396;
1016 396;
1154 396;
1287 396;
1368 396;
510 395;
537 395;
855 395;
913 395;
979 395;
1291 395;
1294 395;
828 394;
886 394;
998 394;
1042 394;
1047 394;
1105 394;
1165 394;
533 393;
681 393;
829 393;
1074 393;
1091 393;
1170 393;
1235 393;
842 392;
851 392;
882 392;
1001 392;
1070 392;
1178 392;
1193 392;
1296 392;
1361 392;
645 391;
657 391;
837 391;
969 391;
1013 391;
1050 391;
1128 391;
1130 391;
1386 391;
1397 391;
641 390;
1047 390;
1245 390;
1253 390;
1264 390;
1275 390;
527 389;
852 389;
865 389;
981 389;
1011 389;
1131 389;
1242 389;
1251 389;
1299 389;
1310 389;
1378 389;
908 388;
1106 388;
867 387;
1021 387;
1047 387;
1080 387;
1090 387;
1254 387;
1264 387;
1385 387;
1022 386;
1039 386;
1084 386;
1139 386;
849 385;
859 385;
1088 385;
1102 385;
539 384;
1159 384;
1267 384;
536 383;
890 383;
1092 383;
1158 383;
1257 383;
1386 383;
523 382;
537 382;
865 382;
1266 382;
1314 382;
515 381;
550 381;
905 381;
1079 381;
1157 381;
1255 381;
1313 381;
1386 381;
1392 381;
668 380;
890 380;
901 380;
973 380;
988 380;
1050 380;
1138 380;
1398 380;
851 379;
859 379;
878 379;
988 379;
996 379;
1127 379;
1174 379;
1178 379;
1042 378;
1089 378;
1117 378;
1165 378;
1131 377;
1151 377;
1283 377;
1306 377;
528 376;
556 376;
652 376;
966 376;
1079 376;
1091 376;
1170 376;
1266 376;
1389 376;
632 375;
862 375;
1270 375;
1376 375;
545 374;
628 374;
906 374;
907 374;
1185 374;
1371 374;
661 373;
995 373;
545 372;
868 372;
996 372;
1023 372;
1031 372;
1046 372;
1265 372;
1293 372;
557 371;
657 371;
980 371;
1060 371;
1258 371;
890 370;
895 370;
903 370;
1070 370;
1154 370;
1165 370;
1249 370;
1306 370;
1390 370;
576 369;
665 369;
899 369;
1015 369;
1148 369;
1197 369;
603 368;
622 368;
548 367;
580 367;
595 367;
901 367;
1025 367;
1055 367;
1073 367;
1075 367;
1179 367;
1257 367;
1383 367;
597 366;
1006 366;
1080 366;
1097 366;
1151 366;
1384 366;
551 365;
627 365;
1045 365;
1050 365;
1092 365;
1182 365;
541 364;
910 364;
997 364;
1121 364;
1258 364;
525 363;
594 363;
888 363;
900 363;
914 363;
979 363;
1306 363;
644 362;
902 362;
974 362;
1003 362;
1091 362;
1252 362;
1303 362;
1370 362;
1381 362;
624 361;
973 361;
1003 361;
1121 361;
1127 361;
1286 361;
1075 360;
1111 360;
1276 360;
1393 360;
517 359;
539 359;
577 359;
591 359;
619 359;
629 359;
1076 359;
1253 359;
1270 359;
1296 359;
1395 359;
527 358;
545 358;
635 358;
893 358;
1032 358;
1132 358;
1240 358;
1244 358;
1269 358;
1295 358;
1391 358;
557 357;
560 357;
1007 357;
1030 357;
1040 357;
1095 357;
1180 357;
1394 357;
771 356;
1012 356;
1129 356;
1135 356;
1174 356;
1263 356;
1268 356;
1290 356;
1391 356;
594 355;
990 355;
995 355;
1179 355;
1292 355;
1012 354;
1123 354;
1301 354;
533 353;
549 353;
1046 353;
1152 353;
1246 353;
761 352;
1047 352;
1057 352;
1164 352;
622 351;
778 351;
1013 351;
1015 351;
1021 351;
1066 351;
1094 351;
1139 351;
1171 351;
1368 351;
574 350;
1038 350;
1298 350;
582 349;
611 349;
1033 349;
1104 349;
1142 349;
1177 349;
1238 349;
1266 349;
1267 349;
777 348;
1008 348;
1108 348;
1390 348;
789 347;
976 347;
1138 347;
1174 347;
1294 347;
552 346;
892 346;
1030 346;
1097 346;
1238 346;
1241 346;
1245 346;
778 345;
990 345;
1026 345;
1065 345;
780 344;
893 344;
1021 344;
1032 344;
1172 344;
1227 344;
1290 344;
526 343;
530 343;
783 343;
985 343;
997 343;
1023 343;
1038 343;
1119 343;
1145 343;
979 342;
982 342;
1150 342;
1259 342;
537 341;
1057 341;
1135 341;
1189 341;
1246 341;
1133 340;
1143 340;
1166 340;
1258 340;
541 339;
884 339;
537 338;
554 338;
1037 338;
1046 338;
1113 338;
1252 338;
1257 338;
1309 338;
529 337;
534 337;
782 337;
1291 337;
1122 336;
1164 336;
1177 336;
1027 335;
1134 335;
1141 335;
538 334;
784 334;
988 334;
1019 334;
1159 334;
1173 333;
1185 333;
1291 333;
543 332;
1003 332;
1018 332;
1104 332;
1053 331;
1130 331;
773 330;
1035 330;
1272 330;
1367 330;
541 329;
736 329;
741 329;
791 329;
792 329;
1014 329;
1047 329;
1177 329;
734 328;
764 328;
1055 328;
1102 328;
545 327;
784 327;
990 327;
1018 327;
1096 327;
1098 327;
1142 327;
1177 327;
1291 327;
1002 326;
1025 326;
1028 326;
1369 326;
781 325;
1029 325;
1261 325;
1379 325;
781 324;
786 324;
1094 324;
1107 324;
1167 324;
1185 324;
774 323;
1028 323;
1145 323;
1191 323;
1261 323;
1377 323;
758 322;
788 322;
794 322;
1021 322;
1023 322;
1045 322;
1276 322;
725 321;
1051 320;
1130 320;
1267 320;
1272 320;
1297 320;
800 319;
1105 319;
1137 319;
1152 319;
1184 319;
1025 318;
1286 318;
567 317;
779 317;
1002 317;
1094 317;
1102 317;
1286 317;
1052 316;
1056 316;
1095 316;
1171 316;
1107 315;
991 314;
1037 314;
1139 314;
1162 314;
1192 314;
1304 314;
556 313;
749 312;
750 312;
766 312;
769 312;
1002 312;
1146 312;
1149 312;
1157 312;
1144 311;
1147 311;
1156 311;
1286 311;
1295 311;
759 310;
796 310;
801 310;
999 310;
1298 310;
1358 310;
1192 309;
783 308;
1022 308;
1108 308;
1119 308;
812 307;
1111 307;
1140 307;
552 306;
555 306;
1165 306;
773 305;
1173 305;
1273 305;
565 304;
577 304;
1192 304;
1269 304;
785 303;
1036 303;
1132 303;
1136 303;
574 302;
785 302;
1135 302;
1165 302;
1293 302;
762 301;
990 301;
1149 301;
1160 301;
1257 301;
572 300;
794 300;
799 300;
1366 300;
546 299;
1141 299;
1158 299;
1278 299;
826 298;
1190 298;
573 297;
1113 297;
1144 297;
578 296;
769 296;
828 296;
1022 296;
1034 296;
1043 296;
1138 295;
1298 295;
780 294;
1029 294;
1036 294;
1143 294;
1261 294;
1007 293;
1013 293;
1193 293;
567 292;
985 292;
1153 292;
1294 292;
1281 291;
566 290;
779 290;
1045 290;
1107 290;
1247 290;
553 289;
579 289;
734 289;
1025 289;
1053 289;
1056 289;
1094 289;
1141 289;
776 288;
590 287;
1258 287;
773 286;
776 286;
841 286;
1037 286;
1153 286;
772 285;
848 285;
1029 285;
1041 285;
1153 285;
1158 285;
1253 285;
862 284;
959 284;
1104 284;
1259 284;
735 283;
939 283;
966 283;
985 283;
1104 283;
927 282;
1095 282;
1102 282;
1167 282;
1015 281;
1269 281;
1354 281;
887 280;
965 280;
1013 280;
1017 280;
1022 280;
1059 280;
1087 280;
1243 280;
804 279;
866 279;
1123 279;
1157 279;
1357 279;
725 278;
798 278;
1253 278;
1275 278;
834 277;
1046 277;
1261 277;
1348 277;
585 276;
729 276;
760 276;
854 276;
901 276;
914 276;
970 276;
1150 276;
592 275;
793 275;
800 275;
837 275;
976 275;
977 275;
1032 275;
831 274;
835 274;
1027 274;
1147 274;
1249 274;
1253 274;
583 273;
814 273;
928 273;
998 273;
742 272;
744 272;
845 272;
867 272;
1040 272;
1260 272;
848 271;
915 271;
958 271;
988 271;
1017 271;
1144 271;
1244 271;
1281 271;
906 270;
1048 270;
1252 270;
584 269;
812 269;
923 269;
996 269;
1026 269;
1029 269;
579 268;
588 268;
930 268;
989 268;
1028 268;
1035 268;
1155 268;
1350 268;
1125 267;
1145 267;
1231 267;
566 266;
585 266;
606 266;
944 266;
952 266;
955 266;
988 266;
599 265;
822 265;
936 265;
1108 265;
1114 265;
1225 265;
821 264;
928 264;
939 264;
943 264;
967 264;
1045 264;
1115 264;
1244 264;
721 263;
968 263;
1011 263;
1094 263;
1113 263;
1138 263;
1248 263;
608 262;
609 262;
1010 262;
1018 262;
1024 262;
1101 262;
854 260;
871 260;
921 260;
929 260;
937 260;
961 260;
1089 260;
1151 260;
944 259;
1263 259;
815 258;
1276 258;
857 257;
876 257;
920 257;
1019 257;
1054 257;
1094 257;
1170 257;
1338 257;
833 256;
955 256;
969 256;
1045 256;
1119 256;
1237 256;
1248 256;
1264 256;
836 255;
871 255;
1014 255;
1088 255;
1170 255;
931 254;
942 254;
1237 254;
920 253;
955 253;
598 252;
923 252;
936 251;
939 251;
968 251;
1126 251;
969 250;
1130 250;
616 249;
1095 249;
1137 249;
1138 249;
1236 249;
1245 249;
1276 249;
1334 249;
1110 248;
834 247;
1160 247;
601 246;
1083 246;
1225 245;
1324 245;
837 244;
1121 244;
1142 244;
1218 244;
1094 243;
608 242;
1113 242;
1218 242;
1228 242;
1253 242;
1321 242;
1207 240;
1077 239;
1216 239;
1239 238;
1137 237;
838 236;
1217 236;
1265 236;
1130 235;
1245 235;
1317 235;
1043 234;
1070 234;
1220 234;
1260 234;
1316 234;
1068 233;
1198 233;
635 232;
1087 232;
624 231;
630 231;
844 231;
1076 231;
616 230;
622 230;
1093 230;
1313 230;
641 229;
1234 229;
1118 228;
1122 228;
1125 227;
1070 226;
1252 226;
1264 226;
828 225;
1076 225;
1071 224;
1074 224;
1094 224;
1196 224;
1076 223;
1122 223;
600 222;
603 222;
842 222;
603 221;
1076 221;
1108 221;
1125 221;
1240 221;
647 220;
1081 220;
1082 220;
1093 220;
843 219;
1088 219;
1101 218;
644 217;
1057 217;
830 216;
831 216;
1112 216;
1211 215;
1231 215;
1305 215;
1120 214;
1209 214;
1212 214;
634 213;
824 213;
1054 213;
830 212;
1096 212;
1129 212;
1248 212;
623 211;
1108 211;
1209 211;
658 210;
1055 210;
1079 210;
1128 210;
1249 210;
840 209;
1085 208;
1233 208;
1116 207;
1212 207;
1243 207;
622 206;
1090 206;
1219 206;
617 205;
620 205;
845 205;
662 204;
1107 204;
1238 204;
657 203;
1083 203;
1061 202;
1202 202;
1230 202;
1237 202;
1291 201;
680 200;
1069 200;
1242 200;
649 198;
674 198;
680 197;
1078 197;
1079 197;
1058 196;
1209 195;
841 194;
1102 194;
1037 193;
1071 193;
1231 193;
1243 193;
1033 192;
1280 192;
1034 191;
1085 191;
659 190;
1046 190;
824 189;
1061 189;
1279 189;
685 188;
1057 187;
1056 186;
1209 186;
686 185;
700 185;
1025 185;
1075 185;
1080 185;
1094 185;
654 184;
680 184;
686 184;
1047 184;
1053 184;
1076 184;
1090 184;
678 183;
1068 183;
1228 183;
702 182;
1214 182;
1054 180;
1216 180;
1232 180;
1093 179;
1197 179;
699 178;
1228 178;
1230 178;
1231 178;
1190 177;
1213 177;
1223 177;
686 176;
719 176;
1059 176;
685 175;
700 175;
1075 175;
1190 175;
662 174;
1192 174;
1263 174;
688 172;
708 172;
1197 172;
1064 171;
1063 170;
670 169;
690 169;
699 169;
1007 169;
1032 169;
1085 169;
1190 169;
703 168;
729 167;
985 167;
1195 167;
705 166;
982 166;
1212 166;
1182 165;
726 164;
1182 164;
1217 164;
683 163;
750 163;
1015 163;
1037 163;
1178 163;
1216 163;
689 162;
1015 161;
1173 161;
1212 161;
694 160;
737 160;
760 160;
986 160;
1009 160;
1044 160;
700 159;
765 159;
964 158;
1026 158;
1165 158;
960 157;
1044 157;
1174 157;
1180 157;
714 155;
737 155;
709 154;
724 154;
703 153;
984 153;
1014 153;
1187 153;
978 152;
993 152;
1037 152;
1179 152;
717 151;
954 151;
942 150;
1158 150;
1162 150;
775 149;
981 149;
719 148;
752 148;
1034 148;
717 147;
749 147;
765 147;
789 147;
1018 147;
1182 147;
730 146;
970 146;
995 146;
1131 146;
1158 146;
732 145;
924 145;
995 145;
1009 145;
716 144;
778 144;
1140 144;
713 143;
731 143;
848 143;
1102 143;
1148 143;
770 142;
980 142;
1097 142;
1144 142;
1159 142;
726 141;
755 141;
770 141;
780 141;
967 141;
1006 141;
1095 141;
1137 141;
1162 141;
778 140;
1131 140;
732 139;
827 139;
847 139;
1004 139;
1148 139;
1150 139;
737 138;
754 138;
824 138;
835 138;
920 138;
999 138;
1107 138;
1137 138;
806 137;
808 137;
826 137;
850 137;
1012 137;
1037 137;
753 136;
774 136;
900 136;
953 136;
957 136;
1162 136;
925 135;
997 135;
1108 135;
1122 135;
1157 135;
741 134;
856 134;
871 134;
903 134;
1031 134;
1100 134;
1144 134;
738 133;
776 133;
918 133;
962 133;
971 133;
1141 133;
762 132;
829 132;
951 132;
817 131;
944 131;
773 130;
861 130;
867 130;
969 130;
790 129;
816 129;
916 129;
925 129;
779 128;
836 128;
756 127;
806 127;
823 127;
987 127;
1130 127;
798 126;
900 126;
769 125;
856 125;
905 125;
959 125;
1112 125;
761 124;
835 124;
930 124;
1103 124;
1153 124;
867 123;
878 123;
880 123;
969 123;
1116 123;
831 122;
884 122;
777 121;
818 121;
984 120;
1098 120;
837 119;
856 119;
932 119;
977 119;
1125 119;
940 118;
954 118;
975 118;
1120 118;
1124 118;
1144 118;
923 117;
1109 117;
1116 117;
877 116;
1130 115;
1133 115;
769 114;
818 114;
855 114;
1128 114;
795 113;
892 113;
1115 113;
1120 113;
1133 113;
818 112;
897 112;
965 112;
1119 112;
750 111;
835 111;
940 111;
963 111;
1105 111;
850 110;
868 110;
982 110;
788 109;
833 109;
866 109;
882 109;
926 109;
894 108;
967 108;
1108 108;
1117 108;
785 107;
858 107;
877 107;
886 107;
932 107;
1109 107;
780 106;
1130 106;
771 105;
802 105;
923 105;
1115 105;
831 103;
1091 103;
1128 102;
858 101;
846 100;
898 100;
1084 100;
908 99;
912 99;
1116 99;
820 98;
825 98;
858 96;
1079 96;
896 95;
1077 95;
897 94;
1088 94;
920 93;
1076 93;
906 92;
1078 92;
1063 90;
1102 89;
1104 87;
1049 85;
1057 85;
1085 83;
1058 82;
1022 80;
1033 80;
1063 79;
1022 77;
1056 76;
1007 75;
1027 74;
1015 73;
1058 73;
998 69;
990 63;
905 62;
916 62;
925 62;
929 60;
908 58]


text:
type FdmAffineModelTermStructure{B <: BusinessCalendar, DC <: DayCount, A <: AffineModel} <: YieldTermStructure
  settlement_days::Int
  referenceDate::Date
  calendar::B
  dc::DC
  modelReferenceDate::Date
  model::A
  r::Vector{Float64}
  t::Float64
end

FdmAffineModelTermStructure{B <: BusinessCalendar, DC <: DayCount, A <: AffineModel}(referenceDate::Date, cal::B, dc::DC, modelReferenceDate::Date, model::A, r::Vector{Float64}) =
                            FdmAffineModelTermStructure{B, DC, A}(0, referenceDate, cal, dc, modelReferenceDate, model, r, year_fraction(dc, modelReferenceDate, referenceDate))

discount_impl(ts::FdmAffineModelTermStructure, T::Float64) = discount_bond(ts.model, ts.t, T + ts.t, ts.r)

set_variable(ts::FdmAffineModelTermStructure, r::Vector{Float64}) = ts.r = r

type FdmAffineModelSwapInnerValue{M1 <: Model, M2 <: Model, FM <: FdmMesher} <: FdmInnerValueCalculator
  disModel::M1
  fwdModel::M2
  swap::VanillaSwap
  exerciseDates::Dict{Float64, Date}
  mesher::FM
  direction::Int
  disTs::FdmAffineModelTermStructure
  fwdTs::FdmAffineModelTermStructure

  function FdmAffineModelSwapInnerValue{M1, M2, FM}(disModel::M1, fwdModel::M2, swap::VanillaSwap, exerciseDates::Dict{Float64, Date}, mesher::FM, direction::Int)
    idx = swap.iborIndex
    newIbor = IborIndex(idx.familyName, idx.tenor, idx.fixingDays, idx.currency, idx.fixingCalendar, idx.convention, idx.endOfMonth, idx.dc)
    newSwap = VanillaSwap(swap.swapT, swap.nominal, swap.fixedSchedule, swap.fixedRate, swap.fixedDayCount, newIbor, swap.spread, swap.floatSchedule,
                          swap.floatDayCount, swap.pricingEngine, swap.paymentConvention)

    return new{M1, M2, FM}(disModel, fwdModel, newSwap, exerciseDates, mesher, direction)
  end

  FdmAffineModelSwapInnerValue{M1, M2, FM}(disModel::M1, fwdModel::M2, swap::VanillaSwap, exerciseDates::Dict{Float64, Date}, mesher::FM, direction::Int,
                                              disTs:: FdmAffineModelTermStructure, fwdTs::FdmAffineModelTermStructure) = new{M1, M2, FM}(disModel, fwdModel, swap, exerciseDates, mesher, direction, disTs, fwdTs)
end

FdmAffineModelSwapInnerValue{M1 <: Model, M2 <: Model, FM <: FdmMesher}(disModel::M1, fwdModel::M2, swap::VanillaSwap, exerciseDates::Dict{Float64, Date}, mesher::FM, direction::Int) =
                            FdmAffineModelSwapInnerValue{M1, M2, FM}(disModel, fwdModel, swap, exerciseDates, mesher, direction)

FdmAffineModelSwapInnerValue{M1 <: Model, M2 <: Model, FM <: FdmMesher}(disModel::M1, fwdModel::M2, swap::VanillaSwap,
                            exerciseDates::Dict{Float64, Date}, mesher::FM, direction::Int, disTs::FdmAffineModelTermStructure, fwdTs::FdmAffineModelTermStructure) =
                            FdmAffineModelSwapInnerValue{M1, M2, FM}(disModel, fwdModel, swap, exerciseDates, mesher, direction, disTs, fwdTs)

function get_state(::G2, calc::FdmAffineModelSwapInnerValue, coords::Vector{Int}, ::Float64)
  retVal = Vector{Float64}(2)
  retVal[1] = get_location(calc.mesher, coords, calc.direction)
  retVal[2] = get_location(calc.mesher, coords, calc.direction + 1)

  return retVal
end

get_state(model::HullWhite, calc::FdmAffineModelSwapInnerValue, coords::Vector{Int}, t::Float64) = [short_rate(get_dynamics(model), t, get_location(calc.mesher, coords, calc.direction))]

function inner_value(calc::FdmAffineModelSwapInnerValue, coords::Vector{Int}, i::Int, t::Float64)
  iterExerciseDate = get(calc.exerciseDates, t, Date())
  disRate = get_state(calc.disModel, calc, coords, t)
  fwdRate = get_state(calc.fwdModel, calc, coords, t)

  if !isdefined(calc, :disTs) || iterExerciseDate != reference_date(calc.disTs)
    # Build new calc
    disc = calc.disModel.ts
    newDisTs = FdmAffineModelTermStructure(iterExerciseDate, disc.calendar, disc.dc, reference_date(disc), calc.disModel, disRate)
    fwd = calc.fwdModel.ts
    newFwdTs = FdmAffineModelTermStructure(iterExerciseDate, fwd.calendar, fwd.dc, reference_date(fwd), calc.fwdModel, fwdRate)
    # cloining swap
    newSwap = clone(calc.swap, calc.swap.pricingEngine, newFwdTs)

    # new calc
    calc = FdmAffineModelSwapInnerValue(calc.disModel, calc.fwdModel, newSwap, calc.exerciseDates, calc.mesher, calc.direction, newDisTs, newFwdTs)
  else
    # do some updating
    set_variable(calc.disTs, disRate)
    set_variable(calc.fwdTs, fwdRate)
  end

  npv = 0.0
  for j = 1:2
    for i = 1:length(calc.swap.legs[j].coupons)
      cf = calc.swap.legs[j].coupons[i]
      if isa(cf, Coupon)
        npv += accrual_start_date(cf) >= iterExerciseDate ? amount(cf) * discount(calc.disTs, date(cf)) : 0.0
      end
    end
    if j == 1
      npv *= -1.0
    end
  end
  if isa(calc.swap.swapT, Receiver)
    npv *= -1.0
  end
  return max(0.0, npv), calc
end

avg_inner_value(calc::FdmAffineModelSwapInnerValue, coords::Vector{Int}, i::Int, t::Float64) = inner_value(calc, coords, i, t)


text:
using Bits, Test
using Bits: NOTFOUND

x ≜ y = typeof(x) == typeof(y) && x == y

@testset "bitsize" begin
    for T in (Base.BitInteger_types..., Float16, Float32, Float64)
        @test bitsize(T) === sizeof(T) * 8
        @test bitsize(zero(T)) === bitsize(one(T)) === bitsize(T)
    end
    @test bitsize(BigInt) === Bits.INF
    @test bitsize(Bool) === 1
    @test bitsize(Float64) === 64
    @test bitsize(Float32) === 32
    @test bitsize(Float16) === 16
    @test_throws MethodError bitsize(BigFloat)
    @test bitsize(BigFloat(1, 256)) == 321
    @test bitsize(BigFloat(1, 100)) == 165
end

@testset "bit functions" begin
    @testset "weight" begin
        for T = (Base.BitInteger_types..., BigInt)
            @test weight(T(123)) === 6
            T == BigInt && continue
            @test weight(typemax(T)) === bitsize(T) - (T <: Signed)
            @test weight(typemin(T)) === Int(T <: Signed)
        end
        @test weight(big(-1)) === weight(big(-999)) === Bits.INF
    end

    @testset "bit & tstbit" begin
        val(::typeof(bit), x) = x
        val(::typeof(tstbit), x) = x % Bool
        for _bit = (bit, tstbit)
            for T = (Base.BitInteger_types..., BigInt)
                T0, T1 = T(0), T(1)
                @test _bit(T(0), rand(1:bitsize(T))) ≜ val(_bit, T0)
                @test _bit(T(1), 1) ≜ val(_bit, T1)
                @test _bit(T(1), 2) ≜ val(_bit, T0)
                @test _bit(T(5), 1) ≜ val(_bit, T1)
                @test _bit(T(5), 2) ≜ val(_bit, T0)
                @test _bit(T(5), 3) ≜ val(_bit, T1)
                @test _bit(T(5), 4) ≜ val(_bit, T0)
            end
            @test _bit( 1.0, 64) ≜ val(_bit, 0)
            @test _bit(-1.0, 64) ≜ val(_bit, 1)
            @test _bit( Float32(1.0), 32) ≜ val(_bit, Int32(0))
            @test _bit(-Float32(1.0), 32) ≜ val(_bit, Int32(1))
            x = BigFloat(-1.0, 128)
            for i=1:128+65
                @test _bit(x, i) == val(_bit, big(i ∈ [128, 129, 128+65]))
            end
        end
    end

    @testset "mask" begin
        for T = Base.BitInteger_types
            i = rand(0:min(999, bitsize(T)))
            @test mask(T, i) ≜ Bits.mask_2(T, i)
            @test mask(T) ≜ -1 % T
        end
        i = rand(0:bitsize(Bits.Word))
        @test mask(Bits.Word, i) === mask(i)
        @test mask(Bits.Word) === mask()
        @test mask(0) == 0
        @test mask(1) == 1
        @test mask(Sys.WORD_SIZE) === mask() === 0xffffffffffffffff
        @test mask(UInt64, 64) == mask(64)
        @test mask(UInt64, 63) == 0x7fffffffffffffff
        @test mask(BigInt, 0) ≜ big(0)
        @test mask(BigInt, 1) ≜ big(1)
        @test mask(BigInt, 1024) ≜ big(2)^1024 - 1

        # 2-arg mask
        for T = (Base.BitInteger_types..., BigInt)
            j, i = minmax(rand(0:min(999, bitsize(T)), 2)...)
            m = mask(T, j, i)
            @test m ≜ Bits.mask_2(T, j, i)
            if T === Bits.Word
                @test m === mask(j, i)
            end
            @test count_ones(m) == i-j
            @test m >>> j ≜ mask(T, i-j)
            @test mask(T, j, -1) ≜ ~mask(T, j)
        end
        @test mask(UInt64, 2, 4) === 0x000000000000000c
    end

    @testset "masked" begin
        for T = (Base.BitInteger_types..., BigInt)
            j, i = minmax(rand(0:min(999, bitsize(T)), 2)...)
            @test masked(mask(T), j) ≜ mask(T, j)
            @test masked(mask(T), j, i) ≜ mask(T, j, i)
            @test masked(mask(T, i), i) ≜ mask(T, i)
            @test masked(mask(T, i), j, i) ≜ mask(T, j, i)
            T == BigInt && continue
            x = rand(T)
            @test masked(x, j) ≜ x & mask(T, j)
            @test masked(x, j, i) ≜ x & mask(T, j, i)
        end
        @test masked(0b11110011, 1, 5) ≜ 0b00010010
        @test masked(-1.0, 52, 63) === 1.0
    end

    @testset "low0, low1, scan0, scan1" begin
        for T = (Base.BitInteger_types..., BigInt)
            x = T(0b01011010)
            @test low1(x, 0) == NOTFOUND
            @test low1(x)    == 2
            @test low1(x, 1) == 2
            @test low1(x, 2) == 4
            @test low1(x, 3) == 5
            @test low1(x, 4) == 7
            for i = 5:min(128, bitsize(T))+1
                @test low1(x, i) == NOTFOUND
            end
            @test low0(x, 0) == NOTFOUND
            @test low0(x)    == 1
            @test low0(x, 1) == 1
            @test low0(x, 2) == 3
            @test low0(x, 3) == 6
            @test low0(x, 4) == 8
            for i = 5:min(128, bitsize(T))
                @test low0(x, i) == (i+4 <= bitsize(T) ? i+4 : NOTFOUND)
            end

            @test scan1(x, 0) == NOTFOUND
            @test scan1(x, 1) == 2
            @test scan1(x)    == 2
            @test scan1(x, 2) == 2
            @test scan1(x, 3) == 4
            @test scan1(x, 4) == 4
            @test scan1(x, 5) == 5
            @test scan1(x, 6) == 7
            @test scan1(x, 7) == 7
            for i = 8:min(128, bitsize(T))+1
                @test scan1(x, i) == NOTFOUND
            end
            @test scan0(x, 0) == NOTFOUND
            @test scan0(x, 1) == 1
            @test scan0(x)    == 1
            @test scan0(x, 2) == 3
            @test scan0(x, 3) == 3
            @test scan0(x, 4) == 6
            @test scan0(x, 5) == 6
            @test scan0(x, 6) == 6
            @test scan0(x, 7) == 8
            @test scan0(x, 8) == 8
            for i = 9:min(128, bitsize(T))
                @test scan0(x, i) == i
            end
            T === BigInt && continue
            @test scan0(x, bitsize(T)+1) == NOTFOUND
        end
    end
end

@testset "bits" begin
    for T in Base.BitInteger_types
        v = bits(one(T))
        @test length(v) == bitsize(T)
        @test v[1] === true
        @test count(v) == 1
        v = bits(zero(T))
        @test length(v) == bitsize(T)
        @test v[1] === false
        @test count(v) == 0
    end
    v = bits(7)
    @test v[1] === v[2] === v[3] === true
    @test count(v) == 3
    for T = (Int64, UInt64)
        v = bits(T(2)^63)
        @test v[64] === true
        @test count(v) == 1
    end
    @test bits(Float64(-16))[53:end] == [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1]
    @test count(bits(Float64(-16))) == 4

    @testset "array indexing" begin
        v = bits(1234)[1:8]
        @test v == [0, 1, 0, 0, 1, 0, 1, 1]
        @test v isa Bits.BitVector1Mask
        v = bits(1.2)[10:17]
        @test v == [1, 0, 0, 1, 1, 0, 0, 1]
        @test v isa Bits.BitVector1Mask{Int64}
        @test all(bits(123)[[1, 2, 4, 5, 6, 7]])
        @test count(bits(123)) == 6
        # test optimization for v[i:j]
        for T = (Base.BitInteger_types..., BigInt)
            i, j = minmax(rand(1:min(999, bitsize(T)), 2)...)
            v = bits(rand(T == BigInt ? (big(-999):big(999)) : T))
            @test v[i:j] == v[i:1:j]
        end
    end
end


text:
#TODO:
# check_parent
# hash
# elem_type  & friends
#

mutable struct RealField <: Nemo.Field
  K::AnticNumberField
  P::InfPlc

  function RealField(K::AnticNumberField, P::InfPlc)
    if !isreal(P)
      error("place is not real")
    end
    r = new()
    r.K = K
    r.P = P
    return r
  end
end

@doc Markdown.doc"""
    real_field(K::AnticNumberField, i::Int)

The real field using the $i$-th conjugate for evaluation and comparison.
$i$ has to define a real embedding.
"""
function real_field(K::AnticNumberField, i::Int)
  r1, r2 = signature(K)
  if i > r1 || i < 1
    error("$i does not define a real embedding")
  end
  P = infinite_places(K)[i]
  return RealField(K, P)
end

@doc Markdown.doc"""
    real_field(K::AnticNumberField, P::InfPlc)

The real field using the real place $P$ to define the embedding for
evaluation and comparison.
"""
function real_field(K::AnticNumberField, P::InfPlc)
  return RealField(K, P)
end

function show(io::IO, K::RealField)
  println(io, "RealField represented by $(K.parent.pol) and $(K.P)")
end

mutable struct RealFieldElem <: Nemo.FieldElem
  data::nf_elem
  parent::RealField
  function RealFieldElem(K::RealField, a::nf_elem)
    r = new()
    r.data = a
    r.parent = K
    return r
  end
end

function show(io::IO, e::RealFieldElem)
  println(io, "$(e.data) = ", evaluate(e.data, e.parent.P))
end

@doc Markdown.doc"""
    evaluate(a::nf_elem, P::InfPlc, p::Int = 20)

The evaluation of $a$ at the place $P$, i.e. a real or complex value.
$p$ specifies the precision to be returned.
"""
function evaluate(a::nf_elem, P::InfPlc, p::Int = 10)
  return conjugates_arb(a, p)[P.i]
end

@doc Markdown.doc"""
    evaluate(a::RealFieldElem, p::Int = 10)

The value of $a$ as a real number under the chosen embedding.
$p$ specifies the prescision of the result.
"""
function evaluate(a::RealFieldElem, p::Int = 10)
  return evaluate(a.data, a.parent.P, p)
end

function Base.isless(a::RealFieldElem, b::RealFieldElem)
  return sign(a.data-b.data, b.parent.P) < 0
end

Base.isless(a::RealFieldElem, b::Integer) = isless(a, a.parent(b))
Base.isless(a::RealFieldElem, b::fmpz) = isless(a, a.parent(b))
Base.isless(a::RealFieldElem, b::fmpq) = isless(a, a.parent(b))
Base.isless(a::Integer, b::RealFieldElem) = isless(b.parent(a), b)
Base.isless(a::fmpz, b::RealFieldElem) = isless(b.parent(a), b)
Base.isless(a::fmpq, b::RealFieldElem) = isless(b.parent(a), b)

Base.isequal(a::RealFieldElem, b::RealFieldElem) = isequal(a.data, b.data)

parent(a::RealFieldElem) = a.parent

==(a::RealFieldElem, b::RealFieldElem) = a.data == b.data
+(a::RealFieldElem, b::RealFieldElem) = RealFieldElem(a.parent, a.data+b.data)
-(a::RealFieldElem, b::RealFieldElem) = RealFieldElem(a.parent, a.data-b.data)
*(a::RealFieldElem, b::RealFieldElem) = RealFieldElem(a.parent, a.data*b.data)
//(a::RealFieldElem, b::RealFieldElem) = RealFieldElem(a.parent, a.data//b.data)
-(a::RealFieldElem) = RealFieldElem(a.parent, -a.data)
inv(a::RealFieldElem) = RealFieldElem(a.parent, inv(a.data))
*(a::Integer, b::RealFieldElem) = RealFieldElem(b.parent, a*b.data)
*(a::fmpz, b::RealFieldElem) = RealFieldElem(b.parent, a*b.data)
*(a::fmpq, b::RealFieldElem) = RealFieldElem(b.parent, a*b.data)
+(a::Integer, b::RealFieldElem) = RealFieldElem(b.parent, a+b.data)
+(a::fmpz, b::RealFieldElem) = RealFieldElem(b.parent, a+b.data)
+(a::fmpq, b::RealFieldElem) = RealFieldElem(b.parent, a+b.data)

(K::RealField)(a::nf_elem) = RealFieldElem(K, a)
(K::RealField)(a::Integer) = RealFieldElem(K, K.parent(a))
(K::RealField)(a::fmpz) = RealFieldElem(K, K.parent(a))
(K::RealField)(a::fmpq) = RealFieldElem(K, K.parent(a))
(K::RealField)(a::RealFieldElem) = RealFieldElem(K, a.data)


text:
# setup internal paths -
_PATH_TO_SRC = dirname(pathof(@__MODULE__))

# load external packages that we depend upon -
using Reexport
@reexport using DataFrames
@reexport using CSV
using HTTP
using TOML
using Dates
using JSON

# load my codes -
include(joinpath(_PATH_TO_SRC,"Types.jl"))
include(joinpath(_PATH_TO_SRC,"Base.jl"))
include(joinpath(_PATH_TO_SRC,"Network.jl"))
include(joinpath(_PATH_TO_SRC,"Factory.jl"))
include(joinpath(_PATH_TO_SRC,"Handler.jl"))