text:
export rhino

### PLugin
#=
Whenever the plugin is updated, run this function and commit the plugin files.
upgrade_plugin()
=#

julia_khepri = dirname(dirname(abspath(@__FILE__)))

dlls = ["KhepriBase.dll", "KhepriRhinoceros.rhp"]
local_plugin = joinpath(dirname(dirname(abspath(@__FILE__))), "Plugin")

upgrade_plugin() =
  let # 1. The dlls are updated in VisualStudio after compilation of the plugin, and they are stored in the folder.
      # contained inside the Plugins folder, which has a specific location regarding this file itself
      plugin_folder = joinpath(dirname(dirname(dirname(dirname(abspath(@__FILE__))))), "Plugins", "KhepriRhinoceros", "KhepriRhinoceros", "bin")
      # 2. The bundle needs to be copied to the current folder
      local_folder = joinpath(julia_khepri, "Plugin")
      # 3. Now we copy the dlls to the local folder
      for dll in dlls
          src = joinpath(plugin_folder, dll)
          dst = joinpath(local_folder, dll)
          rm(dst, force=true)
          cp(src, dst)
      end
  end

##############################################
# We need some additional Encoders
# RH is a subtype of CS
parse_signature(::Val{:RH}, sig::T) where {T} = parse_signature(Val(:CS), sig)
encode(::Val{:RH}, t::Val{T}, c::IO, v) where {T} = encode(Val(:CS), t, c, v)
decode(::Val{:RH}, t::Val{T}, c::IO) where {T} = decode(Val(:CS), t, c)

# We need some additional Encoders
@encode_decode_as(:RH, Val{:Brep}, Val{:Guid})
@encode_decode_as(:RH, Val{:RhinoObject}, Val{:Guid})

encode(::Val{:RH}, t::Union{Val{:Point3d},Val{:Vector3d}}, c::IO, p) =
  encode(Val(:CS), Val(:double3), c, raw_point(p))
decode(::Val{:RH}, t::Val{:Point3d}, c::IO) =
  xyz(decode(Val(:CS), Val(:double3), c)..., world_cs)
decode(::Val{:RH}, t::Val{:Vector3d}, c::IO) =
  vxyz(decode(Val(:CS), Val(:double3), c)..., world_cs)

encode(ns::Val{:RH}, t::Val{:Plane}, c::IO, p) =
  let o = in_world(p),
      x = in_world(vx(1, p.cs)),
      y = in_world(vy(1, p.cs))
    encode(ns, Val(:Point3d), c, o)
    encode(ns, Val(:Vector3d), c, x)
    encode(ns, Val(:Vector3d), c, y)
  end
decode(ns::Val{:RH}, t::Val{:Plane}, c::IO) =
  loc_from_o_vx_vy(
    decode(ns, Val(:Point3d), c),
    decode(ns, Val(:Vector3d), c),
    decode(ns, Val(:Vector3d), c))

rhino_api = @remote_functions :RH """
public void RunScript(string script)
public void SetView(Point3d position, Point3d target, double lens, bool perspective, string style)
public void View(Point3d position, Point3d target, double lens)
public void ViewTop()
public Point3d ViewCamera()
public Point3d ViewTarget()
public double ViewLens()
public string DefaultMaterialsFolder()
public MatId LoadRenderMaterialFromPath(string path)
public MatId CreateMaterial(string name, Color diffuse, Color specular, Color ambient, Color emission, Color reflection, double reflectivity, double reflectionGlossiness, double refractionIndex, double transparency)
public MatId CreateColorMaterial(Color color)
public Guid Point(Point3d p)
public Point3d PointPosition(Guid ent)
public Guid PolyLine(Point3d[] pts)
public Point3d[] LineVertices(RhinoObject id)
public Guid Spline(Point3d[] pts)
public Guid SplineTangents(Point3d[] pts, Vector3d start, Vector3d end)
public Guid ClosedPolyLine(Point3d[] pts)
public Guid ClosedSpline(Point3d[] pts)
public Guid Circle(Point3d c, Vector3d n, double r)
public Point3d CircleCenter(RhinoObject obj)
public Vector3d CircleNormal(RhinoObject obj)
public double CircleRadius(RhinoObject obj)
public Guid Ellipse(Point3d c, Vector3d n, double radiusX, double radiusY)
public Guid Arc(Point3d c, Vector3d vx, Vector3d vy, double radius, double startAngle, double endAngle)
public Guid Text(string str, Point3d c, Vector3d vx, Vector3d vy, double height)
public Guid JoinCurves(Guid[] objs)
public Guid Mesh(Point3d[] pts, int[][] faces, MatId mat)
public Guid NGon(Point3d[] pts, Point3d pivot, bool smooth, MatId mat)
public Guid SurfaceFrom(Guid[] ids, MatId mat)
public Guid SurfaceFromCurvesPoints(Point3d[][] ptss, bool[] smooths, MatId mat)
public Guid SurfaceCircle(Point3d c, Vector3d n, double r, MatId mat)
public Guid SurfaceEllipse(Point3d c, Vector3d n, double radiusX, double radiusY)
public Guid SurfaceArc(Point3d c, Vector3d vx, Vector3d vy, double radius, double startAngle, double endAngle, MatId mat)
public Guid SurfaceClosedPolyLine(Point3d[] pts, MatId mat)
public Guid Sphere(Point3d c, double r, MatId mat)
public Guid Torus(Point3d c, Vector3d vz, double majorRadius, double minorRadius, MatId mat)
public Guid Cylinder(Point3d bottom, double radius, Point3d top, MatId mat)
public Guid Cone(Point3d bottom, double radius, Point3d top, MatId mat)
public Guid ConeFrustum(Point3d bottom, double bottom_radius, Point3d top, double top_radius, MatId mat)
public Guid Box(Point3d corner, Vector3d vx, Vector3d vy, double dx, double dy, double dz, MatId mat)
public Guid XYCenteredBox(Point3d corner, Vector3d vx, Vector3d vy, double dx, double dy, double dz, MatId mat)
public Guid IrregularPyramid(Point3d[] pts, Point3d apex, MatId mat)
public Guid IrregularPyramidFrustum(Point3d[] bpts, Point3d[] tpts, MatId mat)
public Guid PrismWithHoles(Point3d[][] ptss, bool[] smooths, Vector3d dir, MatId mat)
public Guid SurfaceFromGrid(int nU, int nV, Point3d[] pts, bool closedU, bool closedV, int degreeU, int degreeV, MatId mat)
public Brep[] Thicken(RhinoObject obj, double thickness)
public double[] CurveDomain(RhinoObject obj)
public double CurveLength(RhinoObject obj)
public Plane CurveFrameAt(RhinoObject obj, double t)
public Plane CurveFrameAtLength(RhinoObject obj, double l)
public double[] SurfaceDomain(RhinoObject obj)
public Plane SurfaceFrameAt(RhinoObject obj, double u, double v)
public Point3d[] BoundingBox(ObjectId[] ids)
public Brep Extrusion(RhinoObject obj, Vector3d dir)
public Brep[] SweepPathProfile(RhinoObject path, RhinoObject profile, double rotation, double scale)
public Guid[] Intersect(RhinoObject obj0, RhinoObject obj1)
public Guid[] Subtract(RhinoObject obj0, RhinoObject obj1)
public Guid Move(Guid id, Vector3d v)
public Guid Scale(Guid id, Point3d p, double s)
public Guid Rotate(Guid id, Point3d p, Vector3d n, double a)
public Guid Mirror(Guid id, Point3d p, Vector3d n, bool copy)
public void Render(int width, int height, string path)
public void SaveView(int width, int height, string path)
public bool IsPoint(RhinoObject e)
public bool IsCircle(RhinoObject e)
public bool IsPolyLine(RhinoObject e)
public bool IsSpline(RhinoObject e)
public bool IsClosedPolyLine(RhinoObject e)
public bool IsClosedSpline(RhinoObject e)
public bool IsEllipse(RhinoObject e)
public bool IsArc(RhinoObject e)
public bool IsText(RhinoObject e)
public byte ShapeCode(RhinoObject id)
public Guid ClosedPathCurveArray(Point3d[] pts, double[] angles)
public Brep[] PathWall(RhinoObject obj, double lThickness, double rThickness, double height)
public Brep RectangularTable(Point3d c, double angle, double length, double width, double height, double top_thickness, double leg_thickness)
public Brep[] BaseRectangularTable(double length, double width, double height, double top_thickness, double leg_thickness)
public int CreateRectangularTableFamily(double length, double width, double height, double top_thickness, double leg_thickness)
public Guid Table(Point3d c, double angle, int family)
public Brep[] BaseChair(double length, double width, double height, double seat_height, double thickness)
public int CreateChairFamily(double length, double width, double height, double seat_height, double thickness)
public Guid Chair(Point3d c, double angle, int family)
public GeometryBase[] BaseRectangularTableAndChairs(int tableFamily, int chairFamily, double tableLength, double tableWidth, int chairsOnTop, int chairsOnBottom, int chairsOnRight, int chairsOnLeft, double spacing)
public int CreateRectangularTableAndChairsFamily(int tableFamily, int chairFamily, double tableLength, double tableWidth, int chairsOnTop, int chairsOnBottom, int chairsOnRight, int chairsOnLeft, double spacing)
public Guid TableAndChairs(Point3d c, double angle, int family)
public Guid CreateLayer(String name, bool active, byte r, byte g, byte b)
public Guid CurrentLayer()
public void SetCurrentLayer(Guid id)
public void SetLayerActive(Guid id, bool active) {
public Guid ShapeLayer(RhinoObject objId)
public void SetShapeLayer(RhinoObject objId, Guid layerId)
public Point3d[] CurvePointsAt(RhinoObject obj, double[] ts)
public Vector3d[] CurveTangentsAt(RhinoObject obj, double[] ts)
public Vector3d[] CurveNormalsAt(RhinoObject obj, double[] ts)
public int DeleteAll()
public int DeleteAllInLayer(Guid id)
public void Delete(Guid id)
public void DeleteMany(Guid[] ids)
public void Select(Guid id)
public void SelectMany(Guid[] ids)
public void Highlight(RhinoObject o)
public void Unhighlight(RhinoObject o)
public void UnhighlightAll()
public Guid Clone(RhinoObject e)
public Point3d[] GetPosition(string prompt)
public Guid[] GetPoint(string prompt)
public Guid[] GetPoints(string prompt)
public Guid[] GetCurve(string prompt)
public Guid[] GetCurves(string prompt)
public Guid[] GetSurface(string prompt)
public Guid[] GetSurfaces(string prompt)
public Guid[] GetSolid(string prompt)
public Guid[] GetSolids(string prompt)
public Guid[] GetShape(string prompt)
public Guid[] GetShapes(string prompt)
public Guid[] GetAllShapes()
public Guid[] GetAllShapesInLayer(Guid id)
public void SetSun(DateTime dt, double latitude, double longitude, int meridian)
public void EnableGrasshopperSolver()
public void DisableGrasshopperSolver()
public void RunGrasshopperSolver()
"""
#=
public byte Sync()
public byte Disconnect()
public Entity InterpSpline(Point3d[] pts, Vector3d tan0, Vector3d tan1)
public Entity InterpClosedSpline(Point3d[] pts)
public Guid SurfaceFromCurve(Entity curve)
public Guid IrregularPyramidMesh(Point3d[] pts, Point3d apex)
public Entity SolidFromGrid(int m, int n, Point3d[] pts, bool closedM, bool closedN, int level, double thickness)
public Guid NurbSurfaceFrom(ObjectId id)
public double[] SurfaceDomain(Entity ent)
public Frame3d SurfaceFrameAt(Entity ent, double u, double v)
public Guid Sweep(ObjectId pathId, ObjectId profileId, double rotation, double scale)
public Guid Loft(ObjectId[] profilesIds, ObjectId[] guidesIds, bool ruled, bool closed)
public void Unite(ObjectId objId0, ObjectId objId1)
public Guid Revolve(ObjectId profileId, Point3d p, Vector3d n, double startAngle, double amplitude)
public Point3d[] GetPoint(string prompt)
public void ZoomExtents()
public void SetSystemVariableInt(string name, int value)
public int Command(string cmd)
public void DisableUpdate()
public void EnableUpdate()
public BIMLevel FindOrCreateLevelAtElevation(double elevation)
public BIMLevel UpperLevel(BIMLevel currentLevel, double addedElevation)
public double GetLevelElevation(BIMLevel level)
public Entity MeshFromGrid(int m, int n, Point3d[] pts, bool closedM, bool closedN)

public FloorFamily FloorFamilyInstance(double totalThickness, double coatingThickness)
public Entity LightweightPolyLine(Point2d[] pts, double[] angles, double elevation)
public Entity SurfaceLightweightPolyLine(Point2d[] pts, double[] angles, double elevation)
public ObjectId CreatePathFloor(Point2d[] pts, double[] angles, BIMLevel level, FloorFamily family)
=#

abstract type RHKey end
const RHId = Guid
const RHRef = GenericRef{RHKey, RHId}
const RHRefs = Vector{RHRef}
const RHEmptyRef = EmptyRef{RHKey, RHId}
const RHNativeRef = NativeRef{RHKey, RHId}
const RHUnionRef = UnionRef{RHKey, RHId}
const RHSubtractionRef = SubtractionRef{RHKey, RHId}
const RH = SocketBackend{RHKey, RHId}

KhepriBase.void_ref(::RH) = 0 % UInt128

KhepriBase.failed_connecting(b::RH) =
  begin
	@info("Please, ensure Khepri's Rhinoceros plugin is installed.")
  	@info("Go to Rhinoceros: Tools->Options->Plug-ins->Install")
  	@info("Location: $(joinpath(local_plugin, dlls[2]))")
  end

KhepriBase.after_connecting(b::RH) =
  let mat_folder = @remote(b, DefaultMaterialsFolder()),
  	  m = match(r".+?Rhinoceros\\(.+?)\\Localization", mat_folder)
  	if m === nothing
	  error("Unrecognized Rhino materials folder: $(mat_folder)")
  	else
	  let v = VersionNumber(m[1])
	    if v >= v"6.0"
			# Apparently, materials are not yet working in Rhino 6!!!!
	      set_material(rhino, material_metal, rhino_default_material(raw"Metal\Matte\Matte Silver"))
	      set_material(rhino, material_glass, rhino_default_material(raw"Glass\Clear Glass"))
	      set_material(rhino, material_wood, rhino_default_material(raw"Wood\Pear polished"))
	      set_material(rhino, material_concrete, rhino_default_material(raw"Ceramics\Stoneware"))
	      set_material(rhino, material_plaster, rhino_default_material(raw"White Matte"))
	      set_material(rhino, material_grass, rhino_default_material(raw"Textures\Grass"))
	    elseif v >= v"5.0"
          set_material(rhino, material_metal, rhino_default_material(raw"Metal\Silver"))
          set_material(rhino, material_glass, rhino_default_material(raw"Transparent\Glass"))
          set_material(rhino, material_wood, rhino_default_material(raw"Wood\Scots Pine"))
          set_material(rhino, material_concrete, rhino_default_material(raw"Stone\Granite"))
          set_material(rhino, material_plaster, rhino_default_material(raw"White Matte"))
          set_material(rhino, material_grass, rhino_default_material(raw"Special\Grass"))
	    else
	      error("Unrecognized Rhino version: ")
	    end
	  end
    end
  end

const rhino = RH("Rhino", rhino_port, rhino_api)

KhepriBase.has_boolean_ops(::Type{RH}) = HasBooleanOps{false}()
KhepriBase.backend(::RHRef) = rhino

# Primitives
KhepriBase.b_point(b::RH, p) =
  @remote(b, Point(p))

KhepriBase.b_line(b::RH, ps, mat) =
  @remote(b, PolyLine(ps))

KhepriBase.b_polygon(b::RH, ps, mat) =
  @remote(b, ClosedPolyLine(ps))

KhepriBase.b_spline(b::RH, ps, v0, v1, interpolator, mat) =
  if (v0 == false) && (v1 == false)
    @remote(b, Spline(ps))
  elseif (v0 != false) && (v1 != false)
    @remote(b, SplineTangents(ps, v0, v1))
  else
    @remote(b, SplineTangents(ps,
                 v0 == false ? ps[2]-ps[1] : v0,
                 v1 == false ? ps[end-1]-ps[end] : v1))
  end

KhepriBase.b_closed_spline(b::RH, ps, mat) =
  @remote(b, ClosedSpline(ps))

KhepriBase.b_circle(b::RH, c, r, mat) =
  @remote(b, Circle(c, vz(1, c.cs), r))

KhepriBase.b_arc(b::RH, c, r, α, Δα, mat) =
  if r == 0
    @remote(b, Point(c))
  elseif Δα == 0
    @remote(b, Point(c + vpol(r, α, c.cs)))
  elseif abs(Δα) >= 2*pi
    @remote(b, Circle(c, vz(1, c.cs), r))
  else
	let β = α + amplitude
  	  if β > α
  	  	@remote(b, Arc(c, vx(1, c.cs), vy(1, c.cs), r, α, β))
  	  else
  	  	@remote(b, Arc(c, vx(1, c.cs), vy(1, c.cs), r, β, α))
  	  end
    end
  end

b_ellipse() =
  @remote(b, Ellipse(s.center, vz(1, s.center.cs), s.radius_x, s.radius_y))

KhepriBase.b_trig(b::RH, p1, p2, p3, mat) =
  @remote(b, Mesh([p1, p2, p3], [[0, 1, 2, 2]], mat))

KhepriBase.b_quad(b::RH, p1, p2, p3, p4, mat) =
 	@remote(b, Mesh([p1, p2, p3, p4], [[0, 1, 2, 2], [0, 2, 3, 3]], mat))

KhepriBase.b_ngon(b::RH, ps, pivot, smooth, mat) =
 	@remote(b, NGon(ps, pivot, smooth, mat))

############################################################
# Second tier: surfaces

KhepriBase.b_surface_polygon(b::RH, ps, mat) =
  @remote(b, SurfaceClosedPolyLine([ps...,ps[1]], mat))

KhepriBase.b_surface_polygon_with_holes(b::RH, ps, qss, mat) =
  @remote(b, SurfaceFromCurvesPoints([ps, qss...], falses(1 + length(qss)), mat))

KhepriBase.b_surface_circle(b::RH, c, r, mat) =
	@remote(b, SurfaceCircle(c, vz(1, c.cs), r, mat))

KhepriBase.b_surface_arc(b::RH, c, r, α, Δα, mat) =
  @remote(b, SurfaceArc(c, vx(1, c.cs), vy(1, c.cs), r, α, α + Δα, mat))

#=
KhepriBase.b_quad(b::RH, p1, p2, p3, p4, mat) =
 [b_trig(b, p1, p2, p3, mat),
  b_trig(b, p1, p3, p4, mat)]

KhepriBase.b_quad_strip(b::RH, ps, qs, smooth, mat) =
 [b_quad(b, ps[i], ps[i+1], qs[i+1], qs[i], mat)
  for i in 1:size(ps,1)-1]

KhepriBase.b_quad_strip_closed(b::RH, ps, qs, smooth, mat) =
 b_quad_strip(b, [ps..., ps[1]], [qs..., qs[1]], smooth, mat)

KhepriBase.b_surface_grid(b::RH, ptss, closed_u, closed_v, smooth_u, smooth_v, interpolator, mat) =
  let (nu, nv) = size(ptss)
 if smooth_u && smooth_v
  	ptss = [location_at_grid_interpolator(interpolator, u, v)
  			for u in division(0, 1, 2 #=4*nu-7=#), v in division(0, 1, 4*nv-7)] # 2->1, 3->5, 4->9, 5->13
  	(nu, nv) = size(ptss)
 end
 closed_u ?
  	vcat([b_quad_strip_closed(b, ptss[i,:], ptss[i+1,:], smooth_u, mat) for i in 1:nu-1],
  		 closed_v ? [b_quad_strip_closed(b, ptss[end,:], ptss[1,:], smooth_u, mat)] : []) :
  	vcat([b_quad_strip(b, ptss[i,:], ptss[i+1,:], smooth_u, mat) for i in 1:nu-1],
  		 closed_v ? [b_quad_strip(b, ptss[end,:], ptss[1,:], smooth_u, mat)] : [])
  end
=#
KhepriBase.b_surface_grid(b::RH, ptss, closed_u, closed_v, smooth_u, smooth_v, mat) =
  let (nu, nv) = size(ptss),
      order(n) = min(2*ceil(Int,n/16) + 1, 5)
    @remote(b, SurfaceFromGrid(nu, nv,
                               reshape(permutedims(ptss), :),
                               closed_u, closed_v,
                               smooth_u ? order(nu) : 1,
                               smooth_v ? order(nv) : 1,
							   mat))
  end

KhepriBase.b_surface_mesh(b::RH, vertices, faces, mat) =
  let ensure_4(v) = length(v) == 4 ? v : [v..., v[3]]
    @remote(b, Mesh(vertices, map(face->ensure_4(face.-1), faces), mat))
  end

# This is wrong, for sure!
b_surface_ellipse(b::RH, c, rx, ry) =
   @remote(b, SurfaceEllipse(c, vz(1, c.cs), rx, ry))

############################################################
# Third tier: solids

KhepriBase.b_generic_prism(b::RH, bs, smooth, v, bmat, tmat, smat) =
	@remote(b, PrismWithHoles([bs], [smooth], v, tmat))

KhepriBase.b_generic_prism_with_holes(b::RH, bs, smooth, bss, smooths, v, bmat, tmat, smat) =
  @remote(b, PrismWithHoles([bs, bss...], [smooth, smooths...], v, tmat))

KhepriBase.b_pyramid_frustum(b::RH, bs, ts, bmat, tmat, smat) =
  @remote(b, IrregularPyramidFrustum(bs, ts, tmat))

KhepriBase.b_pyramid(b::RH, bs, t, bmat, smat) =
  @remote(b, IrregularPyramid(bs, t, smat))

KhepriBase.b_cylinder(b::RH, c, r, h, bmat, tmat, smat) =
  @remote(b, Cylinder(c, r, c + vz(h, c.cs), smat))

KhepriBase.b_box(b::RH, c, dx, dy, dz, mat) =
  @remote(b, Box(c, vx(1, c.cs), vy(1, c.cs), dx, dy, dz, mat))

KhepriBase.b_sphere(b::RH, c, r, mat) =
	@remote(b, Sphere(c, r, mat))

KhepriBase.b_cone(b::RH, cb, r, h, bmat, smat) =
  @remote(b, Cone(add_z(cb, h), r, cb, smat))

KhepriBase.b_cone_frustum(b::RH, cb, rb, h, rt, bmat, tmat, smat) =
  @remote(b, ConeFrustum(cb, rb, cb + vz(h, cb.cs), rt, smat))

KhepriBase.b_text(b::RH, str, c, size, mat) =
  @remote(b, Text(str, c, vx(1, c.cs), vy(1, c.cs), size))

## Materials
# To encode materials, we use a Guid but convert it to an int,
# taking into consideration that Guid(int) = int + 1
encode(ns::Val{:RH}, t::Val{:MatId}, c::IO, v) =
  encode(ns, Val(:int), c, v%Int32 - 1)
decode(ns::Val{:RH}, t::Val{:MatId}, c::IO) =
  (decode(ns, Val(:int), c) + 1)%Guid

KhepriBase.b_get_material(b::RH, ref) =
  get_rhino_material(b, ref)

get_rhino_material(b, n::Nothing) =
  void_ref(b)
get_rhino_material(b, path::String) =
  @remote(b, LoadRenderMaterialFromPath(path))

struct RhinoDefaultMaterial
  name::String
end

export rhino_default_material
rhino_default_material = RhinoDefaultMaterial

get_rhino_material(b, mat::RhinoDefaultMaterial) =
  get_rhino_material(b, joinpath(@remote(b, DefaultMaterialsFolder()), mat.name*".rmtl"))

KhepriBase.b_new_material(b::RH, path, color, specularity, roughness, transmissivity, transmitted_specular) =
  Guid(@remote(b, new_material(path, convert(RGBA, color), specularity, roughness)) + 1)

realize(b::RH, s::EmptyShape) =
  RHEmptyRef()
realize(b::RH, s::UniversalShape) =
  RHUniversalRef()

backend_map_division(b::RH, f::Function, s::Shape1D, n::Int) =
  let conn = connection(b),
      r = ref(s).value,
      (t1, t2) = @remote(b, CurveDomain(r)),
      ti = division(t1, t2, n),
      ps = @remote(b, CurvePointsAt(r, ti)),
      ts = @remote(b, CurveTangentsAt(r, ti)),
      #ns = ACADCurveNormalsAt(conn, r, ti),
      frames = rotation_minimizing_frames(@remote(b, CurveFrameAt(r, t1)), ps, ts)
    map(f, frames)
  end

realize(b::RH, s::Surface) =
  let ids = @remote(b, SurfaceFrom(collect_ref(s.frontier)))
    foreach(mark_deleted, s.frontier)
    ids
  end
#=
realize(b::RH, s::Text) =
  @remote(b, Text(s.str, s.c, vx(1, s.c.cs), vy(1, s.c.cs), s.h))
=#

backend_surface_domain(b::RH, s::Shape2D) =
    tuple(@remote(b, SurfaceDomain(ref(s).value)...))

backend_map_division(b::RH, f::Function, s::Shape2D, nu::Int, nv::Int) =
    let conn = connection(b)
        r = ref(s).value
        (u1, u2, v1, v2) = @remote(b, SurfaceDomain(r))
        map_division(u1, u2, nu) do u
            map_division(v1, v2, nv) do v
                f(@remote(b, SurfaceFrameAt(r, u, v)))
            end
        end
    end

KhepriBase.b_torus(b::RH, c, ra, rb, mat) =
  @remote(b, Torus(c, vz(1, c.cs), ra, ra, mat))

backend_extrusion(b::RH, s::Shape, v::Vec) =
    and_mark_deleted(
        map_ref(s) do r
            @remote(b, Extrusion(r, v))
        end,
        s)

KhepriBase.b_sweep(b::RH, path, profile, rotation, scale, mat) =
  map_ref(b, profile) do profile_r
    map_ref(b, path) do path_r
      @remote(b, SweepPathProfile(path_r, profile_r, rotation, scale))
    end
  end

#=
realize(b::RH, s::Revolve) =
  and_delete_shape(
    map_ref(s.profile) do r
      @remote(b, Revolve(r, s.p, s.n, s.start_angle, s.amplitude))
    end,
    s.profile)

backend_loft_points(b::Backend, profiles::Shapes, rails::Shapes, ruled::Bool, closed::Bool) =
  let f = (ruled ? (closed ? polygon : line) : (closed ? closed_spline : spline))
    and_delete_shapes(ref(f(map(point_position, profiles), backend=b)),
                      vcat(profiles, rails))
  end

backend_loft_curves(b::RH, profiles::Shapes, rails::Shapes, ruled::Bool, closed::Bool) =
  and_delete_shapes(RHLoft(connection(b),
                             collect_ref(profiles),
                             collect_ref(rails),
                             ruled, closed),
                    vcat(profiles, rails))

backend_loft_surfaces(b::RH, profiles::Shapes, rails::Shapes, ruled::Bool, closed::Bool) =
  backend_loft_curves(b, profiles, rails, ruled, closed)

backend_loft_curve_point(b::RH, profile::Shape, point::Shape) =
  and_delete_shapes(RHLoft(connection(b),
                             vcat(collect_ref(profile), collect_ref(point)),
                             [],
                             true, false),
                    [profile, point])

backend_loft_surface_point(b::RH, profile::Shape, point::Shape) =
  backend_loft_curve_point(b, profile, point)
=#
unite_refs(b::RH, refs::Vector{<:RHRef}) =
    RHUnionRef(tuple(refs...))

unite_refs(b::RH, r::RHUnionRef) =
  r

unite_ref(b::RH, r0::RHNativeRef, r1::RHNativeRef) =
  RHUnionRef((r0, r1))

intersect_ref(b::RH, r0::RHNativeRef, r1::RHNativeRef) =
    let refs = @remote(b, Intersect(r0.value, r1.value))
        n = length(refs)
        if n == 0
            RHEmptyRef()
        elseif n == 1
            RHNativeRef(refs[1])
        else
            RHUnionRef(map(RHNativeRef, tuple(refs...)))
        end
    end

subtract_ref(b::RH, r0::RHNativeRef, r1::RHNativeRef) =
  let refs = try @remote(b, Subtract(r0.value, r1.value)); catch e; [nothing] end,
      n = length(refs)
    if n == 0
        RHEmptyRef()
    elseif n == 1
        if isnothing(refs[1]) # failed
            RHSubtractionRef(r0, tuple(r1))
        else
            RHNativeRef(refs[1])
        end
    else
        RHUnionRef(map(RHNativeRef, tuple(refs...)))
    end
  end

subtract_ref(b::RH, r0::SubtractionRef, r1::RHNativeRef) =
  subtract_ref(b,
    subtract_ref(b, r0.value, r1),
    length(r0.values) == 1 ? r0.values[1] : RHUnionRef(r0.values))

subtract_ref(b::RH, r0::RHRef, r1::RHUnionRef) =
  foldl((r0,r1)->subtract_ref(b,r0,r1), r1.values, init=r0)

#=
slice_ref(b::RH, r::RHNativeRef, p::Loc, v::Vec) =
  (@remote(b, Slice(r.value, p, v); r))

slice_ref(b::RH, r::RHUnionRef, p::Loc, v::Vec) =
  map(r->slice_ref(b, r, p, v), r.values)

=#

#=
realize(b::RH, s::IntersectionShape) =
  foldl((r0,r1)->intersect_ref(b,r0,r1), UniversalRef{RHId}(), map(ref, s.shapes))
  =#

realize(b::RH, s::Slice) =
  slice_ref(b, ref(s.shape), s.p, s.n)

realize(b::RH, s::Move) =
  and_mark_deleted(b,
	map_ref(b, s.shape) do r
      @remote(b, Move(r, s.v))
      r
  end,
  s.shape)

realize(b::RH, s::Scale) =
  and_mark_deleted(b,
  	map_ref(b, s.shape) do r
      @remote(b, Scale(r, s.p, s.s))
      r
    end,
	s.shape)

realize(b::RH, s::Rotate) =
  and_mark_deleted(b,
  	map_ref(b, s.shape) do r
      @remote(b, Rotate(r, s.p, s.v, s.angle))
      r
	end,
	s.shape)

realize(b::RH, s::Mirror) =
  and_delete_shape(map_ref(b, s.shape) do r
                    @remote(b, Mirror(r, s.p, s.n, false))
                   end,
                   s.shape)

realize(b::RH, s::UnionMirror) =
  let r0 = ref(s.shape),
      r1 = map_ref(b, r0) do r
            @remote(b, Mirror(r, s.p, s.n, true))
          end
    UnionRef((r0,r1))
  end

realize(b::RH, s::Thicken) =
  and_mark_deleted(b,
    map_ref(b, s.shape) do r
      @remote(b, Thicken(r, s.thickness))
    end,
    s.shape)

backend_frame_at(b::RH, s::Shape2D, u::Real, v::Real) =
  @remote(b, SurfaceFrameAt(ref(b, s).value, u, v))

# BIM

# Families

realize(b::RH, f::TableFamily) =
    @remote(b, CreateRectangularTableFamily(f.length, f.width, f.height, f.top_thickness, f.leg_thickness))
realize(b::RH, f::ChairFamily) =
    @remote(b, CreateChairFamily(f.length, f.width, f.height, f.seat_height, f.thickness))
realize(b::RH, f::TableChairFamily) =
    @remote(b, CreateRectangularTableAndChairsFamily(
        realize(b, f.table_family), realize(b, f.chair_family),
        f.table_family.length, f.table_family.width,
        f.chairs_top, f.chairs_bottom, f.chairs_right, f.chairs_left,
        f.spacing))

KhepriBase.b_table(b::RH, c, angle, family) =
    @remote(b, Table(c, angle, realize(b, family)))

KhepriBase.b_chair(b::RH, c, angle, family) =
    @remote(b, Chair(c, angle, realize(b, family)))

KhepriBase.b_table_and_chairs(b::RH, c, angle, family) =
    @remote(b, TableAndChairs(c, angle, realize(b, family)))

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

# KhepriBase.b_bounding_box(b::RH, shapes::Shapes) =
#   @remote(b, BoundingBox(collect_ref(shapes)))

KhepriBase.b_set_view(b::RH, camera::XYZ, target::XYZ, lens::Real, aperture::Real) =
  @remote(b, View(camera, target, lens))

KhepriBase.b_get_view(b::RH) =
  @remote(b, ViewCamera()), @remote(b, ViewTarget()), @remote(b, ViewLens())

KhepriBase.b_zoom_extents(b::RH) =
  @remote(b, ZoomExtents())

KhepriBase.b_set_view_top(b::RH) =
  @remote(b, ViewTop())

KhepriBase.b_delete_refs(b::RH, refs::Vector{RHId}) =
  @remote(b, DeleteMany(refs))

KhepriBase.b_all_refs(b::RH) =
  @remote(b, GetAllShapes())

KhepriBase.b_delete_all_refs(b::RH) =
  @remote(b, DeleteAll())

KhepriBase.b_highlight_ref(b::RH, r::RHId) =
  @remote(b, Highlight(r))

KhepriBase.b_unhighlight_ref(b::RH, r::RHId) =
  @remote(b, Unhighlight(r))

KhepriBase.b_unhighlight_all_refs(b::RH) =
  @remote(b, UnhighlightAll())

# Layers
KhepriBase.b_current_layer(b::RH) =
  @remote(b, CurrentLayer())

KhepriBase.b_current_layer(b::RH, layer) =
  @remote(b, SetCurrentLayer(layer))

KhepriBase.b_layer(b::RH, name, active, color) =
  let to255(x) = round(UInt8, x*255)
    @remote(b, CreateLayer(name, true, to255(red(color)), to255(green(color)), to255(blue(color))))
  end

KhepriBase.b_delete_all_shapes_in_layer(b::RH, layer) =
  @remote(b, DeleteAllInLayer(layer))

KhepriBase.b_shape_from_ref(b::RH, r) =
  let code = @remote(b, ShapeCode(r)),
      ref = DynRefs(b=>RHNativeRef(r))
    if code == 1
      point(@remote(b, PointPosition(r)),
            backend=b, ref=ref)
    elseif code == 2
      circle(maybe_loc_from_o_vz(@remote(b, CircleCenter(r)), @remote(b, CircleNormal(r))),
             @remote(b, CircleRadius(r)),
             backend=b, ref=ref)
    elseif 3 <= code <= 6
      line(@remote(b, LineVertices(r)), ref=ref)
    elseif code == 7
      spline([xy(0,0)], false, false, #HACK obtain interpolation points
             ref=ref)
    elseif 103 <= code <= 106
      polygon(@remote(b, LineVertices(r)), ref=ref)
    elseif code == 40
      # FIXME: frontier is missing
      surface(frontier=[], ref=ref)
    elseif code == 41
      surface(frontier=[], ref=ref)
    elseif code == 81
      sphere(backend=b, ref=ref)
    elseif code == 82
      cylinder(backend=b, ref=ref)
    elseif code == 83
      cone(backend=b, ref=ref)
    elseif code == 84
      torus(backend=b, ref=ref)
    else
      error("Unknown shape")
    end
  end

#

KhepriBase.b_select_position(b::RH, prompt::String) =
  begin
    @info "$(prompt) on the $(b) backend."
    let ans = @remote(b, GetPosition(prompt))
      length(ans) > 0 ? ans[1] : nothing
    end
  end

KhepriBase.b_select_positions(b::RH, prompt::String) =
  let ps = Loc[]
      p = nothing
    @info "$(prompt) on the $(b) backend."
    while ((p = @remote(b, GetPosition(prompt)) |> length)) > 0
        push!(ps, p...)
    end
    ps
  end

# HACK: The next operations should receive a set of shapes to avoid re-creating already existing shapes

KhepriBase.b_select_point(b::RH, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetPoint)
KhepriBase.b_select_points(b::RH, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetPoints)

KhepriBase.b_select_curve(b::RH, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetCurve)
KhepriBase.b_select_curves(b::RH, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetCurves)

KhepriBase.b_select_surface(b::RH, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetSurface)
KhepriBase.b_select_surfaces(b::RH, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetSurfaces)

KhepriBase.b_select_solid(b::RH, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetSolid)
KhepriBase.b_select_solids(b::RH, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetSolids)

KhepriBase.b_select_shape(b::RH, prompt::String) =
  select_one_with_prompt(prompt, b, @get_remote b GetShape)
KhepriBase.b_select_shapes(b::RH, prompt::String) =
  select_many_with_prompt(prompt, b, @get_remote b GetShapes)

backend_captured_shape(b::RH, handle) =
  backend_shape_from_ref(b, handle)
#
backend_captured_shapes(b::RH, handles) =
  map(handles) do handle
      backend_shape_from_ref(b, handle)
  end

backend_generate_captured_shape(b::RH, s::Shape) =
    println("captured_shape(rhino, $(ref(s).value))")
backend_generate_captured_shapes(b::RH, ss::Shapes) =
  begin
    print("captured_shapes(rhino, [")
    for s in ss
      print(ref(s).value)
      print(", ")
    end
    println("])")
  end

backend_all_shapes_in_layer(b::RH, layer) =
  Shape[backend_shape_from_ref(b, r) for r in @remote(b, GetAllShapesInLayer(layer))]

KhepriBase.b_render_view(b::RH, path::String) =
  @remote(b, Render(render_width(), render_height(), path))

backend_save_view(b::RH, path::String) =
  @remote(b, SaveView(render_width(), render_height(), path))

#=
BIM families for Rhino
=#

abstract type RhinoFamily <: Family end

struct RhinoLayerFamily <: RhinoFamily
  name::String
  color::RGB
  ref::Parameter{Any}
end

rhino_layer_family(name, color::RGB=rgb(1,1,1)) =
  RhinoLayerFamily(name, color, Parameter{Any}(nothing))

backend_get_family_ref(b::RH, f::Family, af::RhinoLayerFamily) =
  backend_create_layer(b, af.name, true, af.color)

KhepriBase.b_realistic_sky(b::RH, date, latitude, longitude, elevation, meridian, turbidity, withsun) =
  @remote(b, SetSun(date, latitude, longitude, meridian))

export run_grasshopper_solver
run_grasshopper_solver() =
  @remote(rhino, RunGrasshopperSolver())


text:
"""
$(TYPEDSIGNATURES)

Resurrect linspace despite https://github.com/JuliaLang/julia/pull/25896#issuecomment-363769368
"""
linspace(a,b,n)=collect(range(a,b,length=n))

"""
$(TYPEDSIGNATURES)

(Try to) create a subdivision of interval (a,b) stored in the
returned array X such that
  - `X[1]==a, X[end]==b`
  - `(X[2]-X[1])<=ha+tol*(b-a)`
  - `(X[end]-X[end-1])<=hb+tol*(b-a)`
  - There is a number q such that  `X[i+1]-X[i] == q*(X[i]-X[i-1])`
  - X is the array with the minimal possible number of points with the above property

Caveat: the algorithm behind this is  tested for many cases but unproven.

Returns an Array containing the points of the subdivision.
"""
function geomspace(a, b, ha, hb ; tol=1.0e-10, maxiterations=100)
    a=Float64(a)
    b=Float64(b)
    ha=Float64(ha)
    hb=Float64(hb)

    function _geomspace0(l,h0, hl, tol=1.0e-10)
        @assert (l>0.0)
        @assert (h0>0.0)
        @assert (hl>=h0)
        @assert((hl+h0)<l)

        #  We need to  adjust two things:

        # The sum of the geometric progression must
        # match the length of the interval, so lmismatch
        # should be zero:
        function lmismatch(q,k)
            return l - h0*(1-q^k)/(1-q)
        end

        # The claim from experimenral evidence (Wolfram)
        # is that, if written as a polynomial,
        # it has two real zeros: one and the value searched for which
        # is slightly larger than one. All other zeros are on one circle.

        # The size of the last interval should be close to
        # to hl, so hmismatch should be close to one and not larger than one
        function  hmismatch(q,k)
            return  h0*q^(k-1)/hl
        end

        # define initial number of intervals from
        # average of minmal and maximal h
        n=Int(ceil((2.0*l/(h0+hl))))

        if n==1
            n=2
        end

        # define initial q such that hmismatch is one.
        q=(hl/h0)^(1.0/(n-1.0))

        # Iteration until both mismatches are satisfactory
        # Outer loop runs until hmismatch is less than 1
        hmiss=10.0 # some initial value >1 just to run the loop at least once
        if abs(q-1.0)<tol
            hmiss=1.0
        end

        while  hmiss>1.0
            # increase number of intervals until
            # lmismatch becomes less than zero
            while  lmismatch(q,n)>tol
                n+=1
            end

            # find initial interval for q containing
            # value with zero lmismatch
            ns=0

            while lmismatch(q,n)<tol &&  ns<maxiterations
                q*=0.99
                ns+=1
            end

            ql=q*0.9
            qr=q*1.1
            @assert ns<maxiterations "Unable to determine geomspace data after $(maxiterations) iterations"

            # bisection to define q with zero lmismatch
            ns=0
            qm=0.5*(ql+qr)
            while abs(qr-ql)>tol && ns<maxiterations
                ns+=1
                mmm=lmismatch(qm,n)
                if abs(mmm)<0
                    break
                elseif lmismatch(ql,n)*mmm<0
                    qr=qm
                else
                    ql=qm
                end
                qm=0.5*(ql+qr)
            end
            # increase q slightly to increase probability
            # for last interval to be <=hl
            q=qm*(1.0+tol)
            @assert ns<maxiterations "Unable to determine geomspace data after $(maxiterations) iterations"
            hmiss=hmismatch(q,n)
            if hmiss>1.0+tol
                n=n+1
            end
        end
        #  printf("%d %g %g %g\n",n,q,lmismatch(q,n),hmismatch(q,n))

        X = Array{Float64,1}(undef,n+1)
        X[1]=0
        h=h0
        for i=1:n
            X[i+1]=X[i]+h
            h*=q
        end
        X[n+1]=l

        # Fix last interval and remove "overshoots"
        if X[end]>l
            X[end]=l
        end

        while X[end-1]+tol>X[end]
            pop!(X)
            X[end]=l
        end

        return X
    end

    @assert (ha>0.0) "Start step size $(ha) should be positive"
    @assert (hb>0.0) "End step size $(hb) should be positive"
    @assert (a<b)    "Interval ends $(a), $(b) should be increasing"
    @assert ((ha+hb)<b-a) "Sum of step sizes $(ha)+$(hb) should not exceed interval size $(b)-$(a)"

    # Map things to [0,b-a]
    tol=tol*(b-a)
    if ha<hb-tol
        X=_geomspace0(b-a,ha,hb,tol)
        X.+=a
    elseif ha>hb+tol
        X=-reverse(_geomspace0(b-a,hb,ha,tol))
        X.+=b
    else
        n=Int(ceil((b-a)/ha))
        X=collect(range(a,b,length=n+1))
    end
#    @show X[2]-X[1],ha, X[end]-X[end-1],hb

    @assert (X[2]-X[1])<=ha+tol  "First interval turned out to be larger than $(ha)"
    @assert (X[2]-X[1])>ha/2  "First interval turned out to be less  than $(ha)/10"
    @assert (X[end]-X[end-1])<=hb+tol  "Last interval turned out to be larger than $(hb)"
    @assert (X[end]-X[end-1])>=hb/2  "Last interval turned out to be less than $(hb)/10"
    @assert abs(X[1]-a)<tol  "Range start $(X[1]) differs from $(a)"
    @assert abs(X[end]-b)<tol "Range end $(X[end])  differs from $(b)"

    return X
end

collect_or_assign(X)=X
collect_or_assign(X::AbstractRange)=collect(X)
is_monotone(X)=all(X[2:end]-X[1:end-1].>0)

"""
    c=glue(a,b)

Glue together two vectors `a` and b resulting in a vector c. They last element
of `a` shall be equal (up to tol) to the first element of b.
The result fulfills `length(c)=length(a)+length(b)-1`
"""
function glue(_a::AbstractVector, _b::AbstractVector; tol=1.0e-10)
    a=collect_or_assign(_a)
    b=collect_or_assign(_b)

    is_monotone(a) || error("non-monotonous first argument of glue")
    is_monotone(b) || error("non-monotonous second argument of glue")
    Tv=promote_type(eltype(a),eltype(b))
    na=length(a)
    nb=length(b)

    d=b[1]-a[na-1]
    @assert(d>0)
    d=b[1]-a[na]
    @assert(d>-tol)
    @assert(d<tol)

    c=Vector{Tv}(undef,na+nb-1)
    ic=0
    for ia=1:na
        ic+=1
        c[ic]=a[ia]
    end
    for ib=2:nb
        ic+=1
        c[ic]=b[ib]
    end
    return c
end

ext: jl


text:
include("entrez_db.jl")
using .DB
using ..DBUtils
using SQLite

"""
save_efetch_sqlite(efetch_dict, db_config, verbose)

Save the results (dictionary) of an entrez fetch to a SQLite database.

###Arguments:

* `efetch_dict`: Response dictionary from efetch
* `db_config::Dict{Symbol, T}`: Configuration dictionary for initialitizing SQLite
database. Must contain symbols `:db_path` and `:overwrite`
* `verbose`: Boolean to turn on extra print statements

###Example

```julia
db_config =  Dict(:db_path=>"test_db.slqite", :overwrite=>true)
db = save_efetch_sqlite(efetch_dict, db_config)
```

"""
function save_efetch_sqlite{T}(efetch_dict, db_config::Dict{Symbol, T}, verbose=false)
    if haskey(efetch_dict, "PubmedArticle")
        db = init_pubmed_db_sqlite(db_config)
        return pubmed_save_efetch!(efetch_dict, db, verbose)
    else
        error("Unsupported efetch save. Responses must be searches to: PubMed")
        return nothing
    end
end

"""
save_efetch_mysql(efetch_dict, db_config, verbose)

Save the results (dictionary) of an entrez fetch to a MySQL database.

###Arguments:

* `efetch_dict`: Response dictionary from efetch
* `db_config::Dict{Symbol, T}`: Configuration dictionary for initialitizing SQLite
database. Must contain symbols `:host`, `:dbname`, `:username`. `pswd`,
and `:overwrite`
* `verbose`: Boolean to turn on extra print statements

###Example

```julia
db_config =  Dict(:host=>"localhost", :dbname=>"test", :username=>"root",
:pswd=>"", :overwrite=>true)
db = save_efetch_mysql(efetch_dict, db_config)
```

"""
function save_efetch_mysql{T}(efetch_dict, db_config::Dict{Symbol, T}, verbose=false)

    if haskey(efetch_dict, "PubmedArticle")
        db = init_pubmed_db_mysql(db_config)
        pubmed_save_efetch!(efetch_dict, db, verbose)
        return db
    else
        error("Unsupported efetch save. Responses must be searches to: PubMed")
        return nothing
    end
end

"""
pubmed_save_efetch(efetch_dict, db_path)

Save the results (dictionary) of an entrez-pubmed fetch to the input database.
"""
# function pubmed_save_efetch(efetch_dict, db, verbose=false)
#
#     if !haskey(efetch_dict, "PubmedArticle")
#         error("pubmed_save_efetch - error:PubmedArticleSet not found")
#         return
#     end
#
#     articles = efetch_dict["PubmedArticle"]
#
#     #articles should be an array
#     if !isa(articles, Array{Any, 1})
#         println("Error: Could not save to DB articles should be in an Array")
#         return
#     end
#
#     println("Saving " , length(articles) ,  " articles to database")
#
#     for article in articles
#
#         if !haskey(article,"MedlineCitation")
#             println("Error: Could not save to DB key:MedlineCitation not found")
#             return
#         end
#
#         pmid = nothing;
#         title = nothing;
#         pubYear = nothing;
#
#
#         # PMID is used as primary key - therefore it must be present
#         if haskey(article["MedlineCitation"][1],"PMID")
#             pmid = article["MedlineCitation"][1]["PMID"][1]["PMID"][1]
#         else
#             println("Error: Could not save to DB key:PMID not found - cannot be nothing")
#             return
#         end
#
#         # Retrieve basic article info
#         if haskey(article["MedlineCitation"][1],"Article")
#             if haskey(article["MedlineCitation"][1]["Article"][1], "ArticleTitle")
#                 title = article["MedlineCitation"][1]["Article"][1]["ArticleTitle"][1]
#             end
#             if haskey(article["MedlineCitation"][1]["Article"][1], "ArticleDate")
#                 if haskey(article["MedlineCitation"][1]["Article"][1]["ArticleDate"][1], "Year")
#                     pubYear = article["MedlineCitation"][1]["Article"][1]["ArticleDate"][1]["Year"][1]
#                 end
#             else  #series of attempts to pull a publication year from alternative xml elements
#                 try
#                     pubYear = article["MedlineCitation"][1]["Article"][1]["Journal"][1]["JournalIssue"][1]["PubDate"][1]["Year"][1]
#                 catch
#                     try
#                         pubYear = article["MedlineCitation"][1]["Article"][1]["Journal"][1]["JournalIssue"][1]["PubDate"][1]["MedlineDate"][1]
#                         pubYear = parse(Int64, pubYear[1:4])
#                     catch
#                         println("Warning: No date found")
#                     end
#                 end
#             end
#
#             # Save article data
#             insert_row!(db, "article", Dict(:pmid => pmid,
#             :title=>title,
#             :pubYear=>pubYear), verbose)
#
#             # insert all authors
#             forename = nothing
#             lastname = nothing
#             if haskey(article["MedlineCitation"][1]["Article"][1], "AuthorList")
#                 authors = article["MedlineCitation"][1]["Article"][1]["AuthorList"][1]["Author"]
#                 for author in authors
#
#                     if haskey(author, "ValidYN")
#                         if author["ValidYN"][1] == "N"
#                             println("Skipping Author Valid:N: ", author)
#                             continue
#                         end
#                     end
#
#                     if haskey(author, "ForeName")
#                         forename = author["ForeName"][1]
#                     else
#                         forname = "Unknown"
#                     end
#
#                     if haskey(author, "LastName")
#                         lastname = author["LastName"][1]
#                     else
#                         println("Skipping Author: ", author)
#                         continue
#                     end
#
#                     # Authors must be unique - if inserting fails see if already exists
#
#                     # Save author data
#                     author_id = -1
#                     try
#                         author_id = insert_row!(db, "author",
#                         Dict(:id => nothing,
#                         :forename => forename,
#                         :lastname => lastname), verbose)
#                     catch
#                         sel = db_select(db, ["id"], "author",
#                         Dict(:forename => forename, :lastname => lastname))
#                         if length(sel[1]) > 0
#                             author_id = get_value(sel[1][1])
#                             if verbose
#                                 println("Author already in db: ", lastname, ", ", forename, " - id: ", author_id)
#                             end
#                         else
#                             error("Can't save nor find Author: ", lastname, ", ", forename)
#                         end
#
#                     end
#
#                     if (author_id >= 0 )
#                         insert_row!(db, "author2article",
#                         Dict(:aid =>author_id, :pmid => pmid), verbose)
#                     else
#                         error("Invalid ID for Author: ", lastname, ", ", forename)
#                     end
#
#                 end
#             end
#
#             # Save related "keywords" of MESH Descriptors
#             if haskey(article["MedlineCitation"][1], "MeshHeadingList")
#                 if haskey(article["MedlineCitation"][1]["MeshHeadingList"][1], "MeshHeading")
#                     mesh_headings = article["MedlineCitation"][1]["MeshHeadingList"][1]["MeshHeading"]
#                     for heading in mesh_headings
#
#                         if !haskey(heading,"DescriptorName")
#                             println("Error: MeshHeading must have DescriptorName")
#                             return
#                         end
#
#                         #save descriptor
#                         descriptor_name = heading["DescriptorName"][1]["DescriptorName"][1]
#                         descriptor_name = normalize_string(descriptor_name, casefold=true)
#
#                         did = heading["DescriptorName"][1]["UI"][1]
#                         did_int = parse(Int64, did[2:end])  #remove preceding D
#
#                         #name of mesh descriptor must be unique
#                         try
#                             insert_row!(db, "mesh_descriptor",
#                             Dict(:id=>did_int, :name=>descriptor_name), verbose)
#                         catch
#                             try
#                                 sel = db_select(db, ["id"], "mesh_descriptor",
#                                 Dict(:name => descriptor_name))
#                                 if get_value(sel[1][1]) != did_int
#                                     error("Found matching descriptor but did is inconsistent")
#                                 end
#                                 if verbose
#                                     println("Descriptor already in db: ", descriptor_name, " - did: ", did)
#                                 end
#                             catch
#                                 error("Can't insert nor find duplicate")
#                             end
#                         end
#
#                         heading["DescriptorName"][1]["MajorTopicYN"][1] == "Y" ? dmjr = 1 : dmjr = 0
#
#                         #save the qualifiers
#                         if haskey(heading,"QualifierName")
#                             qualifiers = heading["QualifierName"]
#                             for qual in qualifiers
#                                 qualifier_name = qual["QualifierName"][1]
#                                 qualifier_name = normalize_string(qualifier_name, casefold=true)
#
#                                 qid = qual["UI"][1]
#                                 qid_int = parse(Int64, qid[2:end])  #remove preceding Q
#
#                                 try
#                                     insert_row!(db, "mesh_qualifier",
#                                     Dict(:id=>qid_int, :name=>qualifier_name), verbose )
#                                 catch
#                                     sel = db_select(db, ["id"], "mesh_qualifier",
#                                     Dict(:name => qualifier_name))
#                                     if get_value(sel[1][1]) != qid_int
#                                         error("Found matching qualifier but qid is inconsistent")
#                                     end
#                                     if verbose
#                                         println("Qualifier already in DB: ", qualifier_name, " - qid: ", qid)
#                                     end
#                                 end
#
#                                 qual["MajorTopicYN"][1] == "Y" ? qmjr = 1 : qmjr = 0
#
#                                 #save the heading related to this paper
#                                 insert_row!(db, "mesh_heading",
#                                 Dict(:id=>nothing, :pmid=> pmid, :did=>did_int,
#                                 :qid=>qid_int, :dmjr=>dmjr, :qmjr=>qmjr), verbose )
#
#                             end
#                         else
#                             #save the heading related to this paper
#                             insert_row!(db, "mesh_heading",
#                             Dict(:id=>nothing, :pmid=> pmid, :did=>did_int,
#                             :qid=>nothing, :dmjr=>dmjr, :qmjr=>nothing), verbose )
#                         end
#                     end
#                 end
#             end
#
#         end
#
#     end
#
#     return db
#
# end

function pubmed_save_efetch!(efetch_dict, db, verbose=false)

    #Decide type of article based on structrure of efetch
    articles = nothing
    if haskey(efetch_dict, "PubmedArticle")
        TypeArticle = PubMedArticle
        articles = efetch_dict["PubmedArticle"]
    else
        error("Save efetch is only supported for PubMed searches")
    end

    println("Saving " , length(articles) ,  " articles to database")

    for xml_article in articles
        article = TypeArticle(xml_article)

        db_insert!(db, article, verbose)

        #-------MeshHeadingList
        mesh_heading_list = MeshHeadingList(xml_article)
        db_insert!(db, article.pmid.value, mesh_heading_list, verbose)

    end
    db
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Libcap_Ng"
version = v"0.7.10"

# Collection of sources required to build libcap-ng
sources = [
    ArchiveSource("http://people.redhat.com/sgrubb/libcap-ng/libcap-ng-$(version).tar.gz",
                  "a84ca7b4e0444283ed269b7a29f5b6187f647c82e2b876636b49b9a744f0ffbf"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libcap-ng-*/
./configure --prefix=${prefix} --host=${target} \
    --enable-static=no \
    -with-python=no \
    CFLAGS="${CFLAGS} -pthread"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.  We are manually disabling
# many platforms that do not seem to work.
platforms = [p for p in supported_platforms() if Sys.islinux(p)]

# The products that we will ensure are always built
products = [
    LibraryProduct("libcap-ng", :libcap_ng),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Distributions, PDMats
using ArraysOfArrays
using ParallelProcessingTools
using DensityInterface

if Sys.isunix()
    @testset "external_density" begin
        test_src_dir = @__DIR__
        c_src_file = joinpath(test_src_dir, "external_mvnormal_density.cxx")

        dist = MvNormal([1.0 1.5; 1.5 4.0])

        mktempdir() do tmp_build_dir
            test_density_server_binary = joinpath(tmp_build_dir, "external_mvnormal_density")
            compile_cmd = `g++ -O2 $c_src_file -o $test_density_server_binary`
            @info "Running compile command $compile_cmd"
            run(compile_cmd)
            @info isfile(test_density_server_binary)
            likelihood = BAT.ExternalDensity(`$test_density_server_binary`, 0)

            @test begin
                x = [1.23, 2.34]
                llvalue = DensityInterface.logdensityof(likelihood, x)
                llvalue ≈ logpdf(dist, x)
            end

            @test begin
                n = 100
                params = nestedview(rand(2, n) .* 4 .- 2)
                ll_external = zeros(n)
                @sync for i in eachindex(params, ll_external)
                    @mt_async ll_external[i] = DensityInterface.logdensityof(likelihood, params[i])
                end
                ll_expected = logpdf.((dist,), params)
                ll_external ≈ ll_expected
            end

            close(likelihood)
        end
    end
end


text:
include("solver.jl")

"""
	Test if a solution is unique

	Prove existence of a unique solution by proving a solution is unique.

	Prove a solution is unique by showing that, for all combinations of
	all stars in a solution, if the star is blocked (set as not a possible
	part of the solution), and no solution can be found, the solution is unique.
"""
function unique_solution(puzzle, cells, k)
	d, n = size(puzzle)

	for i in 1:k
		stars = get_stars(cells)
		blocked = collect(combinations(1:k*n, i))
		for block in blocked
			# regions for testing solution
			regions = Dict{String,Array}(
				"rows" => Vector{Region}(undef, d),
				"cols" => Vector{Region}(undef, d),
				"groups" => Vector{Region}(undef, d),
			)

			# copy cells and fill regions based on cells
			c = deepcopy(cells)
			fill_regions!(regions, c)

			# unblock each cell
			for n in 1:d
				for m in 1:d
					c[m, n].is_blocked = false
				end
			end

			# get the idx of each star based on indexes from `block`
			star_idx = view(stars, block)

			# get all the stars based on `block` and remove from solution
			v = view(c, star_idx)
			for j in 1:size(v, 1)
				v[j].is_star = false
			end

			block!(c, regions, k)
			for j in 1:size(v, 1)
				v[j].is_blocked = true
			end
			update_regions!(regions, c)

			# try to solve with input cells and regions
			res = solve(puzzle, k)
			if !isnothing(res.cells)
				# another solution was found, the solution is not unique
				return res.cells
			end
		end
	end

	# no other solution found, the solution is unique
	return true
end


text:
using Documenter, SimSearchManifoldLearning

makedocs(;
    modules=[SimSearchManifoldLearning],
    authors="Eric S. Tellez",
    repo="https://github.com/sadit/SimSearchManifoldLearning.jl/blob/{commit}{path}#L{line}",
    sitename="SimSearchManifoldLearning.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", nothing) == "true",
        canonical="https://sadit.github.io/SimSearchManifoldLearning.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "UMAP" => "apiumap.md",
        "API ManifoldLearning" => "apiml.md"
    ],
)

deploydocs(;
    repo="github.com/sadit/SimSearchManifoldLearning.jl",
    devbranch=nothing,
    branch = "gh-pages",
    versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"]
)


text:
export svdK,rotate,logistic,Lambda,Seed, dropCol

#Spectral decomposition
"""

svdK(K::Union{Array{Float64,3},Matrix{Float64}};LOCO::Bool=true)

Returns a 3d-array (or a matrix) of orthogonal vectors and the corresponding matrix (or vector) of singluar values if `LOCO` is true (false).

# Arguments

- `K` : a 3d-array of a symmetric positive definite kinship, `size(K) = (n,n, # of Chromosomes)` if `LOCO = true` (default).  Or a n x n matrix if false.

## Keyword Arguments

- `LOCO` : Boolean.  Default is true.

# Output

- `T` : a 3d-array (or matrix) of orthogonal vectors
- `S` : a matrix (or vector) of singluar values in descending order

"""
function svdK(K::Union{Array{Float64,3},Matrix{Float64}};LOCO::Bool=true,δ::Float64=0.001)

    if (LOCO)
        nchr = size(K,3);  n = size(K,1);
        S = zeros(n,nchr);
        T = zeros(n,n,nchr);

         for j =1:nchr
            # may consider to svd from the genotype data
            if(!isposdef(K[:,:,j])) # check pdf
                 K[:,:,j]=K[:,:,j]+ (abs(eigmin(K[:,:,j]))+δ)*I
            end

        #cholesky & svd for numerical stability
        #   C=cholesky(K[:,:,j])
          F=svd(K[:,:,j])
          T[:,:,j], S[:,j] = F.Vt, F.S
         end

         return  T, S
     else
           if(!isposdef(K)) # check pdf
                 K=K+ (abs(eigmin(K))+δ)*I
            end
        #  C=cholesky(K)
         F=svd(K)
         return F.Vt, F.S

     end

end

"""

    rotate(y::Vector{Float64},X::Matrix{Float64},X₀::Matrix{Float64},T::Matrix{Float64})
    rotateX(X::Matrix{Float64},T::Matrix{Float64})
    rotateY(y::Vector{Float64},T::Matrix{Float64})

Returns transformed data rotated by eigenvectors.  `rotateX` transforms a matrix, `rotateY` does a vector, and `rotate` does both.

# Arguments

- `y` : a n x 1 vector of binary trait
- `X` : a n x p matrix of genetic markers selected from QTL analysis (per Chromosome for LOCO)
- `X₀` : a n x c matrix of covariates including intercepts
- `T` : a n x n matrix of orthogonal matrix (eigen vectors from kinship, `K`).  See also [`eigenK`](@ref).

# Output

- `yt` : a tranformed `y`
- `Xt` : a transformed `X`
- `Xt₀` : a transformed `X₀`

"""
function rotate(y::Vector{Float64},X::Matrix{Float64},X₀::Matrix{Float64},T::Matrix{Float64})

    Xt= rotateX(X,T)
    Xt₀=rotateX(X₀,T)
    yt= rotateY(y,T)

    return Xt, Xt₀, yt

end

function rotateX(X::Matrix{Float64},T::Matrix{Float64})

    #U'X, U'X₀,U'(y-1/2)

    Xt=BLAS.gemm('N','N',T,X)

    return Xt

end

function rotateY(y::Vector{Float64},T::Matrix{Float64})

    n=length(y)

    #U'X, U'X₀, U'(y-1/2)
    yt= BLAS.gemv('N',T,(y.-0.5))

    return yt

end

#compute X'y : 't': 'T' for transpose, 'N' for no transpose
function getXy(t::Char,Xt::Matrix{Float64},yt::Vector{Float64})

    y1= BLAS.gemv(t,Xt,yt)

    return y1
end

#matrix multiplication AB: 'tA', 'tB': 'T' for transpose, 'N' for no transpose
function getXX(tA::Char,A::Matrix{Float64},tB::Char,B::Matrix)

     X = BLAS.gemm(tA,tB,A,B)

    return X

end

# compute symmetric square matrix, t ='T' for X'X, t='N' for XX'
function symXX(t::Char,X::Array{Float64,2})

   XtX = Symmetric(BLAS.syrk('U',t,1.0,X))

    return XtX

end

function dropCol(B::Matrix{Float64},l::Int64)

    return B[:,1:end.!= l]
  end

"""

    logistic(ξ::Float64)

A logistic function computing the probablity of an event.

"""
function logistic(ξ::Float64)

    return 1.0/(1.0+ exp(-ξ))

end

# 2λ(ξ)
function Lambda(ξ::Float64)
        if(ξ==0.0) #adjust 0/0 case by L'Hopital
            return 0.25
        else
            return 0.5*tanh(ξ/2)/ξ
        end
end

"""

    GenoInfo(snp::Array{String,1},chr::Array{Any,1},pos::Array{Float64,1})

A type struct includes information on genotype data

# Arguements

- `snp` : a vector of strings for genetic marker (or SNP) names
- `chr` : a vector of Chromosome information corresponding to `snp`
- `pos` : a vector of `snp` position in `cM`

"""
struct GenoInfo
    snp::Array{String,1}
    chr::Array{Any,1}
    pos::Array{Float64,1} #positon
end

"""

    Seed(M::Int64)

Reseed random number generators to assign workers (or processes) different number of seeds (increased by 1) for reproducible distributed computing.

# Arguments

- `M` : an integer.  Set the seed number.

# Examples

```julia

julia> using Distributed, Random
julia> addprocs(10)    # generate 10 workers (or processes)
julia> @everywhere using SuSiEGLMM  # access the pkg to all workers
julia> Seed(20)

```

"""
function Seed(M::Int64)

        np=nprocs(); pid=procs()

        seeds = collect(M:1:M+np-1)

        for i =1:np
            remotecall_fetch(()->Random.seed!(seeds[i]),pid[i])
        end

end


text:
const DimArrayOrStack = Union{AbstractDimArray,AbstractDimStack}

"""
    set(x, val)
    set(x, args::Pairs...) => x with updated field/s
    set(x, args...; kw...) => x with updated field/s
    set(x, args::Tuple{Vararg{<:Dimension}}; kw...) => x with updated field/s

    set(dim::Dimension, index::AbstractArray) => Dimension
    set(dim::Dimension, lookup::LookupArray) => Dimension
    set(dim::Dimension, lookupcomponent::LookupArrayTrait) => Dimension
    set(dim::Dimension, metadata::AbstractMetadata) => Dimension

Set the properties of an object, its internal data or the traits of its dimensions
and lookup index.

As DimensionalData is so strongly typed you do not need to specify what field
of a [`LookupArray`](@ref) to `set` - there is no ambiguity.

To set fields of a `LookupArray` you need to specify the dimension. This can be done
using `X => val` pairs, `X = val` keyword arguments, or `X(val)` wrapped arguments.

When a `Dimension` or `LookupArray` is passed to `set` to replace the
existing ones, fields that are not set will keep their original values.

## Notes:

Changing a lookup index range/vector will also update the step size and order where applicable.

Setting the [`Order`](@ref) like `ForwardOrdered` will *not* reverse the array or
dimension to match. Use `reverse` and [`reorder`](@ref) to do this.

## Examples

```jldoctest set
julia> using DimensionalData; const DD = DimensionalData
DimensionalData

julia> da = DimArray(zeros(3, 4), (custom=10.0:010.0:30.0, Z=-20:010.0:10.0));

julia> set(da, ones(3, 4))
3×4 DimArray{Float64,2} with dimensions:
  Dim{:custom} Sampled 10.0:10.0:30.0 ForwardOrdered Regular Points,
  Z Sampled -20.0:10.0:10.0 ForwardOrdered Regular Points
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
```

Change the `Dimension` wrapper type:

```jldoctest set
julia> set(da, :Z => Ti, :custom => Z)
3×4 DimArray{Float64,2} with dimensions:
  Z Sampled 10.0:10.0:30.0 ForwardOrdered Regular Points,
  Ti Sampled -20.0:10.0:10.0 ForwardOrdered Regular Points
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
```

Change the lookup `Vector`:

```jldoctest set
julia> set(da, Z => [:a, :b, :c, :d], :custom => [4, 5, 6])
3×4 DimArray{Float64,2} with dimensions:
  Dim{:custom} Sampled Int64[4, 5, 6] ForwardOrdered Regular Points,
  Z Sampled Symbol[a, b, c, d] ForwardOrdered Regular Points
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
```

Change the `LookupArray` type:

```jldoctest set
julia> set(da, Z=DD.NoLookup(), custom=DD.Sampled())
3×4 DimArray{Float64,2} with dimensions:
  Dim{:custom} Sampled 10.0:10.0:30.0 ForwardOrdered Regular Points,
  Z
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
```

Change the `Sampling` trait:

```jldoctest set
julia> set(da, :custom => DD.Irregular(10, 12), Z => DD.Regular(9.9))
3×4 DimArray{Float64,2} with dimensions:
  Dim{:custom} Sampled 10.0:10.0:30.0 ForwardOrdered Irregular Points,
  Z Sampled -20.0:10.0:10.0 ForwardOrdered Regular Points
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
```
"""
function set end
set(A::DimArrayOrStack, x::T) where {T<:Union{LookupArray,LookupArrayTrait}} = _onlydimerror(x)
set(x::DimArrayOrStack, ::Type{T}) where T = set(x, T())

set(A::AbstractDimStack, x::LookupArray) = LookupArrays._cantseterror(A, x)
set(A::AbstractDimArray, x::LookupArray) = LookupArrays._cantseterror(A, x)
set(A, x) = LookupArrays._cantseterror(A, x)
set(A::DimArrayOrStack, args::Union{Dimension,DimTuple,Pair}...; kw...) =
    rebuild(A, data(A), _set(dims(A), args...; kw...))
set(A::AbstractDimArray, newdata::AbstractArray) = begin
    axes(A) == axes(newdata) || _axiserr(A, newdata)
    rebuild(A; data=newdata)
end
set(s::AbstractDimStack, newdata::NamedTuple) = begin
    dat = data(s)
    keys(dat) === keys(newdata) || _keyerr(keys(dat), keys(newdata))
    map(dat, newdata) do d, nd
        axes(d) == axes(nd) || _axiserr(d, nd)
    end
    rebuild(s; data=newdata)
end

@noinline _onlydimerror(x) = throw(ArgumentError("Can only set $(typeof(x)) for a dimension. Specify which dimension you mean with `X => property`"))
@noinline _axiserr(a, b) = throw(ArgumentError("passed in axes $(axes(b)) do not match the currect axes $(axes(a))"))
@noinline _keyerr(ka, kb) = throw(ArgumentError("keys $ka and $kb do not match"))


text:
"""
	gmt2kml(cmd0::String="", arg1=nothing, kwargs...)

Convert GMT data tables to KML files for Google Earth

Full option list at [`gmt2kml`]($(GMTdoc)gmt2kml.html)

Parameters
----------

- **A** | **altitude_mode** :: [Type => Str]       ``Arg = a|g|s[alt|xscale]``

    Select one of three altitude modes recognized by Google Earth that determines the altitude (in m)
    of the feature: ``a`` absolute altitude, ``g`` altitude relative to sea surface or ground,
    ``s`` altitude relative to seafloor or ground.
    ($(GMTdoc)gmt2kml.html#a)
- $(GMT.opt_C)
- **D** | **descript** :: [Type => Str]   ``Arg = descriptfile``

    File with HTML snippets that will be included as part of the main description content for the KML file.
    ($(GMTdoc)gmt2kml.html#d)
- **E** | **extrude** :: [Type => Str | []]  ``Arg = [altitude]``

    Extrude feature down to ground level.
    ($(GMTdoc)gmt2kml.html#e)
- **F** | **feature_type** :: [Type => Str]  ``Arg = e|s|t|l|p|w``

    Sets the feature type. Choose from points (event, symbol, or timespan), line, polygon, or wiggle.
    ($(GMTdoc)gmt2kml.html#f)
- **G** | **fill** :: [Type => Str]  ``Arg = f|nfill``

    Sets color fill (G=:f) or label font color (G=:n).
    ($(GMTdoc)gmt2kml.html#g)
- **I** | **icon** :: [Type => Str]      ``Arg = icon``

    Specify the URL to an alternative icon that should be used for the symbol
    [Default is a Google Earth circle].
    ($(GMTdoc)gmt2kml.html#i)
- **K** | **not_finished** :: [Type => Bool]

    Allow more KML code to be appended to the output later [finalize the KML file].
    ($(GMTdoc)gmt2kml.html#k)
- **L** | **extended_data** :: [Type => Str]      ``Arg = name1,name2,…``

    Extended data given. Append one or more column names separated by commas.
    ($(GMTdoc)gmt2kml.html#l)
- **N** | **feature_name** :: [Type => Str | Number]      ``Arg = [t|col |name_template|name]``

    By default, if segment headers contain a -L”label string” then we use that for the name of the KML feature.
    ($(GMTdoc)gmt2kml.html#n)
- **O** | **overlay** :: [Type => Bool]

    Append KML code to an existing KML file [initialize a new KML file].
    ($(GMTdoc)gmt2kml.html#o)
- **Qa** | **wiggles** :: [Type => Str]      ``Arg =  azimuth``

    Option in support of wiggle plots (requires F=:w).
    ($(GMTdoc)gmt2kml.html#q)
- **Qs** | **wiggle_scale** :: [Type => Str | Number]      ``Arg =  scale[unit]``

    Required setting for wiggle plots (i.e., it requires F=:w). Sets a wiggle scale
    in z-data units per the user’s units
    ($(GMTdoc)gmt2kml.html#q)
- $(GMT.opt_R)
- **S** | **scale** :: [Type => Str]      ``Arg =  c|nscale``

    Scale icons or labels. Here, S=:c sets a scale for the symbol icon, whereas S=:n sets
    a scale for the name labels
    ($(GMTdoc)gmt2kml.html#s)
- **T** | **title** :: [Type => Str]    ``Arg = title[/foldername]``

    Sets the document title [default is unset]. Optionally, append /FolderName;
    ($(GMTdoc)gmt2kml.html#t)
- **W** | **pen** :: [Type => Str | []]      ``Arg =  [pen][attr]``

    Set pen attributes for lines, wiggles or polygon outlines.
    ($(GMTdoc)gmt2kml.html#w)
- **Z** | **attrib** :: [Type => Str]      ``Arg =  args``

    Set one or more attributes of the Document and Region tags.
    ($(GMTdoc)gmt2kml.html#z)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_swap_xy)
"""
function gmt2kml(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("gmt2kml", cmd0, arg1)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode

	cmd, = parse_common_opts(d, "", [:R :V_params :bi :di :e :f :h :i :yx])
	cmd  = parse_these_opts(cmd, d, [[:A :altitude_mode], [:D :descript], [:E :extrude], [:F :feature_type],
	                                 [:I :icon], [:K :not_finished], [:L :extended_data], [:N :feature_name], [:O :overlay], [:Qa :wiggles], [:Qs :wiggle_scale], [:S :scale], [:T :title], [:Z :attrib]])

	cmd = add_opt(d, cmd, "G", [:G :fill])
	cmd *= add_opt_pen(d, [:W :pen], "W", true)     # TRUE to also seek (lw,lc,ls)

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, arg1)
	N_used = got_fname == 0 ? 1 : 0			# To know whether a cpt will go to arg1 or arg2
	cmd, arg1, arg2, = add_opt_cpt(d, cmd, CPTaliases, 'C', N_used, arg1)
	common_grd(d, "gmt2kml " * cmd, arg1, arg2)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
gmt2kml(arg1, cmd0::String=""; kw...) = gmt2kml(cmd0, arg1; kw...)


text:
using SearchablePDFs
using Documenter

makedocs(; modules=[SearchablePDFs], authors="Eric P. Hanson",
         repo="https://github.com/ericphanson/SearchablePDFs.jl/blob/{commit}{path}#L{line}",
         sitename="SearchablePDFs.jl",
         format=Documenter.HTML(; prettyurls=get(ENV, "CI", "false") == "true",
                                assets=String[]), pages=["Home" => "index.md"])


text:
function process_icr!(;collection, openable_cards, account, rarity, parameters, logging_vector)
    if isempty(openable_cards.by_rarity[rarity])
        @goto vault_increment
    end
    card = rand(openable_cards.by_rarity[rarity])
    if collection[card] < 4
        collection[card] += 1
        push!(logging_vector, card)
        if collection[card] >= 4 && parameters.prevent_duplicates
            remove_fully_collected!(card, rarity; openable_cards = openable_cards)
        end
    else
        @label vault_increment
        account.vault_pct += parameters.duplicate_vpct[rarity]
    end
    collection
end

max_stars_repo_name: asieev/ArenaSim

max_issues_repo_name: asieev/ArenaSim

max_forks_repo_name: asieev/ArenaSim


text:
using Bridge
using StaticArrays
import Bridge: b, σ, B, β, a, constdiff
const ℝ = SVector{N, T} where {N, T}

"""
    JRNeuralDiffusion3n <: ContinuousTimeProcess{ℝ{6, T}}
structure defining the Jansen and Rit Neural Mass Model described in
https://mathematical-neuroscience.springeropen.com/articles/10.1186/s13408-017-0046-4 and
https://arxiv.org/abs/1903.01138
"""
struct JRNeuralDiffusion3n{T} <: ContinuousTimeProcess{ℝ{6, T}}
    A::T
    a::T
    B::T
    b::T
    C::T
    νmax::T
    v0::T
    r::T
    μx::T
    μy::T
    μz::T
    σx::T
    σy::T
    σz::T
    # constructor given assumption statistical paper
    function JRNeuralDiffusion3n(A::T, a::T, B::T, b::T, C::T,
            νmax::T, v0::T ,r::T, μx::T, μy::T, μz::T, σx::T, σy::T, σz::T) where T
        new{T}(A, a, B, b, C, νmax, v0, r, μx, μy, μz, σx, σy, σz)
    end
end
#C1 = C, C2 = 0.8C, c3 = 0.25C, c4 =  0.25C,

# in the statistical paper they set μ's to be constant and not function of time.
function μx(t, P::JRNeuralDiffusion3n)
    P.μx
end

function μy(t, P::JRNeuralDiffusion3n)
    P.μy
end

function μz(t, P::JRNeuralDiffusion3n)
    P.μz
end

"""
    sigm(x, P::JRNeuralDiffusion3n)
definition of sigmoid function
"""
function sigm(x, P::JRNeuralDiffusion3n{T}) where T
    P.νmax / (1 + exp(P.r*(P.v0 - x)))
end

function b(t, x, P::JRNeuralDiffusion3n{T}) where T
    ℝ{6, T}(x[4], x[5], x[6],
    P.A*P.a*(μx(t, P) + sigm(x[2] - x[3], P)) - 2P.a*x[4] - P.a*P.a*x[1],
    P.A*P.a*(μy(t, P) + 0.8P.C*sigm(P.C*x[1], P)) - 2P.a*x[5] - P.a*P.a*x[2],
    P.B*P.b*(μz(t, P) + 0.25P.C*sigm(0.25P.C*x[1], P)) - 2P.b*x[6] - P.b*P.b*x[3])
end

#6x3 matrix
function σ(t, x, P::JRNeuralDiffusion3n{T}) where T
    @SMatrix    [0.0  0.0  0.0;
                0.0  0.0  0.0;
                0.0 0.0  0.0;
                P.σx  0.0  0.0;
                0.0  P.σy 0.0;
                0.0  0.0  P.σz]
end

constdiff(::JRNeuralDiffusion3n) = true
clone(::JRNeuralDiffusion3n, θ) = JRNeuralDiffusion3n(θ...)
#Static vector
params(P::JRNeuralDiffusion3n) = [P.A, P.a, P.B, P.b, P.C, P.νmax,
    P.v0, P.r, P.μx, P.μy, P.μz, P.σx, P.σy, P.σz]
param_names(::JRNeuralDiffusion3n) = (:A, :a, :B, :b, :C, :νmax,
    :v0, :r, :μx, :μy, :μz, :σx, :σy, :σz)

#### TODO AGAIN ####
phi(::Val{0}, t, x, P::JRNeuralDiffusion3n) = (P.A*P.a*0.8P.C*sigm(P.C*x[1], P) - 2P.a*x[5] - P.a*P.a*x[2],)
phi(::Val{10}, t, x, P::JRNeuralDiffusion3n) = (P.A*P.a,)
phi(::Val{1}, t, x, P::JRNeuralDiffusion3n) = (0,)
phi(::Val{2}, t, x, P::JRNeuralDiffusion3n) = (0,)
phi(::Val{3}, t, x, P::JRNeuralDiffusion3n) = (0,)
phi(::Val{4}, t, x, P::JRNeuralDiffusion3n) = (0,)
phi(::Val{5}, t, x, P::JRNeuralDiffusion3n) = (0,)
phi(::Val{6}, t, x, P::JRNeuralDiffusion3n) = (0,)
phi(::Val{7}, t, x, P::JRNeuralDiffusion3n) = (0,)
phi(::Val{8}, t, x, P::JRNeuralDiffusion3n) = (0,)
phi(::Val{9}, t, x, P::JRNeuralDiffusion3n) = (0,)
phi(::Val{11}, t, x, P::JRNeuralDiffusion3n) = (0,)
phi(::Val{12}, t, x, P::JRNeuralDiffusion3n) = (0,)

"""
    JRNeuralDiffusion3nAux{T, S1, S2} <: ContinuousTimeProcess{ℝ{6, T}}
structure for a simple auxiliary process defined as integrated Wiener process
"""
struct JRNeuralDiffusion3nAux{R, S1, S2} <: ContinuousTimeProcess{ℝ{6, R}}
    σx::R
    σy::R
    σz::R
    t::Float64
    u::S1
    T::Float64
    v::S2
    # generator given assumptions paper
    function JRNeuralDiffusion3nAux(σx::R, σy::R, σz::R,  t::Float64, u::S1,
                        T::Float64, v::S2) where {R, S1, S2}
        new{R, S1, S2}(σx, σy, σz, t, u, T, v)
    end
end

function B(t, P::JRNeuralDiffusion3nAux{S1, S2}) where {S1, S2}
    @SMatrix [0.0  0.0  0.0  1.0  0.0  0.0;
              0.0  0.0  0.0  0.0  1.0  0.0;
              0.0  0.0  0.0  0.0  0.0  1.0;
              0.0  0.0  0.0  0.0  0.0  0.0;
              0.0  0.0  0.0  0.0  0.0  0.0;
              0.0  0.0  0.0  0.0  0.0  0.0]
end

function β(t, P::JRNeuralDiffusion3nAux{S1, S2}) where {S1, S2}
    ℝ{6}(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
end

#6x3 matrix
function σ(t, P::JRNeuralDiffusion3nAux{T}) where T
    @SMatrix    [0.0  0.0  0.0;
                0.0  0.0  0.0;
                0.0 0.0  0.0;
                P.σx  0.0  0.0;
                0.0  P.σy 0.0;
                0.0  0.0  P.σz]
end

b(t, x, P::JRNeuralDiffusion3nAux) = B(t,P) * x + β(t,P)
a(t, P::JRNeuralDiffusion3nAux) = σ(t,P) * σ(t, P)'

constdiff(::JRNeuralDiffusion3nAux) = true
clone(P::JRNeuralDiffusion3nAux, θ) = JRNeuralDiffusion3nAux(θ[12], θ[13], θ[14], P.t, P.u, P.T, P.v)
clone(P::JRNeuralDiffusion3nAux, θ, v) = JRNeuralDiffusion3nAux(θ[12], θ[13], θ[14], P.t, zero(v), P.T, v)
params(P::JRNeuralDiffusion3nAux) = [P.σx, P.σy, P.σz]
param_names(P::JRNeuralDiffusion3nAux) = (:σx, :σy, :σz )
depends_on_params(::JRNeuralDiffusion3nAux) = (12,13,14)


text:
"""
    arbrs(RR::Vector{<:Real}, SBP::Vector{<:Real}; p::Union{Int, UnitRange{Int}} = 12:18, nfreq::Int = 256, sf::Real = 1, LF::Tuple{<:Real, <:Real} = (0.04, 0.15), HF::Tuple{<:Real, <:Real} = (0.15, 0.4))

Calculate the BRS indices αLF, αHF using AR spectral decomposition.

# Keywords

* p: order of the AR process, if a UnitRange is given, the optimal order will be chosen through AIC, defaults to 12:18.
* nfreq: resolution of the analysed spectra, defaults to 256
* sf: sampling frequency in Hz, defaults to 1
* LF: Tuple that defined the LF band, defaults to (0.04, 0.15)
* HF: Tuple that defined the HF band, defaults to (0.15, 0.4)
* verbose: print some information oder selection, defaults to false
"""
function arbrs(RR::Vector{<:Real}, SBP::Vector{<:Real}; p::Union{Int, UnitRange{Int}} = 12:18, nfreq::Int = 256, sf::Real = 1, LF::Tuple{<:Real, <:Real} = (0.04, 0.15), HF::Tuple{<:Real, <:Real} = (0.15, 0.4), verbose::Bool = false)
    nfreq >= 16 || throw(DomainError("The resolution should be reasonable high."))
    sf > 0 || throw(DomainError("The sampling frequency needs to be positive."))
    all([LF..., HF...] .>= 0) || throw(DomainError("Frequencies in LF/HF band need to be greater 0"))

    SRR, centerFrequenciesRR, _ = arDecomposition(RR, p, nfreq = nfreq, sf = sf, verbose = verbose)
    SSBP, centerFrequenciesSBP, f = arDecomposition(SBP, p, nfreq = nfreq, sf = sf, verbose = verbose)

    LFRR = getSpectralComponent(SRR, centerFrequenciesRR, LF)
    HFRR = getSpectralComponent(SRR, centerFrequenciesRR, HF)
    fullRR = getSpectralComponent(SRR, centerFrequenciesRR)

    LFSBP = getSpectralComponent(SSBP, centerFrequenciesSBP, LF)
    HFSBP = getSpectralComponent(SSBP, centerFrequenciesSBP, HF)
    fullSBP = getSpectralComponent(SSBP, centerFrequenciesSBP)

    αLF, αHF = NaN, NaN
    try # something became negative here
        αLF = sqrt(sum(LFRR) / sum(LFSBP))
        αHF = sqrt(sum(HFRR) / sum(HFSBP))
    catch
        @warn "Some spectral components were negative."
    end

    return arBRS(αLF = αLF, αHF = αHF, LFRR = LFRR, LFSBP = LFSBP, HFRR = HFRR, HFSBP = HFSBP, f = collect(f), fullRR = fullRR, fullSBP = fullSBP)
end

"""
Struct that stores all information regarding the arBRS estimation. The final result is stored in `αLF` & `αHF`. It can be plotted for visual inspection.
"""
@with_kw mutable struct arBRS
    αLF::Real = 0
    αHF::Real = 0
    LFRR::Vector{<:Real} = Float64[]
    LFSBP::Vector{<:Real} = Float64[]
    HFRR::Vector{<:Real} = Float64[]
    HFSBP::Vector{<:Real} = Float64[]
    f::Vector{<:Real} = Float64[]
    fullRR::Vector{<:Real} = Float64[]
    fullSBP::Vector{<:Real} = Float64[]
end

@recipe function f(res::arBRS)
    layout := (2,1)
    freq = res.f

    colHF = :red
    colLF = :green
    cap(spec, lim) = replace(x -> x >= lim ? lim : x , spec)

    @series begin
        subplot := 1
        label := ""
        seriescolor := :black
        linewidth := 2
        fill := (0, :grey)
        seriesalpha := 0.6

        freq, cap(res.fullSBP, 2maximum([res.LFSBP; res.HFSBP]))
    end

    @series begin
        subplot := 1
        label := "HF"
        seriescolor := :black
        fill := (0, colHF)
        seriesalpha := 0.6

        freq, res.HFSBP
    end

    @series begin
        subplot := 1
        yguide := "mmHg²"
        label := "LF"
        seriescolor := :black
        ylims := (0, 2maximum([res.LFSBP; res.HFSBP]))
        fill := (0, colLF)
        seriesalpha := 0.6

        freq, res.LFSBP
    end

    @series begin
        subplot := 2
        yguide := "ms²"
        label := ""
        seriescolor := :black
        linewidth := 2
        fill := (0, :grey)
        seriesalpha := 0.6

        freq, cap(res.fullRR, 2maximum([res.LFRR; res.HFRR]))
    end

    @series begin
        subplot := 2
        yguide := "ms²"
        label := ""
        seriescolor := :black
        fill := (0, colHF)
        seriesalpha := 0.6

        freq, res.HFRR
    end

    @series begin
        subplot := 2
        seriestype := :scatter
        marker := :circ
        label := "αHF: $(round(res.αHF, digits = 2)) ms/mmHg"
        seriescolor := :black

        [0,0], [-10,-10]
    end

    @series begin
        subplot := 2
        seriestype := :scatter
        marker := :circ
        label := "αLF: $(round(res.αLF, digits = 2)) ms/mmHg"
        seriescolor := :black

        [0,0], [-10,-10]
    end

    @series begin
        subplot := 2
        yguide := "ms²"
        label := ""
        seriescolor := :black
        ylims := (0, 2maximum([res.LFRR; res.HFRR]))
        fill := (0, colLF)
        seriesalpha := 0.6
        xguide := "f [Hz]"

        freq, res.LFRR
    end
end


text:
using Match

input = joinpath(@__DIR__, "input")
instructions = readlines(input)

instructions = map(instructions) do ins
    ins = split(ins)
    (first(ins), parse(Int, last(ins)))
end

function navigate(instructions, depth, pos)
    for i in instructions
        (op, val) = i
        @match op begin
            "forward" => begin pos += val end
            "up" => begin depth -= val end
            "down" => begin depth += val end
        end
    end
    return depth * pos
end

function navigate(instructions, depth, pos, aim)
    for i in instructions
        (op, val) = i
        @match op begin
            "forward" => begin
                pos += val
                depth += (val * aim)
            end
            "up" => begin aim -= val end
            "down" => begin aim += val end
        end
    end
    return depth * pos
end

p1 = navigate(instructions, 0, 0)
p2 = navigate(instructions, 0, 0, 0)

println("-----------------------------------------------------------------------")
println("dive! -- part one :: $p1")
println("dive! -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 1714680)
@assert(p2 == 1963088820)


text:
"""
    Heatmap(xedges, yedges; left = true, closed = true)

Create a two dimensional histogram with the bin partition created by `xedges` and `yedges`.
When fitting a new observation, the first value will be associated with X, the second with Y.

- If `left`, the bins will be left-closed.
- If `closed`, the bins on the ends will be closed.  See [Hist](@ref).

# Example

    o = fit!(HeatMap(-5:.1:5, -5:.1:5), eachrow(randn(10^5, 2)))

    using Plots
    plot(o)
"""
mutable struct HeatMap{EX, EY} <: OnlineStat{XY}
    xedges::EX
    yedges::EY
    counts::Matrix{Int}
    out::Int
    left::Bool
    closed::Bool

    function HeatMap(x::T, y::S; left::Bool=true, closed::Bool=true) where {T<:AbstractVector,S<:AbstractVector}
        new{T,S}(x,y, zeros(Int, length(x)-1, length(y)-1), 0, left, closed)
    end
end
nobs(o::HeatMap) = sum(o.counts) + o.out
value(o::HeatMap) = (x=o.xedges, y=o.yedges, z=o.counts)

function _fit!(o::HeatMap, xy)
    i = binindex(o.xedges, xy[1], o.left, o.closed)
    if 1 ≤ i < length(o.xedges)
        j = binindex(o.yedges, xy[2], o.left, o.closed)
        if 1 ≤ j < length(o.yedges)
            o.counts[i, j] += 1
        else
            o.out += 1
        end
    else
        o.out += 1
    end
end

function _merge!(a::HeatMap, b::HeatMap)
    if (a.xedges == b.xedges) && (a.yedges == b.yedges) && (a.left == b.left) && (a.closed == b.closed)
        a.out += b.out
        a.counts .+= b.counts
    else
        @warn "HeatMaps differ in edges or direction of bins being closed.  No merging occurred."
    end
end

@recipe function f(o::HeatMap)
    seriestype --> :heatmap
    z = Float64.(o.counts)
    z[z .== 0] .= NaN
    midpoints(o.xedges), midpoints(o.yedges), z
end


text:
#=
Simple example of use of montecarlo.jl
From the system prompt, execute
mpirun -np X julia pi_mpi.jl
setting X to the number of ranks you'd
like to use, subject to  X-1 being an even divisor
of 1e6, e.g., set X=5.
=#

using MPI, Econometrics, LinearAlgebra, Statistics

function pi_wrapper()
    4.0 * (norm(rand(2)) .< 1.0)
end

# this function reports intermediate results during MC runs
function pi_monitor(sofar, results)
    # Examine results every 12.5*10^5 draws
    if mod(sofar, 10^6) == 0
        m = mean(results[1:sofar,:],dims=1)
        println("reps: $sofar, pihat: $m")
    end
end

# do the monte carlo: 10^6 reps of single draws
function main()
    MPI.Init()
    n_evals = 10^7 # desired number of MC reps
    n_returns = 1
    batchsize = 10^5
    montecarlo(pi_wrapper, pi_monitor,
               MPI.COMM_WORLD, n_evals, n_returns, batchsize)
    MPI.Finalize()
end

main()

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
using InvertedIndices
using Test
using OffsetArrays

@testset "0-d" begin
    A = fill(1)
    @test A[Not(fill(A.==1))] == []
    @test A[Not(CartesianIndex())] == []
    @test A[Not(fill(A.==2))] == [1]
    A[Not(fill(A.==2))] = fill(0)
    @test A[] == 0
    A[Not(fill(A.==2))] .= 1
    @test A[] == 1
    @test A[Not([1,1])] == []
    @test A[Not(1,1)] == []
end

@testset "1-d readonly" for A in (-10:13, reshape(-10:13,2,:), reshape(-10:13,3,2,:))
    @test A[Not(1)] == A[Not(1:1)] == A[Not(A.==-10)] == collect(-9:13)
    @test @views A[Not(1)] == A[Not(1:1)] == A[Not(A.==-10)] == collect(-9:13)
    @test A[Not(end)] == A[Not(end:end)] == A[Not(A.==13)] == collect(-10:12)
    @test @views A[Not(end)] == A[Not(end:end)] == A[Not(A.==13)] == collect(-10:12)
    @test A[Not(iseven.(A))] == A[isodd.(A)] == collect(-9:2:13)
    @test A[Not([])] == A[collect(1:end)] == collect(-10:13)
    @test A[Not(1:end)] == A[Not(:)] == A[[]] == []
    @test A[Not([1,1,1,2,2])] == A[3:end]
    @test A[Not([end,end,end,end-1,end-1])] == A[1:end-2]
    @test A[Not(3,2,1)] == A[Not(3,2,3,1,2,1)] == A[4:end]

    @test_throws BoundsError A[Not(0)]
    @test_throws BoundsError A[Not(end+1)]
    @test_throws BoundsError A[Not(0:end)]
    @test_throws BoundsError A[Not(1:end+1)]
end

@testset "1-d offset readonly" for A in (OffsetArray(-10:13, -3), OffsetArray(reshape(-10:13,2,:), -5, 32), OffsetArray(reshape(-10:13,3,2,:), 4,-2,20))
    f = first(LinearIndices(A))
    l = last(LinearIndices(A))
    @test A[Not(f)] == A[Not(f:f)] == A[Not(A.==-10)] == collect(-9:13)
    @test @views A[Not(f)] == A[Not(f:f)] == A[Not(A.==-10)] == collect(-9:13)
    @test A[Not(l)] == A[Not(l:l)] == A[Not(A.==13)] == collect(-10:12)
    @test @views A[Not(l)] == A[Not(l:l)] == A[Not(A.==13)] == collect(-10:12)
    @test A[Not(iseven.(A))] == A[isodd.(A)]
    @test A[Not([])] == A[collect(f:l)] == collect(-10:13)
    @test A[Not(f:l)] == A[Not(:)] == A[[]] == []
    @test A[Not([f,f,f,f+1,f+1])] == A[Not(f,f+1,f,f+1,f)] == A[f+2:l]
    @test A[Not([l,l,l,l-1,l-1])] == A[Not(l-1,l,l,l-1,l-1)] == A[f:l-2]

    @test_throws BoundsError A[Not(f-1)]
    @test_throws BoundsError A[Not(l+1)]
end

@testset "1-d read/write" for A in (collect(1:4), reshape(collect(1:4),2,2))
    A[Not(2:3)] = [44, 11]
    @test vec(A) == [44, 2, 3, 11]
    A[Not(2:4)] .= 0
    @test vec(A) == [0, 2, 3, 11]
    A[Not(1:end)] .= 100
    @test vec(A) == [0, 2, 3, 11]
    A[Not(:)] .= 100
    @test vec(A) == [0, 2, 3, 11]

    A[Not([])] .= 0
    @test all(A.== 0)

    B = copy(A)
    @test_throws BoundsError A[Not(0)] = 200
    @test_throws BoundsError A[Not(end+1)] = 300
    @test A == B
end

@testset "2-d readonly" for A in (reshape(-10:13,3,:), reshape(-10:13,3,:,1))
    @test A[Not(2), :] == (@view A[Not(2), :]) == A[[1,3],:]
    @test A[:, Not(2)] == (@view A[:, Not(2)]) == A[:,[1;3:end]]
    @test A[Not(2), Not(2)] == (@view A[Not(2), Not(2)]) == A[[1;3:end],[1;3:end]]
    R = collect(CartesianIndices(size(A)))
    @test A[Not(first(R))] == (@view A[Not(first(R))]) == A[2:end]
    @test A[Not(R[1:2])] == (@view A[Not(R[1:2])]) == A[3:end]
    @test A[Not(iseven.(A))] == (@view A[Not(iseven.(A))]) == A[isodd.(A)] == collect(-9:2:13)
    @test A[Not([1,1,1,2,2]), Not([2,2,1,1,1])] == A[Not(2,1,1,2), Not(2,2,1,1,2)] == A[3:end, 3:end]
    @test A[Not([end,end-1,end,end-1]), Not([end,end-1,end,end-1])] == A[Not(end,end-1,end,end-1), Not(end,end-1,end,end-1)] == A[1:end-2, 1:end-2]
    @test A[Not([2,2,1,1,1]), Not([end,end-1,end,end-1])] == A[Not(2,2,1,1,2), Not(end,end-1,end,end-1)] == A[3:end, 1:end-2]
end

@testset "2-d offset readonly" for A in (OffsetArray(reshape(-10:13,3,:), -1, 0), OffsetArray(reshape(-10:13,3,4,:), -10, 20, -30),)
    inds = axes(A)
    f₁, l₁ = first(inds[1]), last(inds[1])
    f₂, l₂ = first(inds[2]), last(inds[2])
    # TODO: Re-enable these tests for 3-d after PLI deprecation
    if ndims(A) == 2
        @test A[Not(f₁+1), f₂:l₂] == (@view A[Not(f₁+1), f₂:l₂]) == A[[f₁,l₁],f₂:l₂]
        @test A[f₁:l₁, Not(f₂+1)] == (@view A[f₁:l₁, Not(f₂+1)]) == A[f₁:l₁,[f₂;f₂+2:l₂]]
        @test A[Not(f₁+1), Not(f₂+1)] == (@view A[Not(f₁+1), Not(f₂+1)]) == A[[f₁,l₁],[f₂;f₂+2:l₂]]
    end
    R = collect(CartesianIndices(axes(A)))
    @test A[Not(first(R))] == (@view A[Not(first(R))]) == A[LinearIndices(A)[2:end]]
    @test A[Not(R[1:2])] == (@view A[Not(R[1:2])]) == A[LinearIndices(A)[3:end]]
    @test A[Not(iseven.(A))] == (@view A[Not(iseven.(A))]) == A[isodd.(A)] == collect(-9:2:13)
end

@testset "Utilities" begin
    @test sprint(show, InvertedIndices.InvertedIndexIterator([2,4], 1:5)) == "Not([2, 4])"
    A = [1]
    ex = try A[Not(2,3)] catch ex; ex end
    @test occursin("Not([2, 3])", sprint(Base.showerror, ex))
end

@testset "broadcasting" begin
    # these scenarios should be later properly handled by DataFrames.jl by expanding first dimension as appropriate
    @test (Not(1) .=> sin) == (InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin)
    @test (Not(1) .=> [sin cos]) ==
          [InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos]
    @test (Not(1) .=> [sin cos] .=> [:a, :b, :c]) ==
          [InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin => :a InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos => :a
           InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin => :b InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos => :b
           InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin => :c InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos => :c]
    @test (Not(1) .=> [sin cos] .=> Not(1)) ==
          [InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin => InvertedIndices.BroadcastedInvertedIndex(Not(1)) InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos => InvertedIndices.BroadcastedInvertedIndex(Not(1))]
    @test (Not(1) .=> [sin, cos]) ==
          [InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin, InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos]
    @test (Not(1) .=> [sin, cos] .=> [:a, :b]) ==
          [InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin => :a, InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos => :b]
    @test (Not(1) .=> [sin, cos] .=> Not(1)) ==
          [InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin => InvertedIndices.BroadcastedInvertedIndex(Not(1)), InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos => InvertedIndices.BroadcastedInvertedIndex(Not(1))]

    @test identity.(Not(1)) == InvertedIndices.BroadcastedInvertedIndex(Not(1))
    @test (identity.(Not(1)) .=> sin) == (InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin)
    @test (identity.(Not(1)) .=> [sin cos]) ==
          [InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos]
    @test (identity.(Not(1)) .=> [sin cos] .=> [:a, :b, :c]) ==
          [InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin => :a InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos => :a
           InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin => :b InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos => :b
           InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin => :c InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos => :c]
    @test (identity.(Not(1)) .=> [sin cos] .=> Not(1)) ==
          [InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin => InvertedIndices.BroadcastedInvertedIndex(Not(1)) InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos => InvertedIndices.BroadcastedInvertedIndex(Not(1))]
    @test (identity.(Not(1)) .=> [sin, cos]) ==
          [InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin, InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos]
    @test (identity.(Not(1)) .=> [sin, cos] .=> [:a, :b]) ==
          [InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin => :a, InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos => :b]
    @test (identity.(Not(1)) .=> [sin, cos] .=> Not(1)) ==
          [InvertedIndices.BroadcastedInvertedIndex(Not(1)) => sin => InvertedIndices.BroadcastedInvertedIndex(Not(1)), InvertedIndices.BroadcastedInvertedIndex(Not(1)) => cos => InvertedIndices.BroadcastedInvertedIndex(Not(1))]
end

@testset "named tuples" begin
    nt = (a=1, b="2", c=:c)
    @test nt[Not(:a)] === (b="2", c=:c)
    @test nt[Not((:a,))] === (b="2", c=:c)
    @test nt[Not((:a, :a))] === (b="2", c=:c)
    @test nt[Not((:a, :b))] === (c=:c,)
    @test nt[Not([:a, :b])] === (c=:c,)
    @test nt[Not(())] === nt
    @test_throws ErrorException nt[Not(:xxx)] === nt
    @test_throws ErrorException nt[Not((:a, :xxx))] === nt

    if VERSION >= v"1.1"
        @test_broken @inferred nt[Not(:a)]
        @test_broken @inferred nt[Not((:a, :b))]
        f = nt -> nt[Not(:a)]
        @test @inferred(f(nt)) === (b="2", c=:c)
        f = nt -> nt[Not((:a, :b))]
        @test @inferred(f(nt)) === (c=:c,)
    end
end


text:
todisplay = VERSION <= VersionNumber("1.2") ? true : false

using CAGD
using SparseArrays
using LinearAlgebra
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
if todisplay
    using ViewerGL
	GL = ViewerGL
	include("views.jl")
end

V = 2 * [
    0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0
    0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0
    0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0
]
EV = [
    [1, 2], [1, 3], [2, 4], [3, 4],
    [1, 5], [2, 6], [3, 7], [4, 8],
    [5, 6], [5, 7], [6, 8], [7, 8]
]
FE = [
    [1, 2, 3, 4], [9, 10, 11, 12],
    [1, 5, 6, 9], [2, 5, 7, 10], [3, 6, 8, 11], [4, 7, 8, 12]
]
CF = [
    [1, 2, 3, 4, 5, 6]
]

m1 = CAGD.Model(V)
CAGD.addModelCells!(m1, 1, EV, signed = true)
CAGD.addModelCells!(m1, 2, FE, signed = true)
CAGD.addModelCells!(m1, 3, CF, signed = false)

m2 = deepcopy(m1)
m2.G .+= [1.0; 1.0; 1.0]

model = CAGD.uniteModels(m1, m2)

if todisplay  displayModel(model)  end

split_model = CAGD.facesplitting(model)

if todisplay  displayModel(split_model)  end

congr_model = CAGD.mergeModelVertices(split_model, signed_merge=true)

if todisplay  displayModel(congr_model)  end

gift_model = deepcopy(congr_model)

FCcycles, bicon_comps = CAGD.tgw(congr_model, 3)

externals = [bc[CAGD.getExternalBoundary(congr_model, 3, FCcycles, bc)] for bc in bicon_comps]
internals = [setdiff(bicon_comps[i], externals[i]) for i = 1 : length(externals)]

FC = FCcycles[:, internals...]
JK = ones(Int8, length(externals))
FC = [FCcycles * SparseArrays.sparse(externals, JK, JK, size(FCcycles, 2), 1) FC]

CAGD.addModelCells!(gift_model, 3, convert(Lar.ChainOp, FC'))

if todisplay  viewExplode(gift_model)  end

# Each face is percurred two times
sum(gift_model.T[3], dims = 1)

boolmatrix, arranged_model = CAGD.bool3(model)

arranged_model == gift_model


text:
immutable InverseGamma <: ContinuousUnivariateDistribution
    invd::Gamma
    β::Float64

    function InverseGamma(α::Real, β::Real)
        (α > zero(α) && β > zero(β)) || error("Both shape and scale must be positive.")
        new(Gamma(α, 1.0 / β), float64(β))
    end

    InverseGamma(α::Real) = InverseGamma(α, 1.0)
    InverseGamma() = InverseGamma(1.0, 1.0)
end

@distr_support InverseGamma 0.0 Inf

#### Parameters

shape(d::InverseGamma) = shape(d.invd)
scale(d::InverseGamma) = d.β
rate(d::InverseGamma) = scale(d.invd)

params(d::InverseGamma) = (shape(d), scale(d))

#### Parameters

mean(d::InverseGamma) = ((α, β) = params(d); α  > 1.0 ? β / (α - 1.0) : Inf)

mode(d::InverseGamma) = scale(d) / (shape(d) + 1.0)

function var(d::InverseGamma)
    (α, β) = params(d)
    α > 2.0 ? β^2 / ((α - 1.0)^2 * (α - 2.0)) : Inf
end

function skewness(d::InverseGamma)
    α = shape(d)
    α > 3.0 ? 4.0 * sqrt(α - 2.0) / (α - 3.0) : NaN
end

function kurtosis(d::InverseGamma)
    α = shape(d)
    α > 4.0 ? (30.0 * α - 66.0) / ((α - 3.0) * (α - 4.0)) : NaN
end

function entropy(d::InverseGamma)
    (α, β) = params(d)
    α + lgamma(α) - (1.0 + α) * digamma(α) + log(β)
end

#### Evaluation

pdf(d::InverseGamma, x::Float64) = exp(logpdf(d, x))

function logpdf(d::InverseGamma, x::Float64)
    (α, β) = params(d)
    α * log(β) - lgamma(α) - (α + 1.0) * log(x) - β / x
end

cdf(d::InverseGamma, x::Float64) = ccdf(d.invd, 1.0 / x)
ccdf(d::InverseGamma, x::Float64) = cdf(d.invd, 1.0 / x)
logcdf(d::InverseGamma, x::Float64) = logccdf(d.invd, 1.0 / x)
logccdf(d::InverseGamma, x::Float64) = logcdf(d.invd, 1.0 / x)

quantile(d::InverseGamma, p::Float64) = 1.0 / cquantile(d.invd, p)
cquantile(d::InverseGamma, p::Float64) = 1.0 / quantile(d.invd, p)
invlogcdf(d::InverseGamma, p::Float64) = 1.0 / invlogccdf(d.invd, p)
invlogccdf(d::InverseGamma, p::Float64) = 1.0 / invlogcdf(d.invd, p)

function mgf(d::InverseGamma, t::Real)
    (a, b) = params(d)
    t == zero(t) ? one(float(t)) : 2.0*(-b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4.0*b*t))
end

function cf(d::InverseGamma, t::Real)
    (a, b) = params(d)
    t == zero(t) ? complex(one(float(t))) : 2.0*(-im*b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4.0*im*b*t))
end

#### Evaluation

rand(d::InverseGamma) = 1.0 / rand(d.invd)

function _rand!(d::InverseGamma, A::AbstractArray)
    s = sampler(d.invd)
    for i = 1:length(A)
    	v = 1.0 / rand(s)
        @inbounds A[i] = v
    end
    A
end