text:
const mjNGROUP = 6
const mMAXOVERLAY = 500

const mjMAXLINE = 100
const mjMAXLINEPNT = 1000
const mjMAXPLANEGRID = 200

@cenum mjtCatBit::Cint begin  # bitflags for vGeom category
    mjCAT_STATIC = 1    # model elements in body 0
    mjCAT_DYNAMIC = 2    # model elements in all other bodies
    mjCAT_DECOR = 4    # decorative geoms
    mjCAT_ALL = 7    # select all categories
end

@cenum mjtMouse::Cint begin   # mouse interaction mode
    mjMOUSE_NONE = 0   # no action
    mjMOUSE_ROTATE_V           # rotate, vertical plane
    mjMOUSE_ROTATE_H           # rotate, horizontal plane
    mjMOUSE_MOVE_V             # move, vertical plane
    mjMOUSE_MOVE_H             # move, horizontal plane
    mjMOUSE_ZOOM               # zoom
    mjMOUSE_SELECT             # selection
end

@cenum mjtPertBit::Cint begin # mouse perturbations
    mjPERT_TRANSLATE = 1    # translation
    mjPERT_ROTATE = 2    # rotation
end

@cenum mjtCamera::Cint begin  # abstract camera type
    mjCAMERA_FREE = 0   # free camera
    mjCAMERA_TRACKING          # tracking camera; uses trackbodyid
    mjCAMERA_FIXED             # fixed camera; uses fixedcamid
    mjCAMERA_USER              # user is responsible for setting OpenGL camera
end

@cenum mjtLabel::Cint begin   # object labeling
    mjLABEL_NONE = 0    # nothing
    mjLABEL_BODY               # body labels
    mjLABEL_JOINT              # joint labels
    mjLABEL_GEOM               # geom labels
    mjLABEL_SITE               # site labels
    mjLABEL_CAMERA             # camera labels
    mjLABEL_LIGHT              # light labels
    mjLABEL_TENDON             # tendon labels
    mjLABEL_ACTUATOR           # actuator labels
    mjLABEL_CONSTRAINT         # constraint labels
    mjLABEL_SKIN               # skin labels
    mjLABEL_SELECTION          # selected object
    mjLABEL_SELPNT             # coordinates of selection point
    mjLABEL_CONTACTFORCE       # magnitude of contact force

    mjNLABEL                   # number of label types
end

@cenum mjtFrame::Cint begin   # frame visualization
    mjFRAME_NONE = 0    # no frames
    mjFRAME_BODY               # body frames
    mjFRAME_GEOM               # geom frames
    mjFRAME_SITE               # site frames
    mjFRAME_CAMERA             # camera frames
    mjFRAME_LIGHT              # light frames
    mjFRAME_WORLD              # world frame

    mjNFRAME                   # number of visualization frames
end

@cenum mjtVisFlag::Cint begin # flags enabling model element visualization
    mjVIS_CONVEXHULL = 0    # mesh convex hull
    mjVIS_TEXTURE              # textures
    mjVIS_JOINT                # joints
    mjVIS_ACTUATOR             # actuators
    mjVIS_CAMERA               # cameras
    mjVIS_LIGHT                # lights
    mjVIS_TENDON               # tendons
    mjVIS_RANGEFINDER          # rangefinder sensors
    mjVIS_CONSTRAINT           # point constraints
    mjVIS_INERTIA              # equivalent inertia boxes
    mjVIS_SCLINERTIA           # scale equivalent inertia boxes with mass
    mjVIS_PERTFORCE            # perturbation force
    mjVIS_PERTOBJ              # perturbation object
    mjVIS_CONTACTPOINT         # contact points
    mjVIS_CONTACTFORCE         # contact force
    mjVIS_CONTACTSPLIT         # split contact force into normal and tanget
    mjVIS_TRANSPARENT          # make dynamic geoms more transparent
    mjVIS_AUTOCONNECT          # auto connect joints and body coms
    mjVIS_COM                  # center of mass
    mjVIS_SELECT               # selection point
    mjVIS_STATIC               # static bodies
    mjVIS_SKIN                 # skin

    mjNVISFLAG                 # number of visualization flags
end

@cenum mjtRndFlag::Cint begin # flags enabling rendering effects
    mjRND_SHADOW = 0    # shadows
    mjRND_WIREFRAME            # wireframe
    mjRND_REFLECTION           # reflections
    mjRND_ADDITIVE             # additive transparency
    mjRND_SKYBOX               # skybox
    mjRND_FOG                  # fog
    mjRND_HAZE                 # haze
    mjRND_SEGMENT              # segmentation with random color
    mjRND_IDCOLOR              # segmentation with segid color

    mjNRNDFLAG                 # number of rendering flags
end

@cenum mjtStereo::Cint begin  # type of stereo rendering
    mjSTEREO_NONE = 0    # no stereo; use left eye only
    mjSTEREO_QUADBUFFERED      # quad buffered; revert to side-by-side if no hardware support
    mjSTEREO_SIDEBYSIDE        # side-by-side
end

@uninitable mutable struct mjvPerturb
    select::Cint
    skinselect::Cint
    active::Cint
    refpos::SVector{3,mjtNum}
    refquat::SVector{4,mjtNum}
    localpos::SVector{3,mjtNum}
    scale::mjtNum
end

@uninitable mutable struct mjvCamera
    type::Cint
    fixedcamid::Cint
    trackbodyid::Cint
    lookat::SVector{3,mjtNum}
    distance::mjtNum
    azimuth::mjtNum
    elevation::mjtNum
end

@uninitable struct mjvGLCamera
    pos::SVector{3,Cfloat}
    forward::SVector{3,Cfloat}
    up::SVector{3,Cfloat}
    frustum_center::Cfloat
    frustum_bottom::Cfloat
    frustum_top::Cfloat
    frustum_near::Cfloat
    frustum_far::Cfloat
end

@uninitable struct mjvGeom
    type::Cint
    dataid::Cint
    objtype::Cint
    objid::Cint
    category::Cint
    texid::Cint
    texuniform::Cint
    texcoord::Cint
    segid::Cint
    texrepeat::SVector{2,Cfloat}
    size::SVector{3,Cfloat}
    pos::SVector{3,Cfloat}
    mat::SVector{9,Cfloat}
    rgba::SVector{4,Cfloat}
    emission::Cfloat
    specular::Cfloat
    shininess::Cfloat
    reflectance::Cfloat
    label::SVector{100,UInt8}
    camdist::Cfloat
    rbound::Cfloat
    transparent::mjtByte
end

@uninitable struct mjvLight
    pos::SVector{3,Cfloat}
    dir::SVector{3,Cfloat}
    attenuation::SVector{3,Cfloat}
    cutoff::Cfloat
    exponent::Cfloat
    ambient::SVector{3,Cfloat}
    diffuse::SVector{3,Cfloat}
    specular::SVector{3,Cfloat}
    headlight::mjtByte
    directional::mjtByte
    castshadow::mjtByte
end

@uninitable mutable struct mjvOption
    label::Cint
    frame::Cint
    geomgroup::SVector{mjNGROUP,mjtByte}
    sitegroup::SVector{mjNGROUP,mjtByte}
    jointgroup::SVector{mjNGROUP,mjtByte}
    tendongroup::SVector{mjNGROUP,mjtByte}
    actuatorgroup::SVector{mjNGROUP,mjtByte}
    flags::SVector{Int(mjNVISFLAG),mjtByte}
end

@uninitable mutable struct mjvScene
    maxgeom::Cint
    ngeom::Cint
    geoms::Ptr{mjvGeom}
    geomorder::Ptr{Cint}

    nskin::Cint
    skinfacenum::Ptr{Cint}
    skinvertadr::Ptr{Cint}
    skinvertnum::Ptr{Cint}
    skinvert::Ptr{Cfloat}
    skinnormal::Ptr{Cfloat}

    nlight::Cint
    lights::SVector{8,mjvLight}
    camera::SVector{2,mjvGLCamera}
    enabletransform::mjtByte
    translate::SVector{3,Cfloat}
    rotate::SVector{4,Cfloat}
    scale::Cfloat
    stereo::Cint
    flags::SVector{Int(mjNRNDFLAG),mjtByte}
end

@uninitable mutable struct mjvFigure
    flg_legend::Cint
    flg_ticklabel::SVector{2,Cint}
    flg_extend::Cint
    flg_barplot::Cint
    flg_selection::Cint
    flg_symmetric::Cint

    legendoff::Cint
    gridsize::SVector{2,Cint}
    selection::Cint
    highlight::SVector{2,Cint}
    gridrgb::SVector{3,Cfloat}
    gridwidth::Cfloat
    figurergba::SVector{4,Cfloat}
    panergba::SVector{4,Cfloat}
    legendrgba::SVector{4,Cfloat}
    textrgb::SVector{3,Cfloat}
    range::SVector{2,SVector{2,Cfloat}}
    xlabel::SVector{100,UInt8}
    title::SVector{100,UInt8}
    xformat::SVector{20,UInt8}
    yformat::SVector{20,UInt8}
    minwidth::SVector{20,UInt8}

    linepnt::SVector{mjMAXLINE,Cint}
    linergb::SVector{mjMAXLINE,SVector{3,Cfloat}}
    linewidth::SVector{mjMAXLINE,Cfloat}
    linedata::SVector{mjMAXLINE,SVector{2 * mjMAXLINEPNT,Cfloat}}
    linename::SVector{mjMAXLINE,SVector{100,UInt8}}

    xaxispixel::SVector{2,Cint}
    yaxispixel::SVector{2,Cint}
    xaxisdata::SVector{2,Cint}
    yaxisdata::SVector{2,Cint}
end


text:
struct DictPairSet{K,V,T<:AbstractDict{K,V}} <: AbstractSet{Tuple{K,V}}
    dict::T
end
Base.length(x::DictPairSet) = length(x.dict)
Base.iterate(x::DictPairSet) =
    (r = iterate(x.dict); r === nothing ? nothing : (Tuple(r[1]), r[2]))
Base.iterate(x::DictPairSet, st) =
    (r = iterate(x.dict, st); r === nothing ? nothing : (Tuple(r[1]), r[2]))
Base.in(v::Pair, x::DictPairSet) = v in x.dict
Base.in(v::Tuple{Any,Any}, x::DictPairSet) = Pair(v[1], v[2]) in x.dict

pyjldict_iter(x::AbstractDict) = Py(Iterator(keys(x)))

pyjldict_contains(x::AbstractDict, k::Py) = Py(haskey(x, @pyconvert(keytype(x), k, return Py(false))))

pyjldict_clear(x::AbstractDict) = (empty!(x); Py(nothing))

pyjldict_getitem(x::AbstractDict, k::Py) = Py(x[pyconvert(keytype(x), k)])

pyjldict_setitem(x::AbstractDict, k::Py, v::Py) = (x[pyconvertarg(keytype(x), k, "key")] = pyconvertarg(valtype(x), v, "value"); Py(nothing))

pyjldict_delitem(x::AbstractDict, k::Py) = (delete!(x, pyconvert(keytype(x), k)); Py(nothing))

function pyjldict_update(x::AbstractDict, items_::Py)
    for item_ in items_
        (k, v) = pyconvert_and_del(Tuple{keytype(x), valtype(x)}, item_)
        x[k] = v
    end
    Py(nothing)
end

const pyjldicttype = pynew()

function init_jlwrap_dict()
    jl = pyjuliacallmodule
    filename = "$(@__FILE__):$(1+@__LINE__)"
    pybuiltins.exec(pybuiltins.compile("""
    class DictValue(AnyValue):
        __slots__ = ()
        __module__ = "juliacall"
        _jl_undefined_ = object()
        def __iter__(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjldict_iter)))
        def __contains__(self, key):
            return self._jl_callmethod($(pyjl_methodnum(pyjldict_contains)), key)
        def __getitem__(self, key):
            if key in self:
                return self._jl_callmethod($(pyjl_methodnum(pyjldict_getitem)), key)
            else:
                raise KeyError(key)
        def __setitem__(self, key, value):
            return self._jl_callmethod($(pyjl_methodnum(pyjldict_setitem)), key, value)
        def __delitem__(self, key):
            if key in self:
                return self._jl_callmethod($(pyjl_methodnum(pyjldict_delitem)), key)
            else:
                raise KeyError(key)
        def keys(self):
            return self._jl_callmethod($(pyjl_methodnum(Py ∘ keys)))
        def values(self):
            return self._jl_callmethod($(pyjl_methodnum(Py ∘ values)))
        def items(self):
            return self._jl_callmethod($(pyjl_methodnum(Py ∘ DictPairSet)))
        def get(self, key, default=None):
            if key in self:
                return self[key]
            else:
                return default
        def setdefault(self, key, default=None):
            if key not in self:
                self[key] = default
            return self[key]
        def clear(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjldict_clear)))
        def pop(self, key, default=_jl_undefined_):
            if key in self:
                ans = self[key]
                del self[key]
                return ans
            elif default is self._jl_undefined_:
                raise KeyError(key)
            else:
                return default
        def popitem(self):
            if len(self):
                return self._jl_callmethod($(pyjl_methodnum(Py ∘ pop!)))
            else:
                raise KeyError()
        def update(self, other=_jl_undefined_, **kwargs):
            if other is self._jl_undefined_:
                pass
            else:
                if hasattr(other, "keys"):
                    items = ((k, other[k]) for k in other.keys())
                else:
                    items = other
                self._jl_callmethod($(pyjl_methodnum(pyjldict_update)), items)
            if kwargs:
                self.update(kwargs)
    import collections.abc
    collections.abc.MutableMapping.register(DictValue)
    del collections
    """, filename, "exec"), jl.__dict__)
    pycopy!(pyjldicttype, jl.DictValue)
end

pyjl(v::AbstractDict) = pyjl(pyjldicttype, v)

ext: jl


text:
#***********************************************************************
#  Copyright 2006 John A. Trangenstein
#
#  This software is made available for research and instructional use
#  only.
#  You may copy and use this software without charge for these
#  non-commercial purposes, provided that the copyright notice and
#  associated text is reproduced on all copies.
#  For all other uses (including distribution of modified versions),
#  please contact the author at
#    John A. Trangenstein
#    Department of Mathematics
#    Duke University
#    Durham, NC 27708-0320
#    USA
#  or
#    johnt@math.duke.edu
#
#  This software is made available "as is" without any assurance that it
#  is completely correct, or that it will work for your purposes.
#  Use the software at your own risk.
#***********************************************************************
__precompile__()
using OffsetArrays
using DocOpt

include("consts.jl")

using consts

include("riemprob.jl")
include("linearad.jl")
include("consdiff.jl")
include("upwind.jl")

function do_computation(ncells, nsteps, verbose, print_solution)

    u    = OffsetArray(Float64, -2:ncells+1)
    x    = OffsetArray(Float64,  0:ncells)
    flux = OffsetArray(Float64,  0:ncells)
    dfdu = OffsetArray(Float64, -2:ncells+1)

    const tmax     =  0.8
    const cfl      =  0.9

    # array bounds:
    const fc=-2
    const lc=ncells+1
    const fm=0
    const lm=ncells-1
    const fs=0
    const ls=ncells-1
    const ifirst=0
    const ilast=ncells-1

    initsl(ncells,fc,lc,fm,lm,ifirst,ilast, u,x)

    if verbose
        println("x : $x u : $u")
    end

    bcmesh(fm,lm,ncells, x)
    bccells(fc,lc,ncells, u)
    fluxderv(fc,lc,fc,lc, u, dfdu)

    istep=0
    t=0.0
    while istep < nsteps && t < tmax
        bccells(fc,lc,ncells, u)
        fluxderv(fc,lc,fc,lc, u, dfdu)
        dt=cfl*stabledt(fc,lc,fm,lm,ifirst,ilast, u, dfdu,x)
        method(dt,fc,lc,fm,lm,fs,ls,ifirst,ilast, u, flux)
        consdiff(fc,lc,fs,ls,fm,lm,ifirst,ilast, x,flux, u)

        t=t+dt
        istep=istep+1
    end

    # write final results (plot later)
    if print_solution
        @unsafe for ic=0:ncells-1
            xc = (x[ic]+x[ic+1])*0.5
            uc = u[ic]
            @printf("%e %e\n",xc,uc)
        end
    end

end

function main()
    const script_name = basename(Base.source_path())
    const doc = """$script_name

Usage:
  $script_name -h | --help
  $script_name [-v | --verbose] [-p | --print] [--cells=<cells>] [--steps=<steps>] [--runs=<runs>]

Options:
  -h --help                  Show this screen.
  -v --verbose               Adds verbosity.
  -p --print                 Print solution.
  --cells=<cells>            Specify a number of cells [default: 100].
  --steps=<steps>            Specify a number of steps [default: 10000].
  --runs=<runs>              Specify a number of runs  [default: 2].
"""

    arguments = docopt(doc)

    verbose   = arguments["--verbose"]
    print_sol = arguments["--print"]
    ncells    = parse(Int, arguments["--cells"])
    nsteps    = parse(Int, arguments["--steps"])
    nruns     = parse(Int, arguments["--runs"])

    if verbose
        println("ncells: $ncells")
        println("nsteps: $nsteps")
    end

    for r = 1:nruns
        @time do_computation(ncells, nsteps, verbose, print_sol)
    end

end # main

main()


text:
const MP = MultivariatePolynomials
@testset "Polynomial" begin
    Mod.@polyvar x

    @test terms(polynomial([1, x^2, x, 2x^2])) == [3x^2, x, 1]
    @test terms(polynomial([x, 3x^4, 2], MP.UniqState())) == [3x^4, x, 2]
    @test terms(polynomial([x^3, 2x^3, x^2, -2x^2, x^2, x, 2, -2], MP.SortedState())) == [3x^3, x]

    @test polynomial(1 + x) == 1 + x
    @test leadingterm(1 + x) == x
    @test leadingterm(x - x) == 0
    @test leadingmonomial(x - x) == 1
    @test leadingcoefficient(x - x) == 0
    @test one(1 + x) == one(1.0 + x) == 1
    @test zero(1 + x) == zero(1.0 + x) == 0
    @test 1 != 1 + x
    @test 2x == x + x^2 + x - x^2
    @test x + x^2 - x^2 - x != 2x
    @test x^2 + x != x^2 + x + 1
    #@inferred one(1 + x)
    @inferred zero(1 + x)
    #@inferred one(1.0 + x)
    @inferred zero(1.0 + x)

    @test polynomial(1:2, [x, x^2]) == x + 2x^2
    @test polynomial(1:2, monomials(x, 1:2)) == 2x + x^2

    @test terms(polynomial(1 + x + x^2 - x + x^2)) == [2x^2, 1]
    @test terms(CustomPoly(1 + x + x^2 - x + x^2)) == [2x^2, 1]

    @test (1.0 + x) * x == x^2 + x
    @test constantterm(1, x) * (1 - x) == 1 - x
    @test promote_type(typeof(1-x), typeof(x)) <: AbstractPolynomial{Int}
    @test x != 1 - x

    @test term(x + x^2 - x) isa AbstractTerm
    @test term(x + x^2 - x) == x^2
    @test term(x - x) isa AbstractTerm
    @test iszero(term(x - x))
    @test_throws ErrorException term(x + x^2)

    Mod.@polyvar y

    p = 3x^2*y^4 + 2x
    @test +(p) === p
    @test *(p) === p
    @test terms(p)[end] == 2x
    typetests(p)
    typetests([p, x + y])
    @test (@inferred polynomial(p)) isa AbstractPolynomial{Int}
    @test (@inferred polynomial(p, Float64)) isa AbstractPolynomial{Float64}

    @test coefficient(2x + 4y^2 + 3, y^2) == 4
    @test coefficient(2x + 4y^2 + 3, x^2) == 0

    Mod.@polyvar a b
    @test coefficient((2a + b)x^2 + 2y^2 + 3, x^2, (x,y)) == 2a + b
    @test coefficient((2a + b)x^2 + 2y^2 + 3x^2*y, x^2, (x,)) == 2a + b + 3y

    @test (@inferred 2x^2*y + 0.0x*y) == 2x^2*y
    @test (@inferred 0.0x^2*y + 3x*y) == 3x*y

    @test iszero(((x + x) - 2x) * (x * (x ^ 2 + y ^ 2)))

    @test Tuple(variables([x + 1, y^2])) == (x, y)
    @test Tuple(variables([y^2, x + 1])) == (x, y)

    @test maxdegree(x^2 - x^2) == 0
    @test maxdegree(x^2 - x^2, x) == 0
    @test maxdegree(x^2 - x^2, y) == 0
    @test mindegree(x^2 - x^2) == 0
    @test mindegree(x^2 - x^2, x) == 0
    @test mindegree(x^2 - x^2, y) == 0
    @test extdegree(x^2 - x^2) == (0, 0)
    @test extdegree(x^2 - x^2, x) == (0, 0)
    @test extdegree(x^2 - x^2, y) == (0, 0)
    @test maxdegree(x*y + 2 + x^2*y + x + y) == 3
    @test maxdegree(x*y + 2 + x^2*y + x + y, x) == 2
    @test maxdegree(x*y + 2 + x^2*y + x + y, y) == 1
    @test mindegree(x*y + 2 + x^2*y + x + y) == 0
    @test mindegree(x*y + 2 + x^2*y + x + y, x) == 0
    @test mindegree(x*y + 2 + x^2*y + x + y, y) == 0
    @test extdegree(x*y + 2 + x^2*y + x + y) == (0, 3)
    @test extdegree(x*y + 2 + x^2*y + x + y, x) == (0, 2)
    @test extdegree(x*y + 2 + x^2*y + x + y, y) == (0, 1)
    @test extdegree(x*y + x^2*y, x) == (1, 2)
    @test extdegree(x*y + x^2*y, y) == (1, 1)
    @test leadingterm(x*y + 2 + x^2*y + x + y) == x^2*y
    @test nvariables(x + y - x) == 2
    @test nvariables(x + x^2) == 1

    @test coefficients(x*y + 2 + 3x^2*y + 4x + 6y, [x, x*y^2, x*y, x^2*y, y, x^3]) == [4, 0, 1, 3, 6, 0]

    # Doc examples
    @test coefficients(4x^2*y + x*y + 2x) == [4, 1, 2]
    @test coefficients(4x^2*y + x*y + 2x + 3, [x, 1, x*y, y]) == [2, 3, 1, 0]

    p = polynomial([4, 9], [x, x*x])
    @test coefficients(p) == [9, 4]
    @test monomials(p)[1] == x^2
    @test monomials(p)[2] == x
    @test p == dot([4, 9], [x, x*x])

    @inferred polynomial(i -> float(i), [x, x*x])
    @inferred polynomial(i -> 3 - float(i), monovec([x*x, x]))
    for p in (polynomial(i -> float(i), [x, x*x]),
              polynomial(i -> 3 - float(i), monovec([x*x, x])))
        @test coefficients(p) == [2.0, 1.0]
        @test monomials(p) == monovec([x^2, x])
    end

    @test (x + y)' == x + y
    @test transpose(x + y) == x + y
    @test transpose([1 2; 3 4] * x) == [1 3; 2 4] * x

    @test removemonomials(4x^2*y + x*y + 2x, [x*y]) == 4x^2*y + 2x

    @test_throws InexactError push!([1], x+1)

    @test polynomial([1 2; 3 4], [x^2, y]) == x^4 + 5x^2*y + 4y^2
    @test polynomial([1 2; 3 4], [x^2, y], Float64) isa AbstractPolynomial{Float64}
    @test polynomial([1 2; 3 4], [y, x^2]) == y^2 + 5x^2*y + 4x^4
    @test polynomial([1 2; 3 4], [y, x^2], Float64) isa AbstractPolynomial{Float64}
    @test polynomial([1 2; 3 4], monovec([y, x^2])) == x^4 + 5x^2*y + 4y^2
    @test polynomial([1 2; 3 4], monovec([y, x^2]), Float64) isa AbstractPolynomial{Float64}

    @test (@inferred round(2.6x + 1.001x^2)) == 3x + 1x^2
    @test (@inferred round(3.1x*y)) == 3x*y
    @test (@inferred round(2.613x + 1.1051x^2, digits=2)) ≈ 2.61x + 1.11x^2
    @test (@inferred round(3.145x*y, digits=1)) ≈ 3.1x*y

    @testset "Graded Lex Order" begin
        Mod.@polyvar x y z
        p = 3*y^2 + 2*y*x
        @test coefficients(p) == [2, 3]
        @test monomials(p) == monovec([x*y, y^2])
        # Examples from p. 59 of the 4th edition of "Ideals, Varieties, and Algorithms" of Cox, Little and O'Shea
        f = 4*x*y^2*z + 4*z^2 - 5*x^3 + 7*x^2*z^2
        @test coefficients(f) == [7, 4, -5, 4]
        @test monomials(f) == monovec([x^2*z^2, x*y^2*z, x^3, z^2])
    end

    @testset "Convertion" begin
        p = 2.5x + 1 - 2.5x
        @test convert(Int, p) == 1
        @test convert(typeof(p), p) === p
        @test convert(Union{Nothing, typeof(p)}, p) === p
    end

    @testset "Vector" begin
        Mod.@polyvar x y
        v = [x - 1, y + 1]
        @test nvariables(v) == 2
        @test variables(v)[1] == x
        @test variables(v)[2] == y
    end

    @testset "Effective variables" begin
        Mod.@polyvar x y z
        T = variable_union_type(x)
        @test x isa T
        @test y isa T
        @test z isa T
        @test T[x] == @inferred effective_variables(x + y - y)
        @test T[x, y] == @inferred effective_variables(z + x + y - z)
        @test T[y, z] == @inferred effective_variables(z + 0 * x + y)
        @test T[z] == @inferred effective_variables(z + 0 * x + y^0)
    end

    @testset "Complex" begin
        # See https://github.com/JuliaAlgebra/MultivariatePolynomials.jl/issues/128
        Mod.@polyvar x
        p = im * x + 2im * x^2
        @test p * p == -x^2 - 4x^4 - 4x^3
    end
end


text:
module LossFunctions

using Markdown
using Base.Cartesian
using InteractiveUtils: subtypes
using RecipesBase

import Base: *

# aggregation mode
include("aggmode.jl")

# observation dimension
include("obsdim.jl")

# trait functions
include("traits.jl")

# supervised losses
include("supervised.jl")

# IO and plot recipes
include("printing.jl")
include("plotrecipes.jl")

export
    # trait functions
    Loss,
    SupervisedLoss,
    MarginLoss,
    DistanceLoss,
    value, deriv, deriv2,
    isdistancebased, ismarginbased,
    isminimizable, isdifferentiable,
    istwicedifferentiable,
    isconvex, isstrictlyconvex,
    isstronglyconvex, isnemitski,
    isunivfishercons, isfishercons,
    islipschitzcont, islocallylipschitzcont,
    isclipable, isclasscalibrated, issymmetric,

    # relevant submodules
    AggMode, ObsDim,

    # margin-based losses
    ZeroOneLoss,
    LogitMarginLoss,
    PerceptronLoss,
    HingeLoss,
    L1HingeLoss,
    L2HingeLoss,
    SmoothedL1HingeLoss,
    ModifiedHuberLoss,
    L2MarginLoss,
    ExpLoss,
    SigmoidLoss,
    DWDMarginLoss,

    # distance-based losses
    LPDistLoss,
    L1DistLoss,
    L2DistLoss,
    PeriodicLoss,
    HuberLoss,
    EpsilonInsLoss,
    L1EpsilonInsLoss,
    L2EpsilonInsLoss,
    LogitDistLoss,
    QuantileLoss,
    LogCoshLoss,

    # other losses
    MisclassLoss,
    PoissonLoss,
    CrossEntropyLoss,

    # meta losses
    ScaledLoss,
    OrdinalMarginLoss,
    WeightedMarginLoss

end # module


text:
@views function StokesSolver!( f::Fields2D, params::ModelParameters, dom::ModelDomain, mat::MaterialParameters, BC::BoundaryConditions, Kuu::SparseMatrixCSC{Float64, Int64}, Kup::SparseMatrixCSC{Float64, Int64}, Kpu::SparseMatrixCSC{Float64, Int64}) #; dom=0, BC=0, phase_perc=0, nphase=0
    """ Solver """
    if params.JFNK==1
        println("Automatically switching to KSP GCR solver since JFNK is requested")
        solver = 2
    end
    # Get input residual
    fu = zeros(Float64, length(f.Fx) + length(f.Fy))
    fp = zeros(Float64, length(f.Fp))
    fu[1:length(f.Fx)]     .= f.Fx[:]
    fu[length(f.Fx)+1:end] .= f.Fy[:]
    fp .= f.Fp[:]
    #------------------------ PINNED PRESSURE
    if params.solver == -1
        # Pinned pressure case
        # Slightly compressible pressure block
        print("Pinned pressure")
        cP   = zeros(Float64,size(f.etac,1), size(f.etac,2))
        cP[1,1] = 1.0
        Kpu[1,:] .= 0.0
        f.Pc[1,1] = 0.0
        I    = f.NumP[:]
        J    = f.NumP[:]
        V    = cP[:]
        PP   = sparse( I, J, V)
        Kmat = [Kuu Kup; Kpu PP]
        F    = [fu; fp]
        dX   = Kmat\F
        println(size(f.Vx))
        println(size(f.NumVx))
        dVx = dX[ f.NumVx ]
        println( size(f.Vx[:,2:end-1]) )
        f.Vx[:,2:end-1] .= f.Vx[:,2:end-1] .+ dX[f.NumVx]
        f.Vy[2:end-1,:] .= f.Vy[2:end-1,:] .+ dX[f.NumVy]
        f.Pc            .= f.Pc            .+ dX[f.NumP.+maximum(f.NumVy)]
        f.Pc            .-= mean(f.Pc)
    elseif params.solver == 0
        #------------------------ PINNED PRESSURE
        # Slightly compressible pressure block
        if params.comp==1
            print("The solver accounts for physical compressibility")
            cP   = 1.0./(f.Kc.*params.dt).*ones(Float64,size(f.etac,1), size(f.etac,2))
        else
            print("The solver does not account for physical compressibility (artificial only)")
            cP   = 1.0./params.gamma.*ones(Float64,size(f.etac,1), size(f.etac,2))
        end
        I    = f.NumP[:]
        J    = f.NumP[:]
        V    = cP[:]
        PP   = sparse( I, J, V)
        Kmat = [Kuu Kup; Kpu PP]
        F    = [fu; fp]
        dX   = Kmat\F
        println(size(f.Vx))
        println(size(f.NumVx))
        dVx = dX[ f.NumVx ]
        println( size(f.Vx[:,2:end-1]) )
        f.Vx[:,2:end-1] .= f.Vx[:,2:end-1] .+ dX[f.NumVx]
        f.Vy[2:end-1,:] .= f.Vy[2:end-1,:] .+ dX[f.NumVy]
        f.Pc            .= f.Pc            .+ dX[f.NumP.+maximum(f.NumVy)]
        if params.comp==0 f.Pc            .-= mean(f.Pc) end
    elseif params.solver == 1
        DecoupledSolver!( f, params, fu, fp, Kuu, Kup, Kpu )
    elseif params.solver == 2
        nVx = length(f.Fx)
        nVy = length(f.Fy)
        nP  = length(f.Fp)
        if params.comp==1
            @tturbo coef      = 1.0./(f.Kc[:].*params.dt).*ones(length(f.etac))#.*etac[:]
            @tturbo coef_inv  = f.Kc[:].*params.dt.*ones(length(f.etac))#.*etac[:]
        else
            @tturbo coef      =    0.0*ones(length(f.etac))#.*etac[:]
            @tturbo coef_inv  = params.gamma*ones(length(f.etac))#.*etac[:]
        end
        @tturbo Kpp   = spdiagm(coef)
        @tturbo Kppi  = spdiagm(coef_inv)
        @tturbo M     = [Kuu Kup; Kpu Kpp]
        dx    = zeros(Float64, size(M,1))
        x     = zero(dx)
        Vxi   = f.Vx[:,2:end-1]
        Vyi   = f.Vy[2:end-1,:]
        x[1:nVx]         .= Vxi[:]
        x[nVx+1:nVx+nVy] .= Vyi[:]
        x[nVx+nVy+1:end] .= f.Pc[:]
        b     = zero(dx)
        b[1:nVx]         .= f.Fx[:]
        b[nVx+1:nVx+nVy] .= f.Fy[:]
        b[nVx+nVy+1:end] .= f.Fp[:]
        epsi   = 1e-8
        noisy  = 2
        # ACHTUNG: here copy solution arrays since they are overwitten in the JFNK solver
        Vx, Vy, Pc = copy(f.Vx), copy(f.Vy), copy(f.Pc)
        if params.JFNK == 0
            KSP_GCR_Stokes!( dx, M, b, epsi, noisy, Kuu, Kup, Kpu, Kpp, Kppi )
        else
            KSP_GCR_Stokes2!( dx, M, x, b, epsi, noisy, Kuu, Kup, Kpu, Kpp, Kppi, params, dom, mat, f, BC )
        end
        # ACHTUNG: increment solutions from copied array
        f.Vx[:,2:end-1] .= Vx[:,2:end-1] .+ dx[f.NumVx]
        f.Vy[2:end-1,:] .= Vy[2:end-1,:] .+ dx[f.NumVy]
        f.Pc            .= Pc            .+ dx[f.NumP.+maximum(f.NumVy)]
    end
    return
end
export StokesSolver!

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

@views function DecoupledSolver!( f::Fields2D, params::ModelParameters, fu::Vector{Float64}, fp::Vector{Float64}, Kuu::SparseMatrixCSC{Float64, Int64}, Kup::SparseMatrixCSC{Float64, Int64}, Kpu::SparseMatrixCSC{Float64, Int64} )
    # Decoupled solve
    ndofu = size(Kup,1)
    ndofp = size(Kup,2)
    if params.comp==1
        @tturbo coef      = 1.0./(f.Kc[:].*params.dt).*ones(length(f.etac))#.*etac[:]
        @tturbo coef_inv  = f.Kc[:].*params.dt.*ones(length(f.etac))#.*etac[:]
    else
        @tturbo coef      =    0.0*ones(length(f.etac))#.*etac[:]
        @tturbo coef_inv  = params.gamma*ones(length(f.etac))#.*etac[:]
    end
    @tturbo Kpp   = spdiagm(coef)
    @tturbo Kppi  = spdiagm(coef_inv)
    @tturbo Kuusc = Kuu - Kup*(Kppi*Kpu) # OK
    @tturbo PC    =  0.5*(Kuusc + Kuusc')
    t = @elapsed Kf    = cholesky(Hermitian(PC), check = false)
    @printf("Cholesky took = %02.2e s\n", t)
    u     = zeros(Float64,ndofu)
    ru    = zeros(Float64,ndofu)
    fusc  = zeros(Float64,ndofu)
    p     = zeros(Float64,ndofp)
    rp    = zeros(Float64,ndofp)
    # Iterations
    for rit=1:10
        @tturbo ru   .= fu .- Kuu*u .- Kup*p;
        @tturbo rp   .= fp .- Kpu*u .- Kpp*p;
        @printf("  --> Powell-Hestenes Iteration %02d\n  Momentum res.   = %2.2e\n  Continuity res. = %2.2e\n", rit, norm(ru)/sqrt(length(ru)), norm(rp)/sqrt(length(rp)))
        if norm(ru)/(length(ru)) < 1e-10 && norm(rp)/(length(ru)) < 1e-10
            break
        end
        @tturbo fusc .=  fu .- Kup*(Kppi*fp .+ p)
        @tturbo u    .= Kf\fusc
        @tturbo p   .+= Kppi*(fp .- Kpu*u .- Kpp*p)
    end
    f.Vx[:,2:end-1] .= f.Vx[:,2:end-1] .+ u[f.NumVx]
    f.Vy[2:end-1,:] .= f.Vy[2:end-1,:] .+ u[f.NumVy]
    f.Pc            .= f.Pc            .+ p[f.NumP]
end
export DecoupledSolver!

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

function KSP_GCR_Stokes!( x::Vector{Float64}, M::SparseMatrixCSC{Float64, Int64}, b::Vector{Float64}, eps::Float64, noisy::Int64, Kuu::SparseMatrixCSC{Float64, Int64}, Kup::SparseMatrixCSC{Float64, Int64}, Kpu::SparseMatrixCSC{Float64, Int64}, Kpp::SparseMatrixCSC{Float64, Int64}, Kppi::SparseMatrixCSC{Float64, Int64} )
    # KSP GCR solver
    norm_r, norm0 = 0.0, 0.0
    N         = length(x)
    restart   = 25
    maxit     = 1000
    ncyc, its = 0, 0
    i1, i2, success=0,0,0
    # Arrays for coupled problem
    f      = zeros(Float64, N)
    v      = zeros(Float64, N)
    s      = zeros(Float64, N)
    val    = zeros(Float64, restart)
    VV     = zeros(Float64, (restart,N))
    SS     = zeros(Float64, (restart,N))
    # Coupled
    # Initial residual
    f      = b - M*x
    norm_r = norm(f)
    norm0  = norm_r;
    #
    ndofu = size(Kup,1)
    ndofp = size(Kup,2)
    Kuusc = Kuu - Kup*(Kppi*Kpu) # OK
    PC    =  0.5*(Kuusc + Kuusc')
    t = @elapsed Kf    = cholesky(Hermitian(PC),check = false)
    @printf("Cholesky took = %02.2e s\n", t)
    # Arrays for decoupled problem
    su    = zeros(Float64, ndofu)
    fusc  = zeros(Float64, ndofu)
    sp    = zeros(Float64, ndofp)
    fu    = zeros(Float64, ndofu)
    fp    = zeros(Float64, ndofp)
    @tturbo fu     .= f[1:ndofu]
    @tturbo fp     .= f[ndofu+1:end]
    if (noisy > 1) @printf("       %1.4d KSP GCR Residual %1.12e %1.12e\n", 0, norm_r, norm_r/norm0); end
    # Solving procedure
     while ( success == 0 && its<maxit )
        for i1=1:restart
            # Apply preconditioner, s = PC^{-1} f
            # s = PC\f
            @tturbo fusc .= fu  - Kup*(Kppi*fp + sp)
            @tturbo su   .= Kf\fusc
            @tturbo sp   .+= Kppi*(fp - Kpu*su)
            @tturbo s[1:ndofu]     .= su
            @tturbo s[ndofu+1:end] .= sp
            # Action of Jacobian on s: v = J*s
            # JacobianAction!(v, M, s; r,kv,T,fc,TW,TE,dx,n)
            @tturbo v .= M*s
            # Approximation of the Jv product
            for i2=1:i1
                @tturbo val[i2] = v' * VV[i2,:]
            end
            # Scaling
            for i2=1:i1
                @tturbo v .-= val[i2] * VV[i2,:]
                @tturbo s .-= val[i2] * SS[i2,:]
            end
            # -----------------
            @tturbo r_dot_v = f'*v
            nrm     = norm(v)
            r_dot_v = r_dot_v / nrm
            # -----------------
            fact    = 1.0/nrm
            @tturbo v     .*= fact
            @tturbo s     .*= fact
            # -----------------
            fact    = r_dot_v;
            @tturbo x     .+= fact*s
            @tturbo f     .-= fact*v
            # -----------------
            norm_r  = norm(f)
            @tturbo fu     .= f[1:ndofu]
            @tturbo fp     .= f[ndofu+1:end]
            @printf("  --> Powell-Hestenes Iteration %02d\n  Momentum res.   = %2.2e\n  Continuity res. = %2.2e\n", its, norm(fu)/sqrt(length(fu)), norm(fp)/sqrt(length(fp)))
            if norm(fu)/(length(fu)) < 1e-10 && norm(fp)/(length(fu)) < 1e-10 #(norm_r < eps * norm0 )
                success = 1
                println("converged")
                break
            end
            # Store
            @tturbo VV[i1,:] .= v
            @tturbo SS[i1,:] .= s
            its              += 1
        end
        its  += 1
        ncyc += 1
    end
    if (noisy>1) @printf("[%1.4d] %1.4d KSP GCR Residual %1.12e %1.12e\n", ncyc, its, norm_r, norm_r/norm0); end
    return its
end
export KSP_GCR_Stokes!

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

function KSP_GCR_Stokes2!( x::Vector{Float64}, M::SparseMatrixCSC{Float64, Int64}, X::Vector{Float64}, b::Vector{Float64}, eps::Float64, noisy::Int64, Kuu::SparseMatrixCSC{Float64, Int64}, Kup::SparseMatrixCSC{Float64, Int64}, Kpu::SparseMatrixCSC{Float64, Int64}, Kpp::SparseMatrixCSC{Float64, Int64}, Kppi::SparseMatrixCSC{Float64, Int64}, params::ModelParameters, dom::ModelDomain, materials::MaterialParameters, f::Fields2D, BC::BoundaryConditions )
    # KSP GCR solver
    norm_r, norm0 = 0.0, 0.0
    N         = length(x)
    restart   = 30
    maxit     = 10#1000
    ncyc, its = 0, 0
    i1, i2, success=0,0,0
    epsi      = 1e-5
    ncx, ncy = dom.ncx, dom.ncy
    # Allocate tables
    Tau       = Tensor2D( zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+1, ncy+1), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0) )
    Eps       = Tensor2D( zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+1, ncy+1), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0) )
    F1        = zero(x)
    F2        = zero(x)
    # Arrays for coupled problem
    r      = zeros(Float64, N)
    v      = zeros(Float64, N)
    s      = zeros(Float64, N)
    val    = zeros(Float64, restart)
    VV     = zeros(Float64, (restart,N))
    SS     = zeros(Float64, (restart,N))
    # Coupled
    # Initial residual
    r      .= b# - Mx
    norm_r = norm(r)
    norm0  = norm_r;
    #
    ndofu = size(Kup,1)
    ndofp = size(Kup,2)
    Kuusc = Kuu - Kup*(Kppi*Kpu) # OK
    PC    =  0.5*(Kuusc + Kuusc')
    t = @elapsed Kf    = cholesky(Hermitian(PC),check = false)
    @printf("Cholesky took = %02.2e s\n", t)
    # Arrays for decoupled problem
    su    = zeros(Float64, ndofu)
    fusc  = zeros(Float64, ndofu)
    sp    = zeros(Float64, ndofp)
    fu    = zeros(Float64, ndofu)
    fp    = zeros(Float64, ndofp)
    @tturbo fu     .= r[1:ndofu]
    @tturbo fp     .= r[ndofu+1:end]
    if (noisy > 1) @printf("       %1.4d KSP GCR Residual %1.12e %1.12e\n", 0, norm_r, norm_r/norm0); end

    # res!( F1, X, f, Eps, Tau, dom, BC, params, materials, 1 )
    # println("|Fx| = ", norm(f.Fx)/length(f.Fx))
    # println("|Fy| = ", norm(f.Fy)/length(f.Fy))
    # println("|Fp| = ", norm(f.Fp)/length(f.Fp))

    # Solving procedure
     while ( success == 0 && its<maxit )
        for i1=1:restart
            # Apply preconditioner, s = PC^{-1} f
            # s = PC\f
            @tturbo fusc .= fu  - Kup*(Kppi*fp + sp)
            @tturbo su   .= Kf\fusc
            @tturbo sp   .+= Kppi*(fp - Kpu*su)
            @tturbo s[1:ndofu]     .= su
            @tturbo s[ndofu+1:end] .= sp

            # Action of Jacobian on s: v = J*s
            # @tturbo v .= M*s

            # res!( F1, s, f, Eps, Tau, dom, BC, params, materials, 0 )
            # v .= F1

            res!( F1, X .+ epsi.*s, f, Eps, Tau, dom, BC, params, materials, 1 )
            res!( F2, X .- epsi.*s, f, Eps, Tau, dom, BC, params, materials, 1 )
            v .= (F1 .- F2) ./ ( 2epsi )

            # Approximation of the Jv product
            for i2=1:i1
                @tturbo val[i2] = v' * VV[i2,:]
            end
            # Scaling
            for i2=1:i1
                @tturbo v .-= val[i2] .* VV[i2,:]
                @tturbo s .-= val[i2] .* SS[i2,:]
            end
            # -----------------
            @tturbo r_dot_v = r'*v
            nrm      = sqrt(v'*v)
            r_dot_v /= nrm
            # -----------------
            @tturbo v     ./= nrm
            @tturbo s     ./= nrm
            # -----------------
            @tturbo x     .+= r_dot_v*s
            @tturbo r     .-= r_dot_v*v
            # -----------------
            norm_r  = sqrt(r'*r)
            @tturbo fu     .= r[1:ndofu]
            @tturbo fp     .= r[ndofu+1:end]
            # @printf("  --> Powell-Hestenes Iteration %02d\n  Momentum res.   = %2.2e\n  Continuity res. = %2.2e\n", its, norm(fu)/sqrt(length(fu)), norm(fp)/sqrt(length(fp)))
            if norm(fu)/(length(fu)) < 1e-10 && norm(fp)/(length(fu)) < 1e-10 #(norm_r < eps * norm0 )
                success = 1
                # println("converged")
                break
            end
            # Store
            @tturbo VV[i1,:] .= v[:]
            @tturbo SS[i1,:] .= s[:]
            its              += 1
        end
        its  += 1
        ncyc += 1
    end
    if (noisy>1) @printf("[%1.4d] %1.4d KSP GCR Residual %1.12e %1.12e\n", ncyc, its, norm_r, norm_r/norm0); end
    return its
end
export KSP_GCR_Stokes2!

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

function Rheology_local!( f::Fields2D, Eps::Tensor2D, Tau::Tensor2D, params::ModelParameters, materials::MaterialParameters, BC::BoundaryConditions, ncx::Int64, ncy::Int64 )
    # Compute effective viscosity for each phase
    f.etac .= 0.0
    # f.Kc   .= 0.0
    for m=1:materials.nphase
        @tturbo eta    = 2.0* materials.eta0[m] * Eps.II.^(1.0/materials.n[m] - 1.0)
        @tturbo f.etac .+= f.phase_perc[m,:,:] .* eta
        # @tturbo f.Kc   .+= f.phase_perc[m,:,:] .* materials.K[m]
    end
    CentroidsToVertices!( f.etav, f.etac, ncx, ncy, BC )
    # println("min Eii: ", minimum(Eps.II), " --- max Eii: ", maximum(Eps.II))
    # println("min eta: ", minimum(f.etac), " --- max eta: ", maximum(f.etac))
end
export Rheology_local!

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

function res!( r::Vector{Float64}, x::Vector{Float64}, f::Fields2D, Eps::Tensor2D, Tau::Tensor2D, dom::ModelDomain, BC::BoundaryConditions, params::ModelParameters, materials::MaterialParameters, nl::Int64 )
    f.Vx[:,2:end-1] .= x[f.NumVx]
    f.Vy[2:end-1,:] .= x[f.NumVy]
    f.Pc            .= x[f.NumP .+ maximum(f.NumVy)]
    SetBCs( f.Vx, f.Vy, BC )
    StrainRate!(Eps, f.Vx, f.Vy, dom)
    if nl==1
         Rheology_local!( f, Eps, Tau, params, materials, BC, dom.ncx, dom.ncy  )
    end
    Stress!( Tau, Eps, f )
    ResidualsComp!( f, Tau, Eps, BC, dom, params )
    nVx = length(f.Fx)
    nVy = length(f.Fy)
    r[1:nVx]         .= -f.Fx[:]
    r[nVx+1:nVx+nVy] .= -f.Fy[:]
    r[nVx+nVy+1:end] .= -f.Fp[:]
end
export res!


text:
using DelimitedFiles
using Dates

function get_bonds_v(Nx::Int, Ny::Int)
    bonds = zeros(Int, Nx*Ny*3, 4)
    for i in 0:Nx-1
        for j in 0:Ny-1
            ij = i + j * Nx
            imj = (i - 1 + Nx) % Nx + j * Nx
            ijm = i + ((j - 1 + Ny) % Ny) * Nx
            bonds[3*ij+1, :] = [6*ij+1, 6*ij+2, 6*imj+5, 6*ijm+4]
            bonds[3*ij+2, :] = [6*ij+2, 6*ij+3, 6*ij+5, 6*ijm+6]
            bonds[3*ij+3, :] = [6*ij+1, 6*ij+3, 6*ij+4, 6*imj+6]
        end
    end
    return bonds
end

function Z_gen(Nx::Int, Ny::Int)
    if Nx < 2 || Ny < 2
        throw("Nx, Ny must be at least 2")
    end
    bonds = zeros(Int, 12)
    bonds[1] = 1
    bonds[2] = 2
    bonds[3] = (Ny - 1) * Nx * 6 + 6
    bonds[4] = (1 + (Ny - 1) * Nx) * 6 + 4
    bonds[5] = 8
    bonds[6] = 9
    bonds[7] = 10
    bonds[8] = Nx * 6 + 5
    bonds[9] = Nx * 6 + 3
    bonds[10] = Nx * 6 + 1
    bonds[11] = (2 * Nx - 1) * 6 + 5
    bonds[12] = (Nx - 1) * 6 + 6
    return bonds
end

function X_gen(Nx::Int, Ny::Int)
    vert = zeros(Int, 6)
    vert[1] = 2
    vert[2] = 4
    vert[3] = 6
    vert[4] = Nx * 3 + 2
    vert[5] = Nx * 3 + 1
    vert[6] = 3
    return vert
end

function get_bonds_b(Nx::Int, Ny::Int)
    bonds = Vector{Int}[]
    for j in 0:Ny-1
        for i in 0:Nx-1
            ij = 6 * (i + j * Nx)
            imj = 6 * ((i - 1 + Nx) % Nx + j * Nx)
            ipj = 6 * ((i + 1) % Nx + j * Nx)
            ijm = 6 * (i + ((j - 1 + Ny) % Ny) * Nx)
            ijp = 6 * (i + ((j + 1) % Ny) * Nx)
            imjp = 6 * ((i - 1 + Nx) % Nx + ((j + 1) % Ny) * Nx)
            ipjm = 6 * ((i + 1) % Nx + ((j - 1 + Ny) % Ny) * Nx)

            t = [ij+2, ij+3, ij+4, imj+6, imj+5, ijm+4]
            push!(bonds, ∩(setdiff(1:6*Nx*Ny, t), ij+2:6*Nx*Ny))
            t = [ij+1, ij+3, ij+5, imj+5, ijm+4, ijm+6]
            push!(bonds, ∩(setdiff(1:6*Nx*Ny, t), ij+3:6*Nx*Ny))
            t = [ij+1, ij+2, ij+4, ij+5, imj+6, ijm+6]
            push!(bonds, ∩(setdiff(1:6*Nx*Ny, t), ij+4:6*Nx*Ny))
            t = [ij+1, ij+3, imj+6, ijp+1, ijp+2, imjp+5]
            push!(bonds, ∩(setdiff(1:6*Nx*Ny, t), ij+5:6*Nx*Ny))
            t = [ij+2, ij+3, ipj+1, ipj+2, ijm+6, ipjm+4]
            push!(bonds, ∩(setdiff(1:6*Nx*Ny, t), ij+6:6*Nx*Ny))
            t = [ijp+2, ijp+3, ijp+5, ipj+1, ipj+3, ipj+4]
            push!(bonds, ∩(setdiff(1:6*Nx*Ny, t), ij+7:6*Nx*Ny))
        end
    end
    return bonds
end

function basis(bonds::Vector{Int}, bonds_dic::Vector{Vector{Int}},
        ret::Vector{Int})
    #global basis_all
    if bonds == []
        #println(ret)
        #push!(basis_all,ret)
        global basis_count += 1
        if basis_count % 1000000 == 0
            println(basis_count)
            flush(stdout)
        end
        i = length(ret) + 1
        conf_count[i] += 1
        conf = zeros(Int, 6*Nx*Ny)
        conf[ret] .= 1

        zc = 0
        for l in 1:12
            zc += conf[Z[l]]
        end
        conf_zc[i] += 1 - zc % 2

        zo = 0
        for l in 1:6
            zo += conf[Z[l]]
        end
        conf_zo[i] += 1 - zo % 2

        xc = 1
        for v in 1:6
            if (conf[bonds_v[X[v],1]] + conf[bonds_v[X[v],2]]
                + conf[bonds_v[X[v],3]] + conf[bonds_v[X[v],4]]) == 0
                xc = 0
                break
            end
        end
        conf_xc[i] += xc

        xo = 1
        if (conf[bonds_v[2,1]] + conf[bonds_v[2,2]]
            + conf[bonds_v[2,3]] + conf[bonds_v[2,4]] == 1
            && conf[bonds_v[3,1]] + conf[bonds_v[3,2]]
            + conf[bonds_v[3,3]] + conf[bonds_v[3,4]] == 1
            && conf[7] == 0 && conf[8] == 0 && conf[Nx*6+1] == 0
            && conf[Nx*6+2] == 0)
            xo6 = (conf[Z[1]] + conf[Z[2]] + conf[Z[3]] + conf[Z[4]]
                   + conf[Z[11]] + conf[Z[12]])
            if xo6 % 2 == 1
                conf_xo0[i] += 1
            else
                xo3 = conf[3] + conf[4] + conf[5]
                if (xo6 ÷ 2 + xo3) % 2 == 0
                    conf_xom[i] += 1
                else
                    conf_xop[i] += 1
                end
            end
        end
    else
        basis(Int[],bonds_dic,ret)
        for i in bonds
            #println(i)
            #println(vcat(ret,[i]))
            basis(∩(bonds_dic[i], bonds),bonds_dic,vcat(ret,[i]))
        end
    end
end

function main1(Nx::Int, Ny::Int)
    global bonds_v = get_bonds_v(Nx, Ny)
    global Z = Z_gen(Nx, Ny)
    global X = X_gen(Nx, Ny)
    conf_max = Nx * Ny * 6 ÷ 4 + 1
    #global basis_all = Vector{Int}[]
    global basis_count = 0
    global conf_count = zeros(Int, conf_max)
    global conf_zc = zeros(Int, conf_max)
    global conf_zo = zeros(Int, conf_max)
    global conf_xc = zeros(Int, conf_max)
    global conf_xo0 = zeros(Int, conf_max)
    global conf_xop = zeros(Int, conf_max)
    global conf_xom = zeros(Int, conf_max)
    basis(collect(1:6*Nx*Ny), get_bonds_b(Nx, Ny), Int[])
    #println(basis_all[1:20])
    #println(length(basis_all))
    println(conf_count)
    println(conf_zc)
    println(conf_zo)
    println(conf_xc)
    println(conf_xo0)
    println(conf_xop)
    println(conf_xom)
end

Nx = parse(Int,ARGS[1])
Ny = parse(Int,ARGS[2])

println(Dates.format(now(), "u d HH:MM:SS"))
main1(Nx, Ny)
println(Dates.format(now(), "u d HH:MM:SS"))


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#test-warlpiri-addition.jl

###############################################################################
##COMPLETE ADDITION TABLE

#column key
#none_ pfew_ pone_ psome ptwo_ pmany pall_ junk_ nnone n_few n_one nsome n_two nmany n_all
wadd = Utype[
#none_ + line
 none_ pfew_ pone_ psome ptwo_ pmany pall_ junk_ none_ n_few n_one nsome n_two nmany n_all;
#pfew_ + line
 pfew_ pf_ps psome ps_pm pmany pmany pall_ junk_ pfew_ nf_pf n_few ns_nf nsome nm_ns n_all;
#pone_ + line
 pone_ psome ptwo_ pmany pmany pmany pall_ junk_ pone_ pfew_ none_ n_few n_one nm_ns n_all;
#psome + line
 psome ps_pm pmany pmany pmany pmany pall_ junk_ psome pf_ps pfew_ nf_pf n_few nm_nf n_all;
#ptwo_ + line
 ptwo_ pmany pmany pmany pmany pmany pall_ junk_ ptwo_ psome pone_ pfew_ none_ nm_nf n_all;
#pmany + line
 pmany pmany pmany pmany pmany pmany pall_ junk_ pmany ps_pm ps_pm pf_pm pf_pm nm_pm n_all;
#pall_ + line
 pall_ pall_ pall_ pall_ pall_ pall_ pall_ junk_ pall_ pall_ pall_ pall_ pall_ pall_ junk_;
#junk_ + line
 junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_ junk_;
#nnone + line
 none_ pfew_ pone_ psome ptwo_ pmany pall_ junk_ none_ n_few n_one nsome n_two nmany n_all;
#n_few + line
 n_few nf_pf pfew_ pf_ps psome ps_pm pall_ junk_ n_few ns_nf nsome nm_ns nmany nmany n_all;
#n_one + line
 n_one n_few none_ pfew_ pone_ ps_pm pall_ junk_ n_one nsome n_two nmany nmany nmany n_all;
#nsome + line
 nsome ns_nf n_few nf_pf pfew_ pf_pm pall_ junk_ nsome nm_ns nmany nmany nmany nmany n_all;
#n_two + line
 n_two nsome n_one n_few none_ pf_pm pall_ junk_ n_two nmany nmany nmany nmany nmany n_all;
#nmany + line
 nmany nm_ns nm_ns nm_nf nm_nf nm_pm pall_ junk_ nmany nmany nmany nmany nmany nmany n_all;
#n_all + line
 n_all n_all n_all n_all n_all n_all junk_ junk_ n_all n_all n_all n_all n_all n_all n_all
]

testop(+, wadd)

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
"""
  getBitarrays(corpus,inds;ttype=Int64,atype=Int64,dtype=Int64)

Get the bit array representation of the corpus data.
Return initial topic indices `T`, lexicon index `A`, and the document indices `D` of the tokens.
Shape is `Array{type,2}(data,1,length(data))`.
Optionally specify numeric types for the output.

# Examples
```
corp = map(x->Word.(rand(1:100,x),Int8.(rand(Categorical([0.5,0.3,0.2]),x))),[10,50,100])
T,A = getBitArrays(corp,1:10,Int16,Int16)
```
"""
function getBitarrays(corpus,inds;ttype=Int64,atype=Int64,dtype=Int64)
  tokens = [map(w->w.t,d) for d in corpus[inds]]
  topics = [map(w->w.a,d) for d in corpus[inds]]
  docinds = [map(i->d,corpus[d]) for d in inds]
  T = Array{ttype}(reduce(vcat,tokens))
  A = Array{atype}(reduce(vcat,topics))
  D = Array{dtype}(reduce(vcat,docinds))
  T = reshape(T,(1,length(T)))
  A = reshape(A,(1,length(A)))
  D = reshape(D,(1,length(D)))
  (T=T,A=A,D=D)
end

"""
    getZ(topics,k;inds=nothing)

Generate an indicator-representation of the topic-assignment matrix.
`Z_{i,j} = 1` when token `j` is generated from latent topic `i`

Optionally specify a term index filter `inds` to sample a distributed epoch.

# Examples
```
getZ(StatsBase.sample(1:3,10),3))
```
"""
function getZ(topics,k;inds=nothing)
  if !isnothing(inds)
    topics = topics[inds]
  end
  z = zeros(Bool,k,length(topics))
  for t in 1:length(topics)
    z[topics[t],t] += 1
  end
  z
end

"""
    getJ(corpus;inds=nothing)

Generate the indicator document/topic loadings `J` s.t. `Z*J=nj`.
`J_{i,j} = 1` when token `i` is in document `j`.

Optionally specify a term index filter `inds` to sample a distributed epoch.

# Examples
```
corp = map(x->Word.(rand(1:100,x),Int8.(rand(Categorical([0.5,0.3,0.2]),x))),[10,50,100])
c = getCorp(c)
A,T = getBitarrays(c)
J = getJ(T)
```
"""
function getJ(corpus;inds=nothing)
  tokens = reduce(vcat,[map(w->w.t,d) for d in corpus])
  if !isnothing(inds)
    tokens = filter(!isnothing,indexin(tokens,inds))
  end
  j = zeros(Bool,length(tokens),length(corpus))
  ct = 0
  for d in 1:length(corpus)
    for t in 1:length(corpus[d])
      ct+=1
      j[ct,d]+=1
    end
  end
  j
end

"""
    getJ(corpus,drange,ndocs;inds=nothing)

See `getJ(corpus;inds=nothing)`.  Inserts padding columns for vertical concatenation.

Optionally specify a term index filter `inds` to sample a distributed epoch.

# Examples
```
corp = map(x->Word.(rand(1:100,x),Int8.(rand(Categorical([0.5,0.3,0.2]),x))),[10,50,100])
c = getCorp(c)
Z = getZ(c,3)
J = getJ(c[1:2],1:2,3)
```
"""
function getJ(corpus,drange,ndocs;inds=nothing)
  tokens = reduce(vcat,[map(w->w.t,d) for d in corpus])
  if !isnothing(inds)
    tokens = filter(!isnothing,indexin(tokens,inds))
  end
  j = zeros(Bool,length(tokens),ndocs)
  ct = 0
  for d in 1:length(corpus)
    for t in 1:length(corpus[d])
      ct+=1
      j[ct,drange[d]]+=1
    end
  end
  j
end

"""
    getW(corpus,l;inds=nothing)

Generate the indicator term/topic loadings `W` s.t. `Z*W=nw`.
`W_{i,j} = 1` when token `i` corresponds to term `j` in the lexicon with length `l`.

Optionally specify a term index filter `inds` to sample a distributed epoch.

# Examples
```
corp = map(x->Word.(rand(1:100,x),Int8.(rand(Categorical([0.5,0.3,0.2]),x))),[10,50,100])
c = getCorp(c)
Z = getZ(c,3)
W = getW(c)
```
"""
function getW(corpus,l;inds=nothing)
  tokens = reduce(vcat,[map(w->w.t,d) for d in corpus])
  if !isnothing(inds)
    tokens = filter(!isnothing,indexin(tokens,inds))
  end
  w = zeros(Bool,length(tokens),l)
  for t in 1:l
    w[:,t] .= .!isnothing.(indexin(tokens,[t]))
  end
  w
end

"""
  nj2token(nj,J)

Expand `nj` (the topic x document counts) over tokens.
`NJK_{i,j} = n` where `n` is the number of times topic `j` is assigned to some token in the parent document of token `i`.

# Examples
```
corp = map(x->Word.(rand(1:100,x),Int8.(rand(Categorical([0.5,0.3,0.2]),x))),[10,50,100])
c = getCorp(c)
Z = getZ(c,3)
J = getJ(c)
nj2token(Z*J,J)
```
"""
function nj2token(nj,J)
  d = map(i->i[2],findmax(J,dims=2)[2])
  NJK = nj[:,vec(d)]
end

"""
  nw2token(nw,W)

Expand `nw` (the topic x term counts) over tokens.
`NWK_{i,j} = n` where `n` is the number of times topic `j` is assigned to some token in the parent document of token `i` corresponding to the same term in the lexicon.

# Examples
```
corp = map(x->Word.(rand(1:100,x),Int8.(rand(Categorical([0.5,0.3,0.2]),x))),[10,50,100])
c = getCorp(c)
Z = getZ(c,3)
W = getW(c)
nj2token(Z*W,W)
```
"""
function nw2token(nw,W)
  t = map(i->i[2],findmax(W,dims=2)[2])
  NWK = nw[:,vec(t)]
end


text:
# NOTE: This also updates the properties of the elements of the vector
function update_vector_property!(updated_object::IEntity,
                                 updated_property::Symbol,
                                 dbconn::Union{Missing, LibPQ.Connection};
                                 editor::Union{Missing, IEntity} = missing)

    # Missing is ignored, one needs to pass an empty vector to remove all the associations
    if ismissing(getproperty(updated_object, updated_property))
        return updated_object
    end

    orm_module = PostgresORM.get_orm(updated_object)

    counterpart_datatype = util_get_onetomany_counterparts(orm_module)[updated_property][:data_type]
    counterpart_property = util_get_onetomany_counterparts(orm_module)[updated_property][:property]
    action_on_remove = util_get_onetomany_counterparts(orm_module)[updated_property][:action_on_remove]

    # Update the elements in the vector to make sure that they all point to the
    #    updated object.
    #  This is for convenience so that it is sufficient to add an element to the
    #    vector to create the reference to the parent when we persist to the db.
    children = getproperty(updated_object, updated_property)

    for e in children
        setproperty!(e,
                     counterpart_property,
                     updated_object)
    end

    # Retrieve the current state of the database
    filter_object = counterpart_datatype()
    setproperty!(filter_object,counterpart_property,updated_object)
    previous_vector = retrieve_entity(filter_object,
                                      false, # do not retrieve complex properties
                                      dbconn
                                      )

    children = util_compare_and_sync_entities(previous_vector,
                                              children,
                                              counterpart_property,
                                              action_on_remove,
                                              dbconn;
                                              editor = editor)

    # @show children

    setproperty!(updated_object,
                 updated_property,
                 children)

end

# This function is useful when trying to update a property of type vector
# eg.
# country_deal_asso = CountryDealAsso(;deal = mydeal)
# old_countries_of_deals = retrieve_countrydeal_asso(country_deal_asso) # retrieve the previous state from the database
# util_compare_and_sync_entities(old_countries_of_deals, # previous_entities
#                                mydeal.countrydeal_assos, # new_entities
#                                dbconn;
#                                editor = superadmin)
function util_compare_and_sync_entities(previous_entities::Vector{<:IEntity},
                                        new_entities::Vector{<:IEntity},
                                        counterpart_property::Symbol,
                                        action_on_remove::CRUDType.CRUD, # either update of delete
                                        dbconn::Union{Missing,LibPQ.Connection};
                                        editor::Union{Missing, IEntity} = missing)

    if length(previous_entities) == 0 && length(new_entities) == 0
       return previous_entities
    end

    # Put the ids in a vector of dicts
    previous_ids = Vector{Dict{String,Any}}()
    for e in previous_entities
        push!(previous_ids, util_get_ids_cols_names_and_values(e,dbconn))
    end
    new_ids = Vector{Dict{String,<:Any}}()
    for e in new_entities
        push!(new_ids, util_get_ids_cols_names_and_values(e,dbconn))
    end

    # @info "previous_ids[$previous_ids]"
    # @info "new_ids[$new_ids]"

    # Retrieve the data_type
    if length(previous_entities) > 0
        data_type = typeof(previous_entities[1])
    else
        data_type = typeof(new_entities[1])
    end

    # ############################################################################## #
    # PART1. Identify the indexes of the entities for the different types of actions
    #   (i.e. remove the entity from the association, create and add an entity to
    #         the association, add an existing entity to the association)
    # ############################################################################## #

    # Loop over the previous_ids to see which ones need to be deleted/nulled
    indexes_of_entities_to_remove = [] # Refers to indexes in 'previous_entities'
    counter = 0
    ids_to_remove = Vector{Dict{String,<:Any}}()
    for previous_id in previous_ids
        counter+=1
        previous_id_still_here = false # Initialize
        for new_id in new_ids

            # If the new_id has a missing value we skip the comparison because
            #   1. It cannot be the previous_id (because it has a complete id)
            #   2. The comparison may return missing if none of the component of
            #         the dict is not missing
            if any(x -> ismissing(x), collect(values(new_id)))
                continue
            end

            if previous_id == new_id
                previous_id_still_here = true
                break
            end
        end
        # If we went through all the new ids without finding the id then it means
        #   that the entity is no longer part of the associations
        if !previous_id_still_here
            push!(ids_to_remove, previous_id)
            push!(indexes_of_entities_to_remove, counter)
        end
    end

    # Loop over the new_ids to see which ones need to be created/updated
    indexes_of_entities_to_create = [] # Refers to indexes in 'news_entities'
    indexes_of_entities_to_add = [] # Refers to indexes in 'news_entities'
    counter = 0
    for new_id in new_ids

        counter+=1

        # If one of the PK is null, create the entity
        if any(x -> ismissing(x), collect(values(new_id)))
            push!(indexes_of_entities_to_create,
                  counter)
            continue
        end

        # The id has already been added to the list of entities to be removed, skip
        if new_id in ids_to_remove
            continue
        end

        # If we cannot find an entity in the DB with this id, create it
        filter_object = data_type()
        for p in util_get_ids_props_names(filter_object)
            setproperty!(filter_object, p,
                         getproperty(new_entities[counter], p))
        end
        if length(retrieve_entity(filter_object,false, dbconn)) == 0 push!(indexes_of_entities_to_create,
                  counter)
            continue
        end

        # In other cases, update the target entity
        push!(indexes_of_entities_to_add,
              counter)
    end

    # ################################################################### #
    # PART2. Perform on the actions on the entities at the given indexes
    # ################################################################### #

    # Remove some associations
    for i in indexes_of_entities_to_remove
        filter_object = data_type()
        for p in util_get_ids_props_names(filter_object)
            setproperty!(filter_object,p,
                         getproperty(previous_entities[i],p))
        end
        # @info "Remove association to entity[$(util_prettify_id_values_as_str(filter_object))]"

        if action_on_remove == CRUDType.delete
            delete_entity(filter_object, dbconn; editor = editor)
        elseif action_on_remove == CRUDType.update
            # Retrieve the latest version of the database, we ignore the
            #   other possible changes
            removed_object = retrieve_entity(filter_object,
                                             false, # retrieve complex properties
                                             dbconn)[1]
            setproperty!(removed_object,counterpart_property,missing)
            update_entity!(removed_object,dbconn; editor = editor)

        else
            throw(DomainError("You can only choose 'delete' or 'update' when updating a vector property"))
        end

    end

    result = [] # Store the created/updated entities

    # Create new entities
    for i in indexes_of_entities_to_create
        filter_object = data_type()
        for p in util_get_ids_props_names(filter_object)
            setproperty!(filter_object,p,
                         getproperty(new_entities[i],p))
        end
        # @info "Create and add to association entity [$(util_prettify_id_values_as_str(filter_object))]"
        push!(result, create_entity!(new_entities[i],dbconn;creator = editor))
    end

    # Update existing entities
    for i in indexes_of_entities_to_add
        filter_object = data_type()
        for p in util_get_ids_props_names(filter_object)
            setproperty!(filter_object,p,
                         getproperty(new_entities[i],p))
        end
        # @info "Add to association existing entity [$(util_prettify_id_values_as_str(filter_object))]"
        push!(result, update_entity!(new_entities[i],dbconn;editor = editor))
    end

    return result

    # OLD Version with single PK (can be removed in a few months)
    # # Retrieve the id_property from one of the 2 vectors
    # if length(previous_entities) > 0
    #     id_property = PostgresORM.get_orm(previous_entities[1]).id_property
    #     data_type = typeof(previous_entities[1])
    # else
    #     id_property = PostgresORM.get_orm(new_entities[1]).id_property
    #     data_type = typeof(new_entities[1])
    # end
    #
    # # Get the ids (excluding the missing) of the previous and the new entities
    # # NOTE: the previous entities should all have an id
    # existings_ids_in_previous_entities =
    #     filter(x -> !ismissing(x),
    #            getproperty.(previous_entities,id_property))
    # existings_ids_in_new_entities =
    #    filter(x -> !ismissing(x),
    #           getproperty.(new_entities,id_property))
    #
    #
    # # The objects to delete are the ones whose ids do not exist in the
    # #    new list of entities
    # ids_to_remove = filter(x -> ! (x in existings_ids_in_new_entities),
    #                        existings_ids_in_previous_entities)
    #
    # # Delete the objects that must be deleted
    # for id in ids_to_remove
    #     filter_object = data_type()
    #     setproperty!(filter_object,id_property,id)
    #     if action_on_remove == CRUDType.delete
    #         delete_entity(filter_object, dbconn; editor = editor)
    #     elseif action_on_remove == CRUDType.update
    #         # Retrieve the latest version of the database, we ignore the
    #         #   other possible changes
    #         removed_object = retrieve_entity(filter_object,
    #                                          false, # retrieve complex properties
    #                                          dbconn)[1]
    #         setproperty!(removed_object,counterpart_property,missing)
    #         update_entity!(removed_object,dbconn; editor = editor)
    #
    #     else
    #         throw(DomainError("You can only choose 'delete' or 'update' when updating a vector property"))
    #     end
    # end
    #
    # # Create or update the other ones
    # result = []
    # for e in new_entities
    #     # If the entity has no id set it means it needs to be created
    #     if ismissing(getproperty(e,
    #                              id_property))
    #         push!(result, create_entity!(e,dbconn;creator = editor))
    #     else
    #         push!(result, update_entity!(e,dbconn;editor = editor))
    #     end
    # end
    #
    # return(result)

end

# query_string = "UPDATE appuser SET lastname = \S1, login = \$2
#                 WHERE id = \$3"
function update_entity!(updated_object::IEntity,
                       table_name::String,
                       columns_selection_and_mapping::Dict,
                       dbconn::LibPQ.Connection;
                       editor::Union{IEntity,Missing} = missing
                       )

    data_type = typeof(updated_object)
    orm_module = PostgresORM.get_orm(updated_object)

    query_string = "UPDATE " * table_name * " SET "

    props = util_get_entity_props_for_comparison(updated_object, dbconn)

    # Check if something has changed
    # NOTE: This is done for both deciding whether there is a need to do some
    #         changes to the database and for tracking the modifications
    previous_state_filter = data_type()
    idprops = util_get_ids_props_names(previous_state_filter)
    PostgresORMUtil.setpropertiesvalues!(previous_state_filter,
                                           idprops,
                                           PostgresORMUtil.getpropertiesvalues(updated_object,idprops))
    previous_state = retrieve_entity(previous_state_filter, orm_module,
                                     false, # we want the ids for comparison
                                     dbconn)

    if length(previous_state) == 0
        throw(DomainError("Unable to retrieve an object of type[$data_type] "
                        * "with id[$(getproperty(updated_object,
                                     orm_module.id_property))] from the database."
                        *" Remind that only existing existing objects can be updated."))
    end

    previous_state = previous_state[1]

    previous_state_props = util_get_entity_props_for_comparison(previous_state,
                                                                dbconn)

    changes = diff_dict(previous_state_props, props)

    # If no changes were made, we just return the object
    if length(changes) == 0
        return updated_object
    end

    #
    # Enrich the updated object with the creator and creation time if needed
    #
    editor_property = util_get_editor_property(orm_module)
    if !ismissing(editor_property)

        # Only set the editor property if it is empty because the rest of the
        #  application may want to set it (eg. if importing data)
        if ismissing(getproperty(updated_object,editor_property))

            setproperty!(updated_object,
                         editor_property,
                         editor)
        # If the object has the creator property set, use it
        else
            editor = getproperty(updated_object,editor_property)
        end
    end
    update_time_property = util_get_update_time_property(orm_module)
    if !ismissing(update_time_property)
        setproperty!(updated_object,
                     update_time_property,
                     now(Dates.UTC))
    end

    # Compare versions if needed, by retrieving the current state of the object
    #   in the database
    if util_get_track_changes(orm_module)

         #
         # Create the modification entries
         #
         action_id = uuid4()

         # Create a string representation of the ID properties values.
         id_value_as_str = util_prettify_id_values_as_str(updated_object)

         # Get the editor id if any
         editor_id = if ismissing(editor)
           missing
         else
             util_prettify_id_values_as_str(editor)
         end

        for (k,v) in changes

             modif = Modification(entity_type =  string(data_type),
                                  entity_id = id_value_as_str,
                                  attrname = string(k),
                                  oldvalue = ismissing(v.old) ? missing : string(v.old),
                                  newvalue = ismissing(v.new) ? missing : string(v.new),
                                  appuser_id = editor_id,
                                  creation_time = now(Dates.UTC),
                                  action_type = CRUDType.update,
                                  action_id = action_id)
             create_entity!(modif,dbconn)
         end

     end # ENDOF if track_changes

    # Create the dictionnary of properties for insertion in the database
    props = util_get_entity_props_for_db_actions(updated_object,
                                                 dbconn,
                                                 true # include missing for updating columns to null
                                                 )
    properties_names = collect(keys(props))
    properties_values = collect(values(props))

    # The prepared statement has arguments corresponding to both the 'SET' clause
    #   and the 'WHERE' clause. Therefore we must increment the same counter as
    #   we loops over the properties that are updated and the ID properties.
    prepared_statement_args_counter = 0
    prepared_statement_args_values = Any[] # Holds the values of the prepared statement

    ids_cols_names = util_get_ids_cols_names(updated_object)

    # Add the 'SET' clause:
    # Loop over the properties that are updated
    # NOTE: we could remove the ID properties from this list but that's alright
    if length(props) > 0
        for (index,property_name) in enumerate(properties_names)

            # Add the restriction
            # eg. 'id = \$3'

            # NOTE: The mapping may not have the properties because some of the
            #          properties are the columns names of the complex objects
            #          which are added at runtime
            columnname = if haskey(columns_selection_and_mapping, property_name)
                            columns_selection_and_mapping[property_name]
                         else
                             property_name
                         end

             # If property is associated to an ID column, skip
             if string(columnname) in ids_cols_names
                 continue
             end

             prepared_statement_args_counter += 1
             # Add ',' if it's not the first property
             if prepared_statement_args_counter > 1
                 query_string *= ", "
             end

            query_string *= " $columnname = \$$prepared_statement_args_counter "

            # Store the new value of this property in the list of the values to
            #   be passed at execution of the prepared statement
            push!(prepared_statement_args_values,
                  props[Symbol(property_name)])
        end
    else
        # TODO throw_error
    end

    # Add the 'WHERE' clause:
    # Loop over the ID properties
    query_string *= " WHERE "
    where_clause_elts = Vector{String}()
    # An object can have several ID properties
    for id_property in util_get_ids_props_names(updated_object)
        # An ID property can have several columns (eg. when the property is a
        #   complex property that has several IDs)
        id_columnnames = tovector(columns_selection_and_mapping[id_property])
        for id_columnname in id_columnnames
            where_clause_elt =  "$id_columnname = \$$(prepared_statement_args_counter+=1)"
            push!(where_clause_elts, where_clause_elt)
            # If the ID property is a complex prop, then the value is stored
            #   under the name of the column name
            if isa(getproperty(updated_object, id_property),IEntity)
                push!(prepared_statement_args_values,
                      props[Symbol(id_columnname)])
            else
                push!(prepared_statement_args_values,
                      props[id_property])
            end

        end
    end
    query_string *= join(where_clause_elts, " AND ")

    query_string *= " RETURNING *";

    prepared_query = LibPQ.prepare(dbconn,
                           query_string)

    query_result = execute(prepared_query,
                           prepared_statement_args_values
                           ;throw_error=true)

     # Get the result as a dataframe
     data = DataFrame(query_result)

     # Close the result
     close(query_result)

     # Format the result to be mapped to the data model
     result = util_clean_and_rename_query_result(data,columns_selection_and_mapping)

     result = dataframe2vector_of_namedtuples(result)

     result =
        util_convert_namedtuple_to_object.(result,data_type,
                                           true, # retrieve_complex_props
                                           dbconn)

     util_overwrite_props!(result[1], # src entity
                           updated_object, # target entity
                           true # Ignore missing properties like all the
                                #    unmapped properties)
                           )

     return updated_object

end

function update_entity!(updated_object::IEntity,
                       orm_module::Module,
                       dbconn::LibPQ.Connection;
                       editor::Union{IEntity,Missing} = missing)
    update_entity!(updated_object,
                  util_get_table_with_schema_name(orm_module),
                  util_get_columns_selection_and_mapping(orm_module),
                  dbconn;
                  editor = editor)
end

function update_entity!(updated_object::IEntity,
                       dbconn::LibPQ.Connection;
                       editor::Union{IEntity,Missing} = missing
                       )
    update_entity!(updated_object,
                  PostgresORM.get_orm(updated_object),
                  dbconn;
                  editor = editor)
end


text:
projected_positions_ECI2ECEF = [
    [-8.537268174201572e6, -1.2988350094779432e7, 2.1851582648954894e7],
    [7.172876999458652e6, -2.1714421682282157e7, 1.3296862881729363e7],
    [2.401937266433255e7, 1.1645127140291367e7, 214656.1633842468],
    [-5.021627715111267e6, 2.2734776424851105e7, 1.2106563114258626e7],
    [8.06727794696111e6, 1.280426674805988e7, 2.181418989894819e7]
    ]
@testset "Sat Position ECI2ECEF" begin
    for i in 1:length(satellite_states)
        out = sat_position_ECI_2_ECEF(satellite_states[i])
        @test out ≈ projected_positions_ECI2ECEF[i]
    end

    sv_bad = deepcopy(sv1_struct)
    sv_bad.decoder_state.subframes_decoded = [0,1,1,1,1]

    out = false
    try
        sat_position_ECEF(sv_bad)
    catch e
        if typeof(e) == PositionVelocityTime.BadData
            out = true
        end
    end
    @test out == true

    @test out == true
end

projected_positions_ECEF = [
    [-8.537268174201585e6, -1.2988350094779423e7, 2.1851582648954894e7],
    [7.172876999458676e6, -2.171442168228215e7, 1.3296862881729363e7],
    [2.4019372664332554e7, 1.1645127140291365e7, 214656.1633842468],
    [-5.021627715111248e6, 2.273477642485111e7, 1.2106563114258626e7],
    [8.067277946961094e6, 1.280426674805989e7, 2.181418989894819e7]
    ]
@testset "Sat Position ECEF" begin
    for i in 1:length(satellite_states)
        out = sat_position_ECEF(satellite_states[i])
        @test out == projected_positions_ECEF[i]
    end

    sv_bad = deepcopy(sv1_struct)
    sv_bad.decoder_state.subframes_decoded = [0,0,0,0,0]
    out = false
    try
        sat_position_ECEF(sv_bad)
    catch e
        if typeof(e) == PositionVelocityTime.BadData
            out = true
        end
    end
    @test out == true
end


text:
using PSCOPF

using Test
using JuMP
using Dates
using DataStructures

@testset verbose=true "test_energy_market_dmo" begin

    #=
    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (limitable) wind_1_1|load
    Pmin=0, Pmax=100    |  S1:55
    Csta=0, Cprop=1     |
      S1: 20            |
                        |
    (pilotable) prod_1_1|
    Pmin=10, Pmax=100   |
    Csta=0, Cprop=10    |
    DMO => 8h           |
                        |
    (pilotable) prod_1_2|
     Pmin=10, Pmax=100  |
     Csta=0, Cprop=15   |
     DMO => 10h30       |
                        |
    =#
    function create_instance(ech, next_ech, ts,
                            prod_1_state, prod_1_state_definitive, prod_1_level,
                            prod_2_state, prod_2_state_definitive, prod_2_level,
                            )
        ECH = [DateTime("2015-01-01T07:00:00"), DateTime("2015-01-01T09:00:00"), DateTime("2015-01-01T10:35:00")]
        @assert (ech in ECH)
        network = PSCOPF.Networks.Network()
        # Buses
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                10., 100.,
                                                0., 10.,
                                                Dates.Second(3*60*60), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_2", PSCOPF.Networks.PILOTABLE,
                                                10., 100.,
                                                0., 15.,
                                                Dates.Second(30*60), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ECH[1], "wind_1_1", ts, "S1", 20.)
        PSCOPF.add_uncertainty!(uncertainties, ECH[1], "bus_1", ts, "S1", 55.)
        PSCOPF.add_uncertainty!(uncertainties, ECH[2], "wind_1_1", ts, "S1", 20.)
        PSCOPF.add_uncertainty!(uncertainties, ECH[2], "bus_1", ts, "S1", 55.)
        PSCOPF.add_uncertainty!(uncertainties, ECH[3], "wind_1_1", ts, "S1", 20.)
        PSCOPF.add_uncertainty!(uncertainties, ECH[3], "bus_1", ts, "S1", 55.)
        # initial generators state : need to pay starting cost at ts
        generators_init_state = SortedDict(
            "prod_1_1" => PSCOPF.OFF,
            "prod_1_2" => PSCOPF.OFF
        )
        mode = PSCOPF.ManagementMode("mode_5mins", Dates.Minute(5))

        prod_1_state_def = prod_1_state_definitive ? prod_1_state : missing
        prod_2_state_def = prod_2_state_definitive ? prod_2_state : missing
        initial_schedule = PSCOPF.Schedule(PSCOPF.Market(), ech - Second(1), SortedDict(
            "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                # SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(missing,
                #                                                                          SortedDict("S1"=>PSCOPF.ON))),
                SortedDict(ts => PSCOPF.UncertainValue{Float64}(missing,
                                                                SortedDict("S1"=>20.)))
                ),
            "prod_1_1" => PSCOPF.GeneratorSchedule("prod_1_1",
                SortedDict(ts => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(prod_1_state_def,
                                                                SortedDict("S1"=>prod_1_state))),
                SortedDict(ts => PSCOPF.UncertainValue{Float64}(missing,
                                                                SortedDict("S1"=>prod_1_level)))
                ),
            "prod_1_2" => PSCOPF.GeneratorSchedule("prod_1_2",
                SortedDict(ts => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(prod_2_state_def,
                                                                SortedDict("S1"=>prod_2_state))),
                SortedDict(ts => PSCOPF.UncertainValue{Float64}(missing,
                                                                SortedDict("S1"=>prod_2_level)))
                )
            )
        )

        context = PSCOPF.PSCOPFContext(network, [ts], mode,
                                        generators_init_state,
                                        uncertainties, nothing)

        context.tso_schedule = initial_schedule

        market = PSCOPF.EnergyMarket()

        firmness = PSCOPF.compute_firmness(market, ech,
                                            next_ech,
                                            [ts], context)

        return context, market, firmness
    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (limitable) wind_1_1|load
    Pmin=0, Pmax=100    |  S1:55
    Csta=0, Cprop=1     |
      S1: 20            |
                        |
    (pilotable) prod_1_1|
    Pmin=10, Pmax=100   |
    Csta=0, Cprop=10    |
    DMO => 8h           |
                        |
    (pilotable) prod_1_2|
     Pmin=10, Pmax=100  |
     Csta=0, Cprop=15   |
     DMO => 10h30       |
                        |

    We can start an pilotable before its DMO.
    ech = 7h
    TS = 11h
    DMO(prod_1_1) = 3h => 8h
    DMO(prod_1_2) = 30mins => 10h30

    We suppose that in the preceding step we decided :
    prod_1_1 : OFF (non definitive)
    prod_1_2 : ON (non definitive)

    but in terms of cost, it is cheaper to use prod_1_1.
    Since we haven't reached the DMO yet, we can start prod_1_1
    => the new decision :
    prod_1_1 : ON (unit started)
    prod_1_2 : OFF (unit shutdown)
    =#
    @testset "energy_market_can_start_unit_when_commitment_firmness_is_FREE" begin
        ech = DateTime("2015-01-01T07:00:00")
        next_ech = DateTime("2015-01-01T08:00:00")
        TS = [DateTime("2015-01-01T11:00:00")]
        context, market, firmness = create_instance(ech, next_ech, TS[1],
                                                    PSCOPF.OFF, false, 0.,
                                                    PSCOPF.ON, false, 35.,
                                                    )

        # firmness
        expected_firmness = PSCOPF.Firmness(
                SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                            "prod_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), ),
                SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                            "prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                            "prod_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), )
                )
        @test firmness == expected_firmness

        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
        @test 20. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
        # we started prod_1_1
        @test PSCOPF.ON == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_1", TS[1], "S1")
        @test 35. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1")
        # we shut down prod_1_2
        @test PSCOPF.OFF == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_2", TS[1], "S1")
        @test PSCOPF.get_prod_value(context.market_schedule, "prod_1_2", TS[1], "S1") < 1e-09
    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (limitable) wind_1_1|load
    Pmin=0, Pmax=100    |  S1:55
    Csta=0, Cprop=1     |
      S1: 20            |
                        |
    (pilotable) prod_1_1|
    Pmin=10, Pmax=100   |
    Csta=0, Cprop=10    |
    DMO => 8h           |
                        |
    (pilotable) prod_1_2|
     Pmin=10, Pmax=100  |
     Csta=0, Cprop=15   |
     DMO => 10h30       |
                        |

    We can no longer start an pilotable after its DMO.
    ech = 9h
    TS = 11h
    DMO(prod_1_1) = 3h => 8h => We are past DMO!
    DMO(prod_1_2) = 30mins => 10h30

    We suppose that in the preceding step we decided :
    prod_1_1 : OFF (definitive)
    prod_1_2 : ON (non definitive)

    but in terms of cost, it is cheaper to use prod_1_1.
    We would like to start prod_1_1 cause it's cheaper. But we are past DMO => we can no longer start it.
    => we are stuck using prod_1_2 :
    prod_1_1 : OFF
    prod_1_2 : ON
    =#
    @testset "energy_market_cannot_start_unit_after_DMO" begin
        ech = DateTime("2015-01-01T09:00:00")
        next_ech = DateTime("2015-01-01T10:00:00")
        TS = [DateTime("2015-01-01T11:00:00")]
        context, market, firmness = create_instance(ech, next_ech, TS[1],
                                                    PSCOPF.OFF, true, 0.,
                                                    PSCOPF.ON, false, 35.,
                                                    )

        # firmness : prod_1_1 is already decided (DMO > ECH)
        expected_firmness = PSCOPF.Firmness(
                SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED),
                            "prod_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), ),
                SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                            "prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                            "prod_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), )
                )
        @test firmness == expected_firmness

        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
        @test 20. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
        # we could not start prod_1_1 due to DMO even if it is cheaper
        @test PSCOPF.OFF == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_1", TS[1], "S1")
        @test PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1") < 1e-09
        # we use prod_1_2
        @test PSCOPF.ON == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_2", TS[1], "S1")
        @test 35. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_2", TS[1], "S1")
    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (limitable) wind_1_1|load
    Pmin=0, Pmax=100    |  S1:55
    Csta=0, Cprop=1     |
      S1: 20            |
                        |
    (pilotable) prod_1_1|
    Pmin=10, Pmax=100   |
    Csta=0, Cprop=10    |
    DMO => 8h           |
                        |
    (pilotable) prod_1_2|
     Pmin=10, Pmax=100  |
     Csta=0, Cprop=15   |
     DMO => 10h30       |
                        |

    We can still shutdown a unit after DMO.
    We can still change a unit's level after DMO (before DP).

    ech = 10h35
    TS = 11h
    DMO(prod_1_1) = 3h => 8h => We are past DMO!
    DMO(prod_1_2) = 30mins => 10h30 => We are past DMO!

    We suppose that in the preceding step we decided :
    prod_1_1 : ON (definitive) at 15MW
    prod_1_2 : ON (definitive) at 20MW

    but in terms of cost, it is cheaper to only use prod_1_1.
    => we can change the levels to :
    prod_1_1 : ON at 35MW (production level changed)
    prod_1_2 : OFF (unit was shutdown)
    =#
    @testset "energy_market_can_shutdown_unit_after_DMO" begin
        ech = DateTime("2015-01-01T10:35:00")
        next_ech = DateTime("2015-01-01T10:45:00")
        TS = [DateTime("2015-01-01T11:00:00")]
        context, market, firmness = create_instance(ech, next_ech, TS[1],
                                                    PSCOPF.ON, true, 15.,
                                                    PSCOPF.ON, true, 20.,
                                                    )
        # firmness : prod_1_1 is already decided (DMO > ECH)
        expected_firmness = PSCOPF.Firmness(
                SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED),
                            "prod_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED), ),
                SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                            "prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                            "prod_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), )
                )
        @test firmness == expected_firmness

        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
        @test 20. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
        # we use prod_1_1, it was already started, we can change production before DP
        @test PSCOPF.ON == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_1", TS[1], "S1")
        @test 35. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1")
        # we can shutdown prod_1_2 after the DMO
        @test PSCOPF.OFF == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_2", TS[1], "S1")
        @test PSCOPF.get_prod_value(context.market_schedule, "prod_1_2", TS[1], "S1") < 1e-09
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using CuArrays
using LinearAlgebra
using IterativeSolvers

#This file contains the incomplete factorization routine

#This function computes the inner procuct of two sparse vectors that are
#stores using the same nzind and nzval, but different offsets iter1:u1 resp.
#iter2:u2
function _innerprod( iter1::Int64,
                     u1::Int64,
                     iter2::Int64,
                     u2::Int64,
                     nzind::Array{Int64,1},
                     nzval::Array{Float64,1} )
  @inbounds begin
    res = zero(Float64)
    while ( iter1 <= u1 ) && ( iter2 <= u2 )
      while (nzind[iter1] == nzind[iter2])  &&( iter1 <= u1 ) && ( iter2 <= u2 )
        res += nzval[iter1] * nzval[iter2]
        iter1 += 1
        iter2 += 1
      end
      if nzind[iter1] < nzind[iter2]
        iter1 += 1
      else
        iter2 += 1
      end
    end
    return res
  end
end

function _innerprod( range1,
                     range2,
                     nzind::Array{Int64,1},
                     nzval::Array{Float64,1} )
  iter1 = range1.start
  u1 = range1.stop

  iter2 = range2.start
  u2 = range2.stop
  @inbounds begin
    res = zero(Float64)
    while ( iter1 <= u1 ) && ( iter2 <= u2 )
      while (nzind[iter1] == nzind[iter2])  &&( iter1 <= u1 ) && ( iter2 <= u2 )
        res += nzval[iter1] * nzval[iter2]
        iter1 += 1
        iter2 += 1
      end
      if nzind[iter1] < nzind[iter2]
        iter1 += 1
      else
        iter2 += 1
      end
    end
    return res
  end
end

##This function computes the in place incomplete Cholesky factorisation using
##an inner product wrapped in a function.
#function _el_icholRWHL!( s::Array{Float64},
#                          ind::Array{Int64, 1},
#                          jnd::Array{Int64, 1},
#                          )
#  @inbounds begin
#    #piv_i goes from 1 to N and designates the number, which is being treated
#    for piv_i = 1 : ( size( ind, 1 ) - 1 )
#      #piv_j ranges over pointers and designates where to find the present column
#      for piv_j = ind[ piv_i ] : ( ind[ piv_i + 1 ] - 1 )
#        #iter designates the pointer, at which the i-iterator starts
#        iter::Int64 = ind[ piv_i ]
#        #jter designates the pointer, at which the j-iterator starts
#        jter::Int64 = ind[ jnd[ piv_j ] ]
#        #The condition makes sure, that only the columns up to th pivot
#        #are being summed.
#
#        s[ piv_j ] -= _innerprod( iter, ind[ piv_i + 1 ] - 2,
#                                  jter, ind[ jnd[ piv_j ] + 1 ] - 2,
#                                  jnd,
#                                  s)
#
#        if jnd[ piv_j ] < piv_i
#          s[ piv_j ] /= ( s[ ind[ jnd[ piv_j ] + 1 ] - 1 ] )
#        elseif jnd[ piv_j ] == piv_i
#          if ( s[ ind[ piv_i + 1 ] - 1 ] ) <= 0.
#            s[ ind[ piv_i + 1 ] - 1 ] = 0
#            #Treating the case where the matrix is near-low rank and our
#            #factorization returns a low-rank matrix.
#            maxj = piv_i
#            for newpiv_i = ( piv_i + 1 ) : ( size( ind, 1 ) - 1 )
#              newpiv_j::Int64 = ind[ newpiv_i ]
#              while jnd[ newpiv_j ] < maxj
#                #iter designates the pointer, at which the i-iterator starts
#                iter = ind[ newpiv_i ]
#                #jter designates the pointer, at which the j-iterator starts
#                jter = ind[ jnd[ newpiv_j ] ]
#                #The condition makes sure, that only the columns up to th pivot
#                #are being summed.
#
#                s[ newpiv_j ] -= _innerprod( iter, ind[ newpiv_i + 1 ] - 2,
#                                          jter, ind[ jnd[ newpiv_j ] + 1 ] - 2,
#                                          jnd,
#                                          s)
#
#                s[ newpiv_j ] /= ( s[ ind[ jnd[ newpiv_j ] + 1 ] - 1 ] )
#                newpiv_j += 1
#              end
#              while newpiv_j <= ( ind[ newpiv_i + 1 ] - 1 )
#                s[ newpiv_j ] = 0.
#                newpiv_j += 1
#              end
#            end
#            println( "Pivot in i = ", piv_i, " is nonpositive",
#                     " returning low rank approximation" )
#            return 0
#          else
#            s[ piv_j ] = sqrt( s[ piv_j ] )
#          end
#        else
#          println( "Warning: jnd[ piv_j ] > piv_i" )
#        end
#      end
#    end
#  end
#end

#This function computes the in place incomplete Cholesky factorisation using
#an inner product wrapped in a function.
function _el_icholRWHL!( s::Array{Float64},
                          ind::Array{Int64, 1},
                          jnd::Array{Int64, 1},
                          )
  @inbounds begin
    #piv_i goes from 1 to N and designates the number, which is being treated
    for piv_i = 1 : ( size( ind, 1 ) - 1 )
      #piv_j ranges over pointers and designates where to find the present column
      for piv_j = ind[ piv_i ] : ( ind[ piv_i + 1 ] - 1 )
        #iter designates the pointer, at which the i-iterator starts
        iter::Int64 = ind[ piv_i ]
        #jter designates the pointer, at which the j-iterator starts
        jter::Int64 = ind[ jnd[ piv_j ] ]
        #The condition makes sure, that only the columns up to th pivot
        #are being summed.

        s[ piv_j ] -= _innerprod( iter, ind[ piv_i + 1 ] - 2,
                                  jter, ind[ jnd[ piv_j ] + 1 ] - 2,
                                  jnd,
                                  s)

        if ( s[ ind[ jnd[ piv_j ] + 1 ] - 1 ] ) > 0.0
          if jnd[ piv_j ] < piv_i
            s[ piv_j ] /= ( s[ ind[ jnd[ piv_j ] + 1 ] - 1 ] )
            #debugging:
            if isnan( s[ piv_j ] ) || isinf( s[ piv_j ] )
              println( "ALERT! Not a number" )
            end
            #end debugging

          else
            #debugging:
            if ind[ piv_i + 1 ] - 1 != piv_j
              println( "ALERT! Not selecting diagonal element" )
            end
            #end debugging
            s[ piv_j ] = sqrt( s[ piv_j ] )
            #debugging:
            if isnan( s[ piv_j ] ) || isinf( s[ piv_j ] )
              println( "ALERT! Not a number" )
            end
            #end debugging

          end
        else
          s[ piv_j ] = 0.0
        end
      end
    end
  end
end

#In-place Cholesky factorization of an uppper triangular matrix
function icholU_high_level!( U::SparseMatrixCSC{Float64, Int64} )
  @fastmath _el_icholRWHL!( U.nzval, U.colptr, U.rowval )
end

function icholGPU!( U::SparseMatrixCSC{Float64,Int64} )
  #prevent loss of structural zeros:
  U.nzval .+= eps(Float32)
  d_U = CuArrays.CUSPARSE.CuSparseMatrixCSC( U )
  CuArrays.CUSPARSE.ic02!(d_U, 'O' );

  U.nzval .= Array{Float64}( d_U.nzVal )
end

# probably not needed
# function UfromLTL(σ, L)
#     U = sparse(L')
#     for j = 1 : size(U,2)
#       for ind = U.colptr[j] : (U.colptr[j+1] - 1)
#         i = U.rowval[ind]
#         U.nzval[ind] = L[:,i]' *invσ * L[:,j]
#         if i == j
#           U.nzval[ind] += 1
#         end
#       end
#     end
#     return U
# end

struct IChol
  L
  U
  # optional additive diagonal
invσ
  tmp
end

function IChol(L, U)
  return IChol(L, U, zeros(eltype(L), size(L,1)), Vector{eltype(L)}(undef, size(L,1)))
end

function IChol(L, U, invσ)
  return IChol(L, U, invσ, Vector{eltype(L)}(undef, size(L,1)))
end

import Base.\
import Base.*
import LinearAlgebra.ldiv!
import LinearAlgebra.lmul!
import LinearAlgebra.mul!
import LinearAlgebra.size
import LinearAlgebra.eltype
# Since ldiv! presently does not support sparse triangular matrices,
# the optimized non-allocating versions have been replaced with allocating
# versions (which dispatch to BLAS calls).

function size(A::IChol)
  return size(A.L)
end

function size(A::IChol, dim)
  return size(A.L, dim)
end

function eltype(A::IChol)
  return eltype(A.L)
end

function ldiv!(y, A::IChol, x)
  # ldiv!(A.tmp, A.L, x)
  # ldiv!(y, A.U, A.tmp)
  y .= A.U \ (A.L \ x)
end

function ldiv!(A::IChol, x)
  # ldiv!(A.tmp, A.L, x)
  # ldiv!(x, A.U, A.tmp)
  x .= A.U \ (A.L \ x)
end

function \(A::IChol, x)
  # ldiv!(A.tmp, A.L, x)
  # return A.U \ A.tmp
  return A.U \ (A.L \ x)
end

function *(A::IChol, x)
  # y = A.U * x
  # lmul!(A.L, y)
  # y .+= A.invσ .* x
  #return y
  return A.L * (A.U * x) .+ A.invσ .* x
end

function mul!(y, A::IChol, x)
  # y .= x
  # lmul!(A.U, y)
  # lmul!(A.L, y)
  # y .+= A.invσ .* x
  y .= A.L * (A.U * x) .+ A.invσ .* x
end

function lmul!(A::IChol, x)
  # A.tmp .= A.invσ .* x
  # lmul!(A.U, x)
  # lmul!(A.L, x)
  # x .+= A.tmp
  y = A.L * (A.U * x)
  y = y + A.invσ .* x
  x .= y
end

# in place squaring of a lower triangular matrix to an upper triangular matrix

function squareSparse(L)
  UTemp = sparse(L')
  U = similar(UTemp)
  U.nzval .= UTemp.nzval
  for i = 1 : U.n
    for j in nzrange(U, i)
      rowj = U.rowval[j]
      U.nzval[j] = _innerprod(nzrange(UTemp,i), nzrange(UTemp,rowj), UTemp.rowval, UTemp.nzval)
    end
  end
  return U
end

struct NoiseCov
  L
  U
  invσ

  LICholσ
  LIChol
  UIChol

  function NoiseCov(L, U,σ)
    LICholσ = IChol(L, L', 1 ./ σ)
    LIChol = IChol(L, L')
    UIChol = IChol(U', U)
    return new(L, U, 1 ./ σ, LICholσ, LIChol, UIChol)
  end
end

function ldiv!(y, A::NoiseCov, x)
  y .= A.invσ .* cg(A.LICholσ, A.LIChol * x, Pl=A.UIChol, maxiter=5, tol=1e-16, verbose=true)
end

function ldiv!(A::NoiseCov, x)
  x .= A.invσ .* cg(A.LICholσ, A.LIChol * x, Pl=A.UIChol, maxiter=5, tol=1e-16, verbose=true)
end

function \(A::NoiseCov, x)
  return A.invσ .* cg(A.LICholσ, A.LIChol * x, Pl=A.UIChol, maxiter=5, tol=1e-16, verbose=true)
end

function *(A::NoiseCov, x)
  return A.LIChol \ (A.LICholσ *( A.invσ .\ x))
end

function lmul!(y, A::NoiseCov, x)
  lmul!(y, A.LICholσ, A.invσ .\ x)
  y .*= A.LIChol \ y
end

function lmul!(A::NoiseCov, x)
  lmul!(A.LICholσ, A.invσ .\ x)
  x .= A.LIChol \ x
end

import LinearAlgebra.logdet
function logdet(A::NoiseCov)
  return 2 * logdet(UpperTriangular(A.U)) - 2 * logdet(UpperTriangular(A.L)) + sum(log.(A.invσ))
end