text:
export
    revit,
    all_walls,
    all_floors,
    RevitSystemFamily,
    RevitFileFamily,
    RevitInPlaceFamily,
    revit_systems_family,
    revit_file_family

#=
We need to ensure the Revit plugin is properly installed.
For Revit, there are a few places where plugins can be installed:

User Addins:
%appdata%\Autodesk\Revit\Addins\<version>\
%appdata%\Autodesk\ApplicationPlugins\

Machine Addins (for all users of the machine):
C:\ProgramData\Autodesk\Revit\Addins\<version>\

Addins packaged for the Autodesk Exchange store:
C:\ProgramData\Autodesk\ApplicationPlugins\

Autodesk servers and services:
C:\Program Files\Autodesk\Revit <version>\AddIns\

C:\Users\<username>\AppData\Roaming\Autodesk\Revit \Addins\<year>

=#

julia_khepri = dirname(dirname(abspath(@__FILE__)))
# 1. The dlls are updated in VisualStudio after compilation of the plugin.
plugin_name = "KhepriRevit"
khepri_dlls = ["KhepriBase.dll", plugin_name*".dll"]
# 2. Depending on whether we are in Debug mode or Release mode,
development_phase = "Debug" # "Release"
# 3. the dlls are located in a folder
dlls_folder = joinpath("bin", "x64", development_phase)
# 4. contained inside the Plugins folder, which has a specific location regarding this file itself
plugin_folder = joinpath(dirname(dirname(dirname(dirname(abspath(@__FILE__))))), "Plugins", plugin_name, plugin_name)
# 5. Besides the dlls, we also need the bundle folder
bundle_name = plugin_name*".bundle"
bundle_dll_folder = joinpath(bundle_name, "Contents")
# 6. which is contained in the Plugins folder
bundle_path = joinpath(plugin_folder, bundle_name)
pkg_cnts_name = "PackageContents.xml"
local_plugins = joinpath(dirname(dirname(abspath(@__FILE__))), "Plugin")
local_khepri_plugin = joinpath(local_plugins, bundle_name)
local_khepri_plugin_dll_folder = joinpath(local_plugins, bundle_dll_folder)

# This only needs to be done when the Revit plugin is updated
upgrade_plugin(; advance_major_version=false, advance_minor_version=true) =
  begin
    # Update major or minor version
    if advance_major_version || advance_minor_version
      bundle_xml = joinpath(bundle_path, pkg_cnts_name)
      doc = readxml(bundle_xml)
      app_pkg = findfirst("//ApplicationPackage", doc)
      major, minor = map(s -> parse(Int, s), split(app_pkg["AppVersion"], '.'))
      print("Advancing version from $(major).$(minor) ")
      major += advance_major_version ? 1 : 0
      minor += advance_minor_version ? 1 : 0
      println("to $(major).$(minor).")
      app_pkg["AppVersion"] = "$(major).$(minor)"
      write(bundle_xml, doc)
    end
    # 7. The bundle needs to be copied to the current folder
    local_bundle_path = joinpath(julia_khepri, "Plugin", bundle_name)
    # 8. but, before, we remove any previously existing bundle
    mkpath(dirname(local_bundle_path))
    rm(local_bundle_path, force=true, recursive=true)
    # 9. Now we do the copy
    cp(bundle_path, local_bundle_path)
    # 10. and we copy the dlls to the local bundle Contents folder
    local_bundle_contents_path = joinpath(local_bundle_path, "Contents")
    for dll in khepri_dlls
      src = joinpath(plugin_folder, dlls_folder, dll)
      dst = joinpath(local_bundle_contents_path, dll)
      rm(dst, force=true)
      cp(src, dst)
    end
  end

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

env(name) = Sys.iswindows() ? ENV[name] : ""

revit_general_plugins = joinpath(dirname(env("CommonProgramFiles")), "Autodesk", "ApplicationPlugins")
revit_allusers_plugins = joinpath(env("ALLUSERSPROFILE"), "Autodesk", "ApplicationPlugins")
revit_user_plugins = joinpath(env("APPDATA"), "Autodesk", "ApplicationPlugins")
revit_khepri_plugin = joinpath(revit_user_plugins, bundle_name)
revit_khepri_plugin_dll_folder = joinpath(revit_user_plugins, bundle_dll_folder)

revit_version(path) =
  let doc = readxml(path),
      app_pkg = findfirst("//ApplicationPackage", doc)
    VersionNumber(map(s -> parse(Int, s), split(app_pkg["AppVersion"], '.'))...)
  end

update_plugin() =
  let local_path_xml = joinpath(local_khepri_plugin, pkg_cnts_name)
      revit_path_xml = joinpath(revit_khepri_plugin, pkg_cnts_name)
    # Do we have the bundle folder?
    isdir(revit_khepri_plugin) || mkpath(revit_khepri_plugin)
    isdir(revit_khepri_plugin_dll_folder) || mkpath(revit_khepri_plugin_dll_folder)
    # Must we upgrade?
    need_upgrade = ! isfile(revit_path_xml) || revit_version(revit_path_xml) < revit_version(local_path_xml)
    if need_upgrade
      # remove first to avoid loosing the local file
      #isfile(revit_path_xml) && rm(revit_path_xml)
      cp(local_path_xml, revit_path_xml, force=true)
      for dll in khepri_dlls
        let path = joinpath("Contents", dll),
            local_path = joinpath(local_khepri_plugin, path),
            revit_path = joinpath(revit_khepri_plugin, path)
            # remove first to avoid loosing the local file
            #isfile(revit_path_xml) && rm(revit_path_xml)
            cp(local_path, revit_path, force=true)
        end
      end
    end
  end

checked_plugin = false

check_plugin() =
  begin
    global checked_plugin
    if ! checked_plugin
      @info("Checking Revit plugin...")
      for i in 1:10
        try
          update_plugin()
          @info("done.")
          checked_plugin = true
          return
        catch exc
          if isa(exc, Base.IOError) && i < 10
            @error("The Revit plugin is outdated! Please, close Revit.")
            sleep(5)
          else
            throw(exc)
          end
        end
      end
    end
  end

#
const revit_template = Parameter(abspath(@__DIR__, "../Plugin/KhepriTemplate.rte"))

start_revit() =
  run(`cmd /c cd "$(dirname(revit_template()))" \&\& $(basename(revit_template()))`, wait=false)

#
# RVT is a subtype of CS
parse_signature(::Val{:RVT}, sig::T) where {T} = parse_signature(Val(:CS), sig)
encode(::Val{:RVT}, t::Val{T}, c::IO, v) where {T} = encode(Val(:CS), t, c, v)
decode(::Val{:RVT}, t::Val{T}, c::IO) where {T} = decode(Val(:CS), t, c)

#
# We need some additional Encoders
encode(::Val{:RVT}, t::Union{Val{:XYZ},Val{:VXYZ}}, c::IO, p) =
  encode(Val(:CS), Val(:double3), c, raw_point(p))
decode(::Val{:RVT}, t::Val{:XYZ}, c::IO) =
  xyz(decode(Val(:CS), Val(:double3), c)..., world_cs)
decode(::Val{:RVT}, t::Val{:VXYZ}, c::IO) =
  vxyz(decode(Val(:CS), Val(:double3), c)..., world_cs)
encode(ns::Val{:RVT}, t::Union{Val{:ElementId},Val{:Element},Val{:Level},Val{:FloorFamily}}, c::IO, v) =
  encode(ns, Val(:int), c, v)
decode(ns::Val{:RVT}, t::Union{Val{:ElementId},Val{:Element},Val{:Level},Val{:FloorFamily}}, c::IO) =
  decode_or_error(ns, Val(:int), c, Int32(-1234))

@encode_decode_as(:RVT, Val{:Length}, Val{:double})

revit_api = @remote_functions :RVT """
public Element Sphere(XYZ centre, Length radius)
public Element ConeFrustumNamed(string name, XYZ bottom, VXYZ axis, Length bottomRadius, Length height, Length topRadius)
public Element ConeFrustum(XYZ bottom, VXYZ axis, Length bottomRadius, Length height, Length topRadius)
public Element Cylinder(XYZ bottom, VXYZ axis, Length radius, Length height)
public Element Cone(XYZ bottom, VXYZ axis, Length bottomRadius, Length height)
public ElementId SurfaceFromGrid(int m, int n, XYZ[] pts, bool closedM, bool closedN, int level)
public Element PyramidFrustumNamed(String name, XYZ[] ps, XYZ[] qs, ElementId materialId)
public Element PyramidFrustumWithMaterial(XYZ[] ps, XYZ[] qs, ElementId materialId)
public Element PyramidFrustum(XYZ[] ps, XYZ[] qs)
public Element ExtrudedContourNamed(string name, XYZ[] contour, bool smoothContour, XYZ[][] holes, bool[] smoothHoles, XYZ v, ElementId materialId)
public Element ExtrudedContourWithMaterial(XYZ[] contour, bool smoothContour, XYZ[][] holes, bool[] smoothHoles, VXYZ v, ElementId materialId)
public Element ExtrudedContour(XYZ[] contour, bool smoothContour, XYZ[][] holes, bool[] smoothHoles, VXYZ v)
public Element SurfaceGrid(XYZ[] linearizedMatrix, int n, int m)
public Element Union(ElementId idA, ElementId idB)
public Element Intersection(ElementId idA, ElementId idB)
public Element Subtraction(ElementId idA, ElementId idB)
public Level FindOrCreateLevelAtElevation(Length elevation)
public Level UpperLevel(Level level, Length addedElevation)
public Length GetLevelElevation(Level level)
public ElementId LoadFamily(string fileName)
public ElementId FamilyElement(ElementId familyId, string[] names, Length[] values)
public String InstalledLibraryPath(String root)
public void MoveElement(ElementId id, XYZ translation)
public void RotateElement(ElementId id, double angle, XYZ axis0, XYZ axis1)
public ElementId CreatePolygonalFloor(XYZ[] pts, ElementId levelId)
public ElementId CreatePolygonalRoof(XYZ[] pts, ElementId levelId, ElementId famId)
public ElementId CreatePathFloor(XYZ[] pts, double[] angles, ElementId levelId)
public ElementId CreatePathRoof(XYZ[] pts, double[] angles, ElementId levelId, ElementId famId)
public Element InsertDoor(Length deltaFromStart, Length deltaFromGround, Element host, ElementId familyId)
public Element InsertWindow(Length deltaFromStart, Length deltaFromGround, Element host, ElementId familyId, string[] names, object[] values)
public Element InsertRailing(Element host, ElementId familyId)
public void CreateFamily(string familyTemplatesPath, string familyTemplateName, string familyName)
public void CreateFamilyExtrusionTest(XYZ[] pts, double height)
public void InsertFamily(string familyName, XYZ p)
public void CreatePolygonalOpening(XYZ[] pts, Element host)
public void CreatePathOpening(XYZ[] pts, double[] angles, Element host)
public ElementId CreateBeam(XYZ p0, XYZ p1, double rotationAngle, ElementId famId)
public Element CreateColumn(XYZ location, ElementId baseLevelId, ElementId topLevelId, ElementId famId)
public Element CreateColumnPoints(XYZ p0, XYZ p1, Level level0, Level level1, ElementId famId)
public ElementId[] CreateLineWall(XYZ[] pts, ElementId baseLevelId, ElementId topLevelId, ElementId famId)
public ElementId[] CreateUnconnectedLineWall(XYZ[] pts, ElementId baseLevelId, double height, ElementId famId)
public ElementId CreateSplineWall(XYZ[] pts, ElementId baseLevelId, ElementId topLevelId, ElementId famId, bool closed)
public ElementId CreateSplineCurtainWall(XYZ[] pts, ElementId baseLevelId, ElementId topLevelId, ElementId famId, bool closed)
public Element CreateLineRailing(XYZ[] pts, ElementId baseLevelId, ElementId familyId)
public Element CreatePolygonRailing(XYZ[] pts, ElementId baseLevelId, ElementId familyId)
public Level[] DocLevels()
public Element[] DocElements()
public Element[] DocFamilies()
public Element[] DocFloors()
public Element[] DocCeilings()
public Element[] DocWalls()
public Element[] DocWallsAtLevel(Level level)
public XYZ[] LineWallVertices(Element element)
public ElementId ElementLevel(Element element)
public ElementId WallTopLevel(Element element)
public double WallHeight(Element element)
public void HighlightElement(ElementId id)
public ElementId[] GetSelectedElements()
public bool IsProject()
public void DeleteAllElements()
public void SetView(XYZ camera, XYZ target, double focal_length)
public void EnergyAnalysis()
"""

abstract type RVTKey end
const RVTId = Int32
const RVTIds = Vector{RVTId}
const RVTRef = GenericRef{RVTKey, RVTId}
const RVTRefs = Vector{RVTRef}
const RVTEmptyRef = EmptyRef{RVTKey, RVTId}
const RVTUniversalRef = UniversalRef{RVTKey, RVTId}
const RVTNativeRef = NativeRef{RVTKey, RVTId}
const RVTUnionRef = UnionRef{RVTKey, RVTId}
const RVTSubtractionRef = SubtractionRef{RVTKey, RVTId}
const RVT = SocketBackend{RVTKey, RVTId}
const RVTVoidId = RVTId(-1)

KhepriBase.void_ref(b::RVT) = RVTNativeRef(RVTVoidId)

KhepriBase.has_boolean_ops(::Type{RVT}) = HasBooleanOps{true}()

#
KhepriBase.before_connecting(b::RVT) =
  check_plugin()
KhepriBase.after_connecting(b::RVT) =
  begin
    # C:\\ProgramData\\Autodesk\\RVT 2017\\Libraries\\US Metric\\
    set_backend_family(default_wall_family(), revit, revit_system_family())
    set_backend_family(default_window_family(), revit, revit_system_family())
    set_backend_family(default_slab_family(), revit, revit_system_family())
    set_backend_family(default_column_family(), revit, revit_file_family(
          revit_library_path("Imperial Library", raw"../US Metric/Structural Columns/Concrete/M_Concrete-Rectangular-Column.rfa"),
          ["b"=>f->f.profile.dx, "h"=>f->f.profile.dy]))
    set_backend_family(default_beam_family(), revit, revit_file_family(
          revit_library_path("Imperial Library", raw"../US Metric/Structural Framing/Wood/M_Timber.rfa"),
          ["b"=>f->f.profile.dx, "d"=>f->f.profile.dy]))
    set_backend_family(default_truss_bar_family(), revit, revit_file_family(
          revit_library_path("Imperial Library", raw"../US Metric/Structural Framing/Steel/M_W-Wide Flange.rfa")))
    set_backend_family(default_truss_node_family(), revit, revit_file_family(
          revit_library_path("Imperial Library", raw"../US Metric/Structural Framing/Steel/M_W-Wide Flange.rfa")))
  #=
    #set_backend_family(default_column_family(), unity, unity_material_family("Materials/Concrete/Concrete2"))
    #set_backend_family(default_door_family(), unity, unity_material_family("Materials/Wood/InteriorWood2"))

    =#
    set_backend_family(default_panel_family(), revit, revit_system_family())
  end

const revit = RVT("Revit", revit_port, revit_api)

# Levels

realize(b::RVT, s::Level) =
  @remote(b, FindOrCreateLevelAtElevation(s.height))

# Families
#=

Revit families are divided into
1. System Families (for walls, roofs, floors, pipes)
2. Loadable Families (for building components that have an associated file)
3. In-Place Families (for unique elements created just for the current project)

=#

abstract type RevitFamily <: Family end

struct RevitSystemFamily <: RevitFamily
    family_map::Dict{String, Function}
    instance_map::Dict{String, Function}
    instance_ref::Parameter{RVTId}
end

revit_system_family(family_map=(), instance_map=()) =
    RevitSystemFamily(
        Dict(family_map...),
        Dict(instance_map...),
        Parameter(RVTId(0))) # instead of RVTVoidId.  We need to think this carefully.

backend_get_family_ref(b::RVT, f::Family, rvtf::RevitSystemFamily) =
  begin
    if rvtf.instance_ref()===RVTVoidId
      let param_map = rvtf.family_map,
          params = keys(param_map)
        rvtf.instance_ref(
            @remote(b, FamilyElement(
                0,
                collect(keys),
                [param_map[param](f) for param in params])))
      end
    end
    rvtf.instance_ref()
  end

struct RevitFileFamily <: RevitFamily
    path::String
    family_map::Dict{String, Function}
    instance_map::Dict{String, Function}
    family_ref::Parameter{RVTId}
    instance_ref::Parameter{RVTId}
end

revit_file_family(path, family_map=(), instance_map=()) =
    RevitFileFamily(
        path,
        Dict(family_map...),
        Dict(instance_map...),
        Parameter(RVTVoidId),
        Parameter(RVTVoidId))

backend_get_family_ref(b::RVT, f::Family, rvtf::RevitFileFamily) =
  begin
    if true #rvtf.family_ref()===RVTVoidId
      rvtf.family_ref(@remote(b, LoadFamily(rvtf.path)))
    end
    if true #rvtf.instance_ref()===RVTVoidId
      let param_map = rvtf.family_map,
          params = keys(param_map)
        rvtf.instance_ref(
            @remote(b, FamilyElement(
                rvtf.family_ref(),
                collect(params),
                [param_map[param](f) for param in params])))
      end
    end
    rvtf.instance_ref()
  end
#

# This is for future use
struct RevitInPlaceFamily <: RevitFamily
    parameter_map::Dict{Symbol,String}
    ref::Parameter{Int}
end

#=
root should be "Imperial Library" or "Metric Library"
path can be something as "Structural Framing\\Wood\\M_Timber.rfa"
=#
export revit_library_path
revit_library_path(root::String, path::String) =
  joinpath(@remote(revit, InstalledLibraryPath(root)), path)

switch_to_backend(from::Backend, to::RVT) =
    let height = level_height(default_level())
        current_backend(to)
        default_level(level(height))
    end

#=
realize(b::RVT, f::TableFamily) =
    @remote(b, CreateRectangularTableFamily(f.length, f.width, f.height, f.top_thickness, f.leg_thickness))
realize(b::RVT, f::ChairFamily) =
    @remote(b, CreateChairFamily(f.length, f.width, f.height, f.seat_height, f.thickness))
realize(b::RVT, f::TableChairFamily) =
    @remote(b, CreateRectangularTableAndChairsFamily(
        ref(f.table_family), ref(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)))

backend_rectangular_table(b::RVT, c, angle, family) =
    @remote(b, Table(c, angle, ref(family)))

backend_chair(b::RVT, c, angle, family) =
    @remote(b, Chair(c, angle, ref(family)))

backend_rectangular_table_and_chairs(b::RVT, c, angle, family) =
    @remote(b, TableAndChairs(c, angle, ref(family)))
=#

realize(b::RVT, s::EmptyShape) =
  RVTEmptyRef()
realize(b::RVT, s::UniversalShape) =
  RVTUniversalRef()

KhepriBase.b_slab(b::RVT, profile::Region, level, family) =
  let outer = outer_path(profile),
      inners = inner_paths(profile),
      slab_r = b_slab(b, outer, level, family)
    for inner in inners
      create_slab_opening(b, inner, slab_r)
    end
    slab_r
  end

KhepriBase.b_slab(b::RVT, contour::ClosedPolygonalPath, level, family) =
  begin
    @remote(b, CreatePolygonalFloor(convert(ClosedPolygonalPath, contour).vertices, ref(b, level).value))
    # we are not using the family yet
    # ref(b, s.family))
  end

create_slab_opening(b::RVT, contour::ClosedPolygonalPath, slab_r) =
  @remote(b, CreatePolygonalOpening(convert(ClosedPolygonalPath, contour).vertices, slab_r))

create_slab_opening(b::RVT, contour::ClosedPath, slab_r) =
  let (locs, arcs) = locs_and_arcs(contour)
    @remote(b, CreatePathOpening(locs, arcs, slab_r))
  end

locs_and_arcs(arc::ArcPath) =
    ([arc.center + vpol(arc.radius, arc.start_angle)],
     [arc.amplitude])

locs_and_arcs(circle::CircularPath) =
    let (locs1, arcs1) = locs_and_arcs(arc_path(circle.center, circle.radius, 0, pi))
        (locs2, arcs2) = locs_and_arcs(arc_path(circle.center, circle.radius, pi, pi))
        ([locs1..., locs2...], [arcs1..., arcs2...])
    end

KhepriBase.b_slab(b::RVT, contour::ClosedPath, level, family) =
  let (locs, arcs) = locs_and_arcs(contour)
    @remote(b, CreatePathFloor(locs, arcs, ref(b, level).value))
    # we are not using the family yet
    # ref(b, s.family))
  end

KhepriBase.b_roof(b::RVT, contour::ClosedPath, level, family) =
  let (locs, arcs) = locs_and_arcs(contour)
    @remote(b, CreatePathRoof(locs, arcs, ref(b, level).value, family))
  end

#Beams are aligned along the top axis.
KhepriBase.b_beam(b::RVT, c, h, angle, family) =
  @remote(b, CreateBeam(c, add_z(c, h), angle, realize(b, family)))

KhepriBase.b_column(b::RVT, cb, angle, bottom_level, top_level, family) =
  @remote(b, CreateColumn(cb, realize(b, bottom_level), realize(b, top_level), realize(b, family)))

#Columns are aligned along the center axis.
KhepriBase.b_free_column(b::RVT, cb, h, angle, family) =
  let ct = in_world(add_z(cb, h)),
      cb = in_world(cb),
      lb = @remote(b, FindOrCreateLevelAtElevation(cb.z)),
      lt = @remote(b, FindOrCreateLevelAtElevation(ct.z))
    @remote(b, CreateColumnPoints(cb, ct, lb, lt, realize(b, family)))
  end

KhepriBase.realize_wall_no_openings(b::RVT, s::Wall) =
  # Revit also considers unconnected walls. These have a top level with id -1
  if ref(b, s.top_level).value == RVTVoidId
      @remote(b, CreateUnconnectedLineWall(
          convert(OpenPolygonalPath, s.path).vertices,
          ref(b, s.bottom_level).value,
          s.top_level.height - s.bottom_level.height,
          realize(b, s.family)))
  else
      @remote(b, CreateLineWall(
          convert(OpenPolygonalPath, s.path).vertices,
          ref(b, s.bottom_level).value,
          ref(b, s.top_level).value,
          realize(b, s.family)))
  end

realize_wall_openings(b::RVT, w::Wall, w_ref, openings) =
  begin
      for opening in openings
          realize(b, opening)
      end
      w_ref
  end

realize(b::RVT, s::Window) =
  let rvtf = backend_family(b, s.family),
      param_map = rvtf.instance_map,
      params = keys(param_map)
    @remote(b, InsertWindow(
        s.loc.x,
        s.loc.y,
        ref(b, s.wall).value,
        backend_get_family_ref(b, s.family, rvtf),
        collect(params),
        [param_map[param](s.family) for param in params]))
  end

backend_add_door(b::RVT, w::Wall, loc::Loc, family::DoorFamily) = finish_this()
backend_add_window(b::RVT, w::Wall, loc::Loc, family::WindowFamily) =
    let d = window(w, loc, family=family)
        push!(w.windows, d)
        if realized(w) && ! realized(d)
            realize(b, d)
        end
        w
    end

realize(b::RVT, s::TrussNode) =
  @remote(b, CreateBeam(s.p, add_x(s.p, 0.1), 0, realize(b, s.family)))

realize(b::RVT, s::TrussBar) =
  @remote(b, CreateBeam(s.p0, s.p1, s.angle, realize(b, s.family)))

############################################
# Select New Family ...
# Choose Metric Generic Model

# Revit does not use materials!
KhepriBase.material_ref(b::RVT, m::Material) = nothing

KhepriBase.b_pyramid(b::RVT, bs, t, bmat, smat) =
  @remote(b, Pyramid(bs, t))
KhepriBase.b_pyramid_frustum(b::RVT, bs, ts, bmat, tmat, smat) =
  @remote(b, PyramidFrustum(bs, ts))
backend_right_cuboid(b::RVT, cb, width, height, h, material) =
  @remote(b, CenteredBox(cb, width, height, h))
KhepriBase.b_box(b::RVT, c, dx, dy, dz, mat) =
  @remote(b, Box(c, dx, dy, dz))
KhepriBase.b_cone(b::RVT, cb, r, h, bmat, smat) =
  @remote(b, Cone(cb, vz(1, cb.cs), r, h))
KhepriBase.b_cone_frustum(b::RVT, cb, rb, h, rt, bmat, tmat, smat) =
  @remote(b, ConeFrustum(cb, vz(1, cb.cs), rb, h, rt))
KhepriBase.b_cylinder(b::RVT, cb, r, h, bmat, tmat, smat) =
  @remote(b, Cylinder(cb, vz(1, cb.cs), r, h))
#Experiment with private Element Cylinder2(XYZ bottom, VXYZ axis, Length radius, Length height) {
KhepriBase.b_sphere(b::RVT, c, r, mat) =
  @remote(b, Sphere(c, r))
realize(b::RVT, s::Torus) =
  @remote(b, Torus(s.center, vz(1, s.center.cs), s.re, s.ri))

#
realize_prism(b::RVT, top, bot, side, path::PathSet, h::Real) =
  let v = planar_path_normal(path)*h,
      contour = path.paths[1],
      holes = path.paths[2:end]
    @remote(b, ExtrudedContour(
      path_vertices(contour), is_smooth_path(contour),
      path_vertices.(holes), is_smooth_path.(holes),
      v
      ))
  end

backend_surface_grid(b::RVT, points, closed_u, closed_v, smooth_u, smooth_v) =
    @remote(b, SurfaceFromGrid(
        size(points,2),
        size(points,1),
        reshape(points,:),
        closed_u,
        closed_v,
        0))

#
unite_ref(b::RVT, r0::RVTNativeRef, r1::RVTNativeRef) =
    ensure_ref(b, @remote(b, Union(r0.value, r1.value)))

intersect_ref(b::RVT, r0::RVTNativeRef, r1::RVTNativeRef) =
    ensure_ref(b, @remote(b, Intersection(r0.value, r1.value)))

subtract_ref(b::RVT, r0::RVTNativeRef, r1::RVTNativeRef) =
    ensure_ref(b, @remote(b, Subtraction(r0.value, r1.value)))

unite_refs(b::RVT, refs::Vector{<:RVTRef}) =
    RVTUnionRef(tuple(refs...))

realize(b::RVT, s::IntersectionShape) =
  let r = foldl((r0,r1)->intersect_ref(b,r0,r1), map(ref, s.shapes),
                init=RVTUniversalRef())
    mark_deleted(s.shapes)
    r
  end
############################################

# Create 4 (3) reference planes that give the panel outline
# Create extrusion nos reference planes
# Use align tool to lock the extrusion to the reference planes
# Use draw tools to create the holes
# Change to Floor plan to give thickness
# Return to Elevation Front to visualize result
# [Don't (Save the family (requires name))] ->
# Load into project
# Place aligned with a wall

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

# Select Family Metric Mass (inside Conceptual Mass)
# Create Model
# Reference lines
# Create Form -> Solid Form
# Create Form -> Void Form
# Load into project
# Place aligned with a wall

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

KhepriBase.backend_name(b::RVT) = "Revit"

KhepriBase.b_set_view(b::RVT, camera::Loc, target::Loc, lens::Real, aperture::Real) =
  @remote(b, SetView(camera, target, lens))

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

zoom_extents(b::RVT) = @remote(b, ZoomExtents())

view_top(b::RVT) =
    @remote(b, ViewTop())

KhepriBase.b_delete_all_refs(b::RVT) =
  @remote(b, DeleteAllElements())

prompt_position(prompt::String, b::RVT) =
  let ans = @remote(b, GetPoint(prompt))
    length(ans) > 0 && ans[1]
  end

all_levels(b::RVT) =
    [level_from_ref(r, b) for r in @remote(b, DocLevels())]

level_from_ref(r, b::RVT) =
  level(r == RVTVoidId ?
          error("Level for unconnected height") :
          @remote(b, GetLevelElevation(r)),
        backend=b, ref=LazyRef(b, RVTNativeRef(r)))

unconnected_level(h::Real, b::RVT) =
    level(h, backend=b, ref=LazyRef(b, RVTNativeRef(RVTVoidId)))

all_Elements(b::RVT) =
    [element_from_ref(r, b) for r in @remote(b, DocElements())]

all_walls(b::RVT) =
    [wall_from_ref(r, b) for r in @remote(b, DocWalls())]
all_walls_at_level(level::Level, b::RVT) =
    [wall_from_ref(r, b) for r in @remote(b, DocWallsAtLevel(ref(level).value))]

wall_from_ref(r, b::RVT) =
    begin
        path = convert(Path, @remote(b, LineWallVertices(r)))
        bottom_level_id = @remote(b, ElementLevel(r))
        top_level_id = @remote(b, WallTopLevel(r))
        bottom_level = level_from_ref(bottom_level_id, b)
        top_level = top_level_id == RVTVoidId ?
                        unconnected_level(bottom_level.height + @remote(b, WallHeight(r)), b) :
                        level_from_ref(top_level_id, b)
        wall(path,
             bottom_level=bottom_level,
             top_level=top_level,
             backend=b,
             ref=LazyRef(b, RVTNativeRef(r)))
    end

#=

struct revit_family
    path::String
    map::Dict
end

struct archicad_family
    name::String
    map::Dict
end

# for a non-BIM backend
bars_family = beam_family(width=10,height=20,based_on=Dict(
    revit => revit_family(
        "C:\\ProgramData\\Autodesk\\RVT 2017\\Libraries\\US Metric\\Structural Framing\\Steel\\M_HSS-Hollow Structural Section.rfa",
        Dict(:width=>"b", :height=>"d", :angle=>"Cross-Section Rotation"))
#    archicad => archicad_family("SpecialBeam", Dict(:width=>"width", :height=>"height"))
))

=#


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

export optional_calibfor
function optional_calibfor end

export calibfor
function calibfor end

export CalibCatalogEntry

struct CalibCatalogEntry
    key::FileKey
    valid::Dict{Symbol, Int64}
end

function CalibCatalogEntry(p::PropDict)
    key = p.key
    valid = Dict(k => timestamp2unix(v) for (k, v) in p.valid)
    CalibCatalogEntry(key, valid)
end

Base.convert(::Type{CalibCatalogEntry}, p::PropDict) = CalibCatalogEntry(p)

export CalibCatalog

mutable struct CalibCatalog
    entries::Vector{CalibCatalogEntry}
    lut::Dict{Symbol,Vector{Pair{Int64, FileKey}}}
end

function CalibCatalog(entries::Vector{CalibCatalogEntry})
    lut = Dict{Symbol,Vector{Pair{Int64, FileKey}}}()

    for e in entries
        for (system, time) in e.valid
            if !haskey(lut, system)
                lut[system] = Vector{Pair{Int64, FileKey}}()
            end
            push!(lut[system], time => e.key)
        end
    end

    for (system, cals) in lut
        sort!(cals, by = x -> x[2].time, alg = Base.Sort.DEFAULT_STABLE)
        sort!(cals, by = x -> x[1], alg = Base.Sort.DEFAULT_STABLE)
    end

    CalibCatalog(entries, lut)
end

function Base.read(::Type{CalibCatalog}, filename::AbstractString)
    entries = open(filename) do input
        [convert(CalibCatalogEntry, convert(PropDict, l)) for l in eachline(input)]
    end
    CalibCatalog(entries)
end

function calib_available(calibs::CalibCatalog, key::FileKey, system::Symbol)
    if haskey(calibs.lut, system)
        lut = calibs.lut[system]
        cal_found = searchsortedlast(lut, key.time, by = x -> x[1])
        cal_found != 0
    else
        false
    end
end

function calibfor(calibs::CalibCatalog, key::FileKey, system::Symbol)
    lut = calibs.lut[system]
    cal_found = searchsortedlast(lut, key.time, by = x -> x[1])
    if cal_found == 0
        error("No calibration available for $key-$system")
    else
        lut[cal_found][2]
    end
end

ext: jl


text:
# https://github.com/JuliaAcademy/DataScience/blob/master/01.%20Data.ipynb
using BenchmarkTools
using DataFrames
import DataFrames
DF = DataFrames
using DelimitedFiles
using CSV
using XLSX

import JLD
import NPZ
import RData
import MAT

# Download from network
# pl: programming language
pl_csv_filename = download("https://raw.githubusercontent.com/nassarhuda/easy_data/master/programming_languages.csv",
     "programming_languages.csv")
# This style save to a tmp file
#filename = download("https://raw.githubusercontent.com/nassarhuda/easy_data/master/programming_languages.csv")

# Read into arrays
# Body, Header
pl_body, pl_header = readdlm(pl_csv_filename, ','; header=true)

# Write to a text file with "-" as delimiter
writedlm("programming_languages_dlm.txt", pl_body, "-")

# CSV.jl
pl_csv = CSV.read(pl_csv_filename) # => DataFrame
#show(pl_csv; allrows=true)
@show typeof(pl_csv)
@show describe(pl_csv)
# column names
@show names(pl_csv)
# column
pl_csv.year
pl_csv.language

# Timing
# comment out due to slow execution
#@btime P,H = readdlm("programminglanguages.csv", ','; header=true)
#@btime C = CSV.read("programminglanguages.csv")

# Write CSV from arrays
CSV.write("programming_languages_CSV.csv", pl_csv)

# Notice the dot (.) operator below
# Create DataFrame from arrays
P1 = pl_body[:,1] # => Vector{Any}
P1 = Int64.(P1)   # => Vector{Int64}
P2 = pl_body[:,2] # => Vector{Any}
P2 = String.(P2)  # => Vector{String}
H = pl_header[:]
pl_df = DataFrame([P1, P2], H)
P1 = undef
P2 = undef
H = undef

# Excel file
# Excel file is super slow, convert to CSV
city_xlsx = XLSX.readdata("data/zillow_data_download_april2020.xlsx",
    "Sale_counts_city", # sheet name
    "A1:F9" # cell range
    )

# Does not work
# G,L = XLSX.readdata("zillow_data_download_april2020.xlsx",
#     "Sale_counts_city")
    #"A1:ER28760")

# Much faster
city_csv = CSV.read("zillow_data_download_april2020.csv")

# Build DataFrame from arrays
foodvec = ["apple", "cucumber", "tomato", "banana"]
calorievec = [105, 47, 22, 105]
pricevec = [0.85, 1.6, 0.8, 0.6]
df_calorie = DataFrame([foodvec, calorievec], ["Item", "Calorie"])
df_price = DataFrame(Item=foodvec, Price=pricevec)
# inner join
df = innerjoin(df_calorie, df_price, on="Item")
# Write to Excel file
rm("food.xlsx"; force=true) # remove old file
XLSX.writetable("food.xlsx", collect(eachcol(df)), names(df))

# Read other data format

# Julia
jld_data = JLD.load("data/mytempdata.jld")
JLD.save("mywrite.jld", "A", jld_data)
#jld2 = JLD.load("mywrite.jld")

# Numpy
npz_data = NPZ.npzread("data/mytempdata.npz")
NPZ.npzwrite("mywrite.npz", npz_data)

# R
r_data = RData.load("data/mytempdata.rda")

# Matlab
mat_data = MAT.matread("data/mytempdata.mat")
MAT.matwrite("mywrite.mat", mat_data)

# Time to process the data from Julia
function year_created(P, language::String)
    println("Array version")
    filter = P[:,2] .== language
    ind = findfirst(filter)
    if isnothing(ind)
        throw(MissingException("Language not found: $language"))
    end
    #println("loc = $loc")
    return P[ind,1]
end
year_created(pl_body, "Julia")
#year_created(P, "not exist")

function year_created(df::DataFrame, language::String)
    println("DataFrame version")
    loc = findfirst(df.language .== language)

    if isnothing(loc)
        e = MissingException("Language not found: $language")
        throw(e)
    end

    return df.year[loc]
end
year_created(pl_df, "Julia")
#year_created(pl_df, "W#")

function count_by_year(df::DataFrame, year::Int64)
    found_vec = findall(df.year .== year) # => Array{Int64,1}
    yearcount = length(found_vec)

    return yearcount
end
count_by_year(pl_df, 2011)

# Dict with Any type
D = Dict([("A", 1), ("B", 2), (1,[1,2])])

# Dict with types
pl_dic = Dict{Integer, Vector{String}}()
# Populate pl_dic
for i = 1:nrow(pl_df)
    #size(pl_df,1)
    #nrow(pl_df)
    year,lang = pl_df[i,:]
    if haskey(pl_dic, year)
        vec = pl_dic[year]
        push!(vec, lang)
    else
        pl_dic[year] = [lang]
    end
end

# dropmissing


text:
using Test
import ParallelStencil
using ParallelStencil.ParallelKernel
import ParallelStencil.ParallelKernel: @reset_parallel_kernel, @is_initialized, @get_package, @get_numbertype, SUPPORTED_PACKAGES, PKG_CUDA
import ParallelStencil.ParallelKernel: @require, @symbols, longnameof
import ParallelStencil.ParallelKernel: checkargs_init, check_already_initialized, set_initialized, is_initialized, check_initialized
using ParallelStencil.ParallelKernel.Exceptions
@static if (PKG_CUDA in SUPPORTED_PACKAGES) import CUDA end

@static for package in SUPPORTED_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @testset "1. initialization of ParallelKernel" begin
            @require !@is_initialized()
            @init_parallel_kernel($package, Float64)
            @testset "initialized" begin
                @test @is_initialized()
                @test @get_package() == $package
                @test @get_numbertype() == Float64
            end;
            @testset "Data" begin
                @test @isdefined(Data)
                @test length(@symbols($(@__MODULE__), Data)) > 1
                @test Symbol("Number") in @symbols($(@__MODULE__), Data)
                @test Symbol("Array") in @symbols($(@__MODULE__), Data)
                @test Symbol("DeviceArray") in @symbols($(@__MODULE__), Data)
            end;
            @reset_parallel_kernel()
        end;
        @testset "2. Exceptions" begin
            @testset "already initialized" begin
                set_initialized(true)
                @require is_initialized()
                @test_throws IncoherentCallError check_already_initialized()
                set_initialized(false)
            end;
            @testset "arguments" begin
                MyType = Float32
                @test_throws ArgumentError checkargs_init(5, Float64)
                @test_throws ArgumentError checkargs_init($package, MyType)
            end;
            @testset "check_initialized" begin
                @require !is_initialized()
                @test_throws NotInitializedError check_initialized()
            end;
            @reset_parallel_kernel()
        end;
    end;
)) end == nothing || true;


text:
include("solve_SDP.jl")
include("solve_minlp.jl")

struct node
    fixing::Array{Int64}
    father_lb::Float64
end

struct open_node
    fixing::Array{Int64}
    node_lb::Float64
end

function BandB_fixingsome1and0(ROPF, flag, BB_param, max_time)
    output_instance_path = ROPF.output_instance_path
    output_decomposition_path = ROPF.output_decomposition_path
    formulation = ROPF.decomposition
    instance = ROPF.instance_name
    generation = ROPF.generation
    search_strategy = BB_param.search_strategy
    branch_strategy = BB_param.branch_strategy
    seuil_u = BB_param.seuil_u
    seuil_l = BB_param.seuil_l
    #initilization
    start_time = time()
    elapsed_time = 0.0
    #SDP
    instance_dat_file_path = joinpath(output_instance_path, "$instance.dat")
    Pinput_csv_file = joinpath(generation_files_path,"$(instance)_$(generation).csv")
    outsolutionpath = joinpath("Mosek_solutions")
    solution_file = joinpath(outsolutionpath, "$(instance)_$(generation)_$(flag).dat")
    λ, Sgen_var_list, SDP_var_list, Bin_var_list, dict_quad_ctr, dict_bounds_ctr, dict_constants_ctr, dict_Bin_ctr,
     dict_MONO, dict_linear_ctr = read_dat_file(instance_dat_file_path)
    cliques_dict, CLIQUE_TREE = read_blocks(output_decomposition_path, formulation, instance)
    index_var = Dict{String, Int64}()
    nb_bin = length(Bin_var_list)
    for i=1:nb_bin
        var = Bin_var_list[i]
        index_var[var] = i
    end
    output_file = "BandB_fixingsome1_$(seuil_u)_and0_$(seuil_l)_$(search_strategy)_$(branch_strategy)_$(instance)_$(generation)_$(flag).txt"
    isdir("BandB_runs") || mkpath("BandB_runs")
    f = open(joinpath("BandB_runs", output_file), "w")
    #B&B
    mosek_pt = readdlm(joinpath(outsolutionpath, "$(instance)_$(generation)_$(flag).dat"))
    fixing_from_SDP = -ones(nb_bin)
    shunt_values = [(mosek_pt[i,2], mosek_pt[i,1]) for i in 1:size(mosek_pt,1) if mosek_pt[i,1][1:5]=="Shunt"]
    sort!(shunt_values, rev=true)
    nb_bin_unfixed = nb_bin
    for i in 1:nb_bin
        tuple = shunt_values[i]
        varname = tuple[2]
        if tuple[1] > seuil_u
            fixing_from_SDP[index_var[varname]] = 1
            nb_bin_unfixed -= 1
        elseif tuple[1] < seuil_l
            fixing_from_SDP[index_var[varname]] = 0
            nb_bin_unfixed -= 1
        end
    end
    write(f, "NB binary variables : $nb_bin_unfixed \n")
    # println(nb_bin_unfixed)
    best_ub = solve_minlp(ROPF, flag, fixing_from_SDP, index_var)
    println(best_ub)
    node0 = node(fixing_from_SDP, -Inf)
    node_list = Set([node0])
    open_nodes = Set([])
    isdir("solutions") || mkpath("solutions")
    mv("knitro_solution.csv", joinpath("solutions", "BB_solution_$(instance)_$(flag).csv"), force=true)
    write(f, "UB : $best_ub \n")
    write(f, "var=$Bin_var_list \n")
    close(f)
    if nb_bin_unfixed == 0
        println("All binary variables fixed \n")
        write(f, "BEST UB=$(best_ub) \nNb explored nodes : 0 \n")
        write(f, "Time : $(elapsed_time) \n")
        write(f, "Nb open nodes : $(length(open_nodes)) \n")
        close(f)
        return best_ub, 1, 0
    end
    log = open(joinpath("BandB_runs",output_file[1:end-3]*"csv"), "w")
    write(log, "Nb noeuds explores ; Nb noeuds a explorer ; Best UB ; Best LB ; Gap (%) ; Temps ; Nb noeuds gap optimalite \n")
    close(log)
    nb_explored_nodes = 0
    #exploring nodes
    while length(node_list) > 0 && elapsed_time < max_time
        if search_strategy == "bestfirst"
            next_node = extract_bestlb(node_list)
        elseif search_strategy == "deepfirst"
            next_node = extract_deepfirst(node_list)
        end
        nb_explored_nodes +=1
        println("Node $nb_explored_nodes \n")
        lb_father = next_node.father_lb
        fixing = copy(next_node.fixing)
        delete!(node_list, next_node)
        dict_variables_to_fix = Dict(var => fixing[index_var[var]] for var in Bin_var_list)
        value_bins, value_SDP_var, opt_value, primal_status = construct_SDP(cliques_dict, CLIQUE_TREE, Pinput_csv_file, flag, Sgen_var_list, SDP_var_list, Bin_var_list,
           dict_quad_ctr, dict_linear_ctr, dict_bounds_ctr, dict_constants_ctr, dict_Bin_ctr, dict_MONO, dict_variables_to_fix, "")
        if opt_value < (1-10^(-4))*best_ub && (primal_status == MOI.FEASIBLE_POINT || primal_status == MOI.NEARLY_FEASIBLE_POINT)
            if maximum(min(1-b, b) for (name,b) in value_bins) <= 10^(-6)
                #STOP: no nodes after
                #NOTE : vérifier car la relaxation SDP ne donne pas forcément une solution réalisable pour le QCQP non convexe
                ub = solve_minlp(ROPF, flag, fixing, index_var)
                f = open(joinpath("BandB_runs", output_file), "a")
                if (ub-opt_value)/ub > 10^(-4) #gap
                    write(f, "OPEN NODE ")
                    push!(open_nodes, open_node(fixing, opt_value))
                end
                write(f, "Node $nb_explored_nodes : fixing=$fixing ; UB=$ub ; LB=$opt_value \n")
                close(f)
                if ub < best_ub
                    best_ub = ub
                    # log = open(joinpath("BandB_runs",output_file[1:end-3]*"csv"), "a")
                    # write(log, "$nb_explored_nodes ; $(length(node_list)) ; $best_ub ; $(minimum(n.father_lb for n in node_list)) ; $((best_ub-minimum(n.father_lb for n in node_list))/best_ub*100)  ; $(length(open_nodes)) \n")
                    # close(log)
                    for op_node in node_list
                        if op_node.father_lb >= (1-10^(-4))*best_ub
                            write(f, "CUT : node fixing $(op_node.fixing) \n")
                            delete!(node_list, op_node)
                        end
                    end
                    for op_node in open_nodes
                        if op_node.node_lb >= (1-10^(-4))*best_ub
                            delete!(open_nodes, op_node)
                        end
                    end
                    mv("knitro_solution.csv", joinpath("solutions", "BB_solution_$(instance)_$(flag).csv"), force=true)
                else
                    # log = open(joinpath("BandB_runs",output_file[1:end-3]*"csv"), "a")
                    # write(log, "$nb_explored_nodes ; $(length(node_list)) ; $best_ub ; $(minimum(n.father_lb for n in node_list)) ; $((best_ub-minimum(n.father_lb for n in node_list))/best_ub*100)  ; $(length(open_nodes))  \n")
                    # close(log)
                end
            else
                f = open(joinpath("BandB_runs", output_file), "a")
                write(f, "Node $nb_explored_nodes : fixing=$fixing ;  LB=$opt_value \n")
                close(f)
                if branch_strategy == "0.5"
                    var_to_branch = select_05(value_bins)
                elseif branch_strategy == "1"
                    var_to_branch = select_1(value_bins)
                end
                    next_fixing0 = copy(fixing)
                    next_fixing0[index_var[var_to_branch]] = 0
                    node_left = node(next_fixing0, opt_value)
                    push!(node_list, node_left)
                    next_fixing1 = copy(fixing)
                    next_fixing1[index_var[var_to_branch]] = 1
                    node_right = node(next_fixing1, opt_value)
                    push!(node_list, node_right)
            end
        else
            f = open(joinpath("BandB_runs", output_file), "a")
            write(f, "Node $nb_explored_nodes : NON REALISABLE ou borne inf>best_ub fixing=$fixing ;  LB=$opt_value \n")
            close(f)
        end

        best_lb = -Inf
        if length(node_list) > 0
            best_lb = (minimum(n.father_lb for n in node_list))
        else
            best_lb = opt_value
        end
        gap = ((best_ub-best_lb)/best_ub*100)
        elapsed_time = time() - start_time
        log = open(joinpath("BandB_runs",output_file[1:end-3]*"csv"), "a")
        write(log, "$nb_explored_nodes ; $(length(node_list)) ; $best_ub ; $best_lb ; $gap ; $elapsed_time ; $(length(open_nodes)) \n")
        close(log)
    end
    f = open(joinpath("BandB_runs", output_file), "a")
    if length(node_list) == 0
        write(f, "B&B DONE \n")
    else
        write(f, "B&B stopped because of time limit, nb nodes not explored : $(length(node_list)) \n")
    end
    write(f, "BEST UB=$(best_ub) \nNb explored nodes : $nb_explored_nodes \n")
    write(f, "Time : $(elapsed_time) \n")
    write(f, "Nb open nodes : $(length(open_nodes)) \n")
    close(f)
    return best_ub, nb_explored_nodes, open_nodes
end

function extract_deepfirst(node_list)
    max_fixed_var = 0
    selected_node = first(node_list)
    for node in node_list
        nb_fixed_var = 0
        for elem in node.fixing
            if elem == 1 || elem == 0
                nb_fixed_var +=1
            end
        end
        if nb_fixed_var > max_fixed_var
            max_fixed_var = nb_fixed_var
            selected_node = node
        elseif nb_fixed_var == max_fixed_var
            nb_1_node = length([value for value in node.fixing if value==1])
            nb_1_selected_node = length([value for value in selected_node.fixing if value==1])
            if nb_1_node > nb_1_selected_node
                selected_node = node
            end
        end
    end
    return selected_node
end

function extract_bestlb(node_list)
    best_lb = Inf
    selected_node = first(node_list)
    for node in node_list
        lb = node.father_lb
        if lb < best_lb
            selected_node = node
            best_lb = lb
        end
    end
    return selected_node
end

function select_05(value_bins)
    min_distance_to_05 = 1
    var = first(keys(value_bins))
    for (varname, value) in value_bins
        distance_to_05 = abs(value-0.5)
        if distance_to_05 <= min_distance_to_05
            min_distance_to_05 = distance_to_05
            var = varname
        end
    end
    return var
end

function select_1(value_bins)
    u_values = [(value, varname) for (varname,value) in value_bins if abs(1-value) > 10.0^(-6)]
    sort!(u_values, rev=true)
    tuple = u_values[1]
    var = tuple[2]
    return var
end


text:
mutable struct TreeCallableFixedRateEngine{S <: ShortRateModel, T <: ShortRateTree} <: LatticeShortRateModelEngine{S}
  model::S
  timeSteps::Int
  common::LatticeShortRateModelEngineCommon{T}
  latticeBuilt::Bool

  function TreeCallableFixedRateEngine{S, T}(model::S,
                                            timeSteps::Int,
                                            common::LatticeShortRateModelEngineCommon{T},
                                            latticeGen::Bool = true) where {S, T}
    te = new{S, T}(model, timeSteps, common, latticeGen)

    add_observer!(model, te)

    return te
  end
end

function TreeCallableFixedRateEngine(model::S, timeSteps::Int) where {S <: ShortRateModel}
  # create empty tree
  tg = TimeGrid([1.0], 1)
  lattice = tree(model, tg)
  common = LatticeShortRateModelEngineCommon{typeof(lattice)}(tg, lattice)
  TreeCallableFixedRateEngine{S, typeof(lattice)}(model, timeSteps, common, false)
end

function _calculate!(pe::TreeCallableFixedRateEngine, bond::CallableFixedRateBond)
  tsmodel = pe.model

  refDate = reference_date(tsmodel.ts)
  dc = tsmodel.ts.dc

  callableBond = DiscretizedCallableFixedRateBond(bond, refDate, dc, pe.common.lattice.treeLattice)

  if ~pe.latticeBuilt
    times = mandatory_times(callableBond)
    tg = TimeGrid(times, pe.timeSteps)
    pe.common.tg = tg
    update!(pe)
    pe.latticeBuilt = true
  end

  # if isdefined(pe, :common)
  #   lattice = pe.common.lattice
  # else
  #   times = mandatory_times(callableBond)
  #   tg = TimeGrid(times, pe.timeSteps)
  #   lattice = tree(pe.model, tg)
  #   pe.common = LatticeShortRateModelEngineCommon(tg, lattice)
  # end

  lattice = pe.common.lattice

  redemptionTime = year_fraction(dc, refDate, callableBond.args.redemptionDate)

  initialize!(callableBond, lattice.treeLattice, redemptionTime)
  rollback!(callableBond, 0.0)
  bond.settlementValue = present_value(callableBond)

  return pe, bond
end


text:
# This function is used to measure the performance and efficiency of the interpolation functions.

using GridInterpolations
using Statistics

function benchmark(interpType::Type, numDims::Int=6, pointsPerDim::Int=15, numRandomTests::Int=5000; quiet=false, k=1)
    # Calculate interpolation benchmarks
    # create a multi-dimensional test Grid with data, cut points, etc.

    # Set up the data structures:
    cutPointMat = 0:(1/(pointsPerDim-1)):1
    cutPoints = [cutPointMat for i=1:numDims]
    cut_counts = Int[length(cutPoints[i]) for i = 1:length(cutPoints)]

    if interpType == SimplexGrid
        grid = SimplexGrid(tuple(cutPoints...)...)
    elseif interpType == RectangleGrid
        grid = RectangleGrid(tuple(cutPoints...)...)
    end

    data = rand(length(grid))

    # Warmup
    interpolate(grid,data,rand(numDims))

    # Run the evaluation: totalInterpolations attempts
    # at interpolating numRandomTests points
    totalInterpolations = 1000

    timesToInterpolate = zeros(numRandomTests)

    averageTime = 0
    for test=1:numRandomTests
    testPoint = rand(numDims, totalInterpolations)
    elapsedTime = @elapsed begin
        for i=1:totalInterpolations
            interpolate(grid,data,testPoint[:,i])
        end
    end
    timesToInterpolate[test] = elapsedTime
    end
    averageTime = mean(timesToInterpolate)
    plusMinus95CI = std(timesToInterpolate)*2
    #@show numDims,pointsPerDim, length(grid), averageTime

    if !quiet
        println("$(label(grid)) interpolation of $numDims dimensions with $pointsPerDim cut points per dimension:")
        println("  requires $averageTime seconds (+/- $plusMinus95CI for 95% CI) on average to interpolate $totalInterpolations test points (averaged over $numRandomTests attempts)")
    end

    return averageTime, plusMinus95CI
end

function compareBenchmarks(numDims::Int=6, pointsPerDim::Int=15, numRandomTests::Int=1000; quiet=false)
    # Compare rectangular and simplex interpolation benchmarks

    elapsedTimeSimplex, sdSimplex = benchmark(SimplexGrid, numDims, pointsPerDim, numRandomTests, quiet=true);
    elapsedTimeRectangle, sdRectangle = benchmark(RectangleGrid, numDims, pointsPerDim, numRandomTests, quiet=true);

    println("$numRandomTests interpolations of $numDims dimensions with $pointsPerDim cut points per dimension:")
    println("  Rectangle required $elapsedTimeRectangle +/- $sdRectangle sec")
    println("  Simplex   required $elapsedTimeSimplex +/- $sdSimplex sec")

end

function compareSpeedUp(nDims, nPoints; marginoferror=0.1)
    # Holds the number of nPoints constant to within
    # marginoferror, varying the number of dimensions,
    # and returns the speedup of simplex over the
    # rectangle interpolation

    nTests = 30

    # Compile everything
    benchmark(RectangleGrid, 3, 10, nTests, quiet=true);

    speedup = []
    println("How large is the simplex grid speed up over the multilinear grid?")
    for i = 2:nDims
        nPointsPerDim = nPoints^(1/i)
        nPointsPerDim = convert(Int, round(nPointsPerDim))

        if abs(nPointsPerDim^i - nPoints) > marginoferror*nPoints
            continue
        end

        sspeed, ssd = benchmark(SimplexGrid, i, nPointsPerDim, nTests, quiet=true);
        println("  limiting to $i dimensions and therefore $nPointsPerDim points per dim:")
        println("    mean speed: $sspeed, std dev: $ssd")

        #push!(speedup, (i,rspeed/sspeed) )
        GC.gc()
    end
    return speedup
end

compareBenchmarks();
compareSpeedUp(1000,100000)

# Warm the cache & get results
#benchmark(RectangleGrid, quiet=true)
#benchmark(RectangleGrid, quiet=true)
#benchmark(RectangleGrid)

# Warm the cache & get results
#benchmark(SimplexGrid, quiet=true)
#benchmark(SimplexGrid, quiet=true)
#benchmark(SimplexGrid)

#=
Profile.clear()
@profile (for i = 1:100; benchmark(SimplexGrid, k=1); end)
Profile.print()
Profile.print(format=:flat)
=#


text:
using LinearAlgebra, DataFrames,Optim, ForwardDiff, BenchmarkTools,Distributions,Expectations, QuantEcon, Statistics
using Distributions: invsqrt2π, log2π, sqrt2, invsqrt2
using DCDC
using Test
using Distributed, Suppressor

# cd(joinpath(DEPOT_PATH[1],"dev","DCDC"))
begin
    n = 300
    xdata = randn(n,2);
    x = [0.3,0.4]
    h = [2,2];
    y = xdata * [1,2] + randn(n);
    w = zeros(n);
    ekernel4(x,xdata,h,w,n); #Assign value to weight
    w_diag = diagm(0=>w);
    β_kernel = inv(xdata'*w_diag * xdata) * (xdata' * w_diag * y)
    β_OLS = inv(xdata'*xdata) * (xdata'*y);

    af = ApproxFn(xdata,y,:gaussian,2);

end
# Multi-Variable Kernel Regression
@testset "Kernel Regressions" begin
    @test β_kernel != β_OLS;
    @test any(abs.(β_OLS - [1,2]) .< 1)
end

@testset "Kernel bandwidth constant" begin
    @test abs(bw_constant(1,2,:epan) - 2.34) < 0.01;
    @test abs(bw_constant(1,4,:epan) - 3.03) < 0.01;
    @test abs(bw_constant(1,6,:epan) - 3.53) < 0.01;
    @test abs(bw_constant(1,4,:triw) - 3.72) < 0.01;
    @test abs(bw_constant(1,4,:gaussian) - 1.08) < 0.01;
end

param1 = Parameter();
param2 = Parameter([1,2],[1,2]);
@testset "Parameter test" begin
    @test param1.γ == [1]
    @test param2.γ == [1,2]
end
@testset "State" begin
    state = State()
    state = State(rand(),randn(3),rand())
    @test typeof(state) == State;
end
@testset "ProfitFn" begin
    profitfn = ProfitFn(param2);
    @test typeof(profitfn) == ProfitFn;
    @test profitfn() == 0;
    pf = ProfitFn(param1);
    @test pf() == 0;
end

# Plot
using Plots

@testset "ApproxFn" begin
    n=300;
    xdata=4*randn(n);
    y=sin.(xdata)+ 0.3*randn(n);
    scatter(xdata,y)
    af=ApproxFn(xdata,y,:gaussian,2);
    y_hat= af(xdata);
    scatter!(xdata,y_hat);
end

# DDC
@testset "GMM small test" begin
      n = 3000;
      sig = rand(3,3);
      mu = rand(3);
      mn = MvNormal(mu,sig * sig');
      x_m = rand(mn,n)';
      x = x_m[:,[1,2]];
      β₀ = [1,3];
      y = x * β₀ + x_m[:,3];
      @show mn.Σ.mat;
      @show cov(x);
      β_OLS  = inv(x'*x) *(x'*y);
      obj = θ->(begin
            η = (x * θ - y) ;
            η'*η
            end)
      result = optimize(obj,zeros(2),Newton(),autodiff = :forward)
      @test result.minimizer ≈ β_OLS ;
end

@testset "Test Optim" begin
      n = 3000;
      sig = [01.0 0.5 0.5;
             0.0 0.5 0.0 ;
             0.0 0.0 1.0];

      mu = [1,9,0];
      mn = MvNormal(mu,sig * sig')
      W = rand(mn,n)'
      x = W[:,1];z=W[:,2];ϵ=W[:,3];
      y = x * 3 + ϵ;
      β_OLS = inv(x'*x)*(x'*y);
      β_IV = inv(z'*x)*(z'*y);
      w = hcat(x,y,z);
      obj = θ->(begin
            η = (w[:,1] .* θ - w[:,2]) .* w[:,3];
            η'*η
            end)
      result = optimize(obj,zeros(1),Newton(),autodiff = :forward)
      # @show converged(result) || error("Failed to converge in $(iterations(result)) iterations")
      @show xmin = result.minimizer
      @show result.minimum
      @test (β_OLS - result.minimizer[1] )< 0.1
end

@testset "Dynamic Decision Process" begin
    σ₀ = 1;
    β = 0.8;
    nM = 50;
    nT = 5;
    ddc = DynamicDecisionProcess(σ₀,0.8);
    UpdateVal!(ddc);
    plot(ddc.ValueFn.xdata,ddc.ValueFn.y);
    plot(ddc.PolicyFn.xdata,ddc.PolicyFn.y);
    data = simulate_ddc(nM,nT,ddc);
end

# Generate IV regression data
@testset "IV regression" begin
    n = 3000;
    sig = [1.0 0.5 0.5;
         0.0 0.5 0.0 ;
         0.0 0.0 1.0];

    mu = [1,9,0];
    mn = MvNormal(mu,sig * sig')
    W = rand(mn,n)'
    x = W[:,1];z=W[:,2];ϵ=W[:,3];
    y = x * 3 + ϵ;
    β_OLS = inv(x'*x)*(x'*y);
    β_IV = inv(z'*x)*(z'*y);
    w = hcat(x,y,z);
    @test β_OLS != β_IV;
end

@testset begin "Solve equilibrium"
    ddc=DynamicDecisionProcess(2,0.8);
    computeEquilibrium(ddc);
    checked = [check_ee(ddc) for i = 1:100];
end


text:
# Copyright 2014 Mitchell Kember. Subject to the MIT License.
# Google Code Jam Round 1A 2008
# Problem C: Numbers

# Returns the last three digits before the decimal point of the sum 3 + sqrt(5)
# raised to the nth power, as a 3-character string left-padded with zeros.
function last3(n)
	mat = [0 1; -4 6]^(n-1)
	x = (6mat[1,1] + 28mat[1,2] + 1_000_000) % 1000
	x = (x == 0) ? 999 : x - 1
	lpad(x, 3, 0)
end

function main()
	n = int(readline())
	for i = 1:n
		x = int(readline())
		println("Case #$i: ", last3(x))
	end
end

main()


text:
function primefactors(x)
    factors = []
    i=2
    while x != i
        if x%i==0
            push!(factors, i)
            x = x/i
            i = 2
        else
            i += 1
        end
    end
    push!(factors, i)
    return factors
end
print("Number to factor: ")
num = readline()
num = parse(Int64, num)
display(primefactors(num))

max_stars_repo_name: ellygaytor/Projects

max_issues_repo_name: ellygaytor/Projects

max_forks_repo_name: ellygaytor/Projects


text:
module ImagineProcedures

using Unitful, Reexport

using ImagineInterface

@reexport using CachedCalls
import CachedCalls: CachedCall, description

export SignalGenerator,
        ImagineProcedure,
        outputs,
        process

#This is just a CachedCall with a Vector{ImagineSignal} return type
const SignalGenerator{F, A} = CachedCall{F, A, Vector{ImagineSignal}}
SignalGenerator(description::String, f::F, args::A) where {F<:Function, A} = SignalGenerator{F, A}(description, f, args, nothing)
SignalGenerator(f::F, args::A) where {F<:Function, A} = SignalGenerator("", f, args)

#An ImagineProcedure couples a function for generating command signals with a function for analyzing those commands in the context of something else.
#Typically "something else" is one or more results of running those commands on the microscope (analog and digital signal recordings, camera images).
#This design was chosen because currently Imagine can't be operated directly from Julia, but we would still like to store a repeatable set of procedures
#(for example testing and calibration procedures).  It _may_ also prove useful for standardizing storage of scientific experiments and their analyses.
#user-facing api is outputs() for getting a set of output commands and process(args...) for running its analyze function with the outputs as first argument + set of other args
#therefore when creating an ImagineProcedure it is critical that the analyze() function expects a set of output signals as its first argument
struct ImagineProcedure
    description::String
    sig_gen::SignalGenerator #SignalGenerator generates signals to be run with Imagine
    analyze::Function #a function of call!(sig_gen) and one or more args
end

#getters. probably wont' export?
description(ip::ImagineProcedure) = ip.description
signal_generator(ip::ImagineProcedure) = ip.sig_gen
analyze_function(ip::ImagineProcedure) = ip.analyze

outputs(ip::ImagineProcedure) = call!(signal_generator(ip))
process(ip::ImagineProcedure, inputs...) = analyze_function(ip)(outputs(ip), inputs...)

function Base.show(s::IO, p::ImagineProcedure)
    write(s, "ImagineProcedure")
    write(s, "\nDescription:                    $(description(p))")
    write(s, "\nSignalGenerator description:    $(description(signal_generator(p)))")
    write(s, "\nAnalysis function:              $(analyze_function(p))")
end

end # module


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# recipe for base case
@recipe function f(segment::Segment)
  seriestype --> :scatterpath
  seriescolor --> :black
  primary --> false

  Tuple.(coordinates.(vertices(segment)))
end

# plot facets in a recursion
@recipe function f(polytope::Polytope)
  for f in facets(polytope)
    @series begin
      primary --> false
      f
    end
  end
end

@recipe function f(ray::Ray)
  seriestype --> :path
  seriescolor --> :black
  arrow --> true

  a = coordinates(ray(0))
  b = coordinates(ray(1))
  [Tuple(a), Tuple(b)]
end

@recipe function f(sphere::Sphere, nsamples=100)
  seriestype --> :path
  seriescolor --> :black

  points = collect(sample(sphere, RegularSampling(nsamples)))
  [points; first(points)]
end

@recipe function f(chain::Chain)
  seriestype --> :path
  seriescolor --> :black
  label --> "chain"

  xs = coordinates.(vertices(chain))

  isclosed(chain) && push!(xs, xs[begin])

  Tuple.(xs)
end

@recipe function f(polyarea::PolyArea)
  seriestype --> :path
  seriescolor --> :black
  label --> "polyarea"

  pchains = chains(polyarea)

  # plot outer chain
  @series begin
    first(pchains)
  end

  # plot inner chains
  for pchain in pchains[2:end]
    @series begin
      primary --> false
      pchain
    end
  end
end

@recipe function f(geometries::AbstractVector{<:Geometry})
  for geometry in geometries
    @series begin
      geometry
    end
  end
end


text:
module DyldExport

using MachO
using StrPack

# Constants we'll need
import MachO: N_EXT, N_ABS, NO_SECT,
    MH_MAGIC_64, CPU_TYPE_X86_64, MH_DYLIB, MH_NOUNDEFS,
    LC_SYMTAB, LC_DYSYMTAB, LC_SEGMENT_64, MH_BUNDLE,
    MH_PREBOUND, MH_EXECUTE, MH_DYLIB_STUB

# And datastructures we'll need
import MachO: mach_header_64, nlist_64, symtab_command, dysymtab_command,
    segment_command_64

if OS_NAME != :Darwin
    error("Currently only supported on OS X")
end

function create_sections(symbols::Dict{Symbol,Ptr{Void}})
    symtab = IOBuffer()
    strtab = IOBuffer()
    symbols = [("_$k",v) for (k,v) in symbols]
    sort!(symbols,by=x->x[1])
    for (name,value) in symbols
        sym = nlist_64(
            position(strtab), # n_strx
            N_EXT | N_ABS,    # n_type
            NO_SECT,          # n_sect
            0,                # n_desc
            value)            # n_value
        pack(symtab, sym)
        write(strtab, name)
        write(strtab, UInt8('\0'))
    end
    symtab, strtab
end

function compute_size(symbols)
    headersize = sizeof(mach_header_64)
    lcssize = 8*sizeof(UInt32)+sizeof(symtab_command)+sizeof(dysymtab_command)+
        2*sizeof(segment_command_64)
    symsize = 0
    for (k,v) in symbols
        symsize += sizeof(nlist_64)
        symsize += sizeof(string(k))+1
    end
    max(headersize+lcssize+symsize,4096)
end

function create_image(symbols::Dict{Symbol,Ptr{Void}},vmbase)
    totalsize = compute_size(symbols)
    buf = IOBuffer()
    symtab, strtab = create_sections(symbols)

    # Write dummy versions of everything so we know
    # the offsets
    lcssize = 8*sizeof(UInt32)+sizeof(symtab_command)+sizeof(dysymtab_command)+
        2*sizeof(segment_command_64)
    pack(buf, mach_header_64(
        MH_MAGIC_64,     # magic
        CPU_TYPE_X86_64, # cputype
        0,               # cpusubtype
        MH_DYLIB,        # filetype
        4,               # ncmds
        # sizeofcmds
        lcssize,
        MH_NOUNDEFS|MH_PREBOUND,    # flags
    ))

    lcpos = position(buf)

    # Needs a dummy writable segment, which is also named __TEXT and at least
    # one byte large (or dyld crashes - filed as radar 23944790)
    write(buf, UInt32(LC_SEGMENT_64))
    write(buf, UInt32(sizeof(segment_command_64)+2*sizeof(UInt32)))
    pack(buf, segment_command_64(
        MachO.small_fixed_string(reinterpret(UInt128,UInt8['_','_','T','E','X','T',0,0,0,0,0,0,0,0,0,0])[]),
        vmbase,4096,0,4096,0x2,0x2,0,0
    ))

    # Needed, and also needs be be non-zero for dyld not to crash ugh!
    write(buf, UInt32(LC_SEGMENT_64))
    write(buf, UInt32(sizeof(segment_command_64)+2*sizeof(UInt32)))
    pack(buf, segment_command_64(
        MachO.small_fixed_string(reinterpret(UInt128,UInt8['_','_','L','I','N','K','E','D','I','T',0,0,0,0,0,0])[]),
        vmbase+4096,totalsize,0,totalsize,0x1,0x1,0,0
    ))

    write(buf, UInt32(LC_SYMTAB))
    write(buf, UInt32(sizeof(symtab_command)+2*sizeof(UInt32)))
    symtablc_offs = position(buf)
    pack(buf, symtab_command())

    write(buf, UInt32(LC_DYSYMTAB))
    write(buf, UInt32(sizeof(dysymtab_command)+2*sizeof(UInt32)))
    pack(buf, dysymtab_command(
        0,               # ilocalsym
        0,               # nlocalsym
        0,               # iextdefsym
        length(symbols), # nextdefsym
        0,               # iundefsym
        0,               # nundefsym
        0,               # tocoff
        0,               # ntoc
        0,               # modtaboff
        0,               # nmodtab
        0,               # extrefsymoff
        0,               # nextrefsyms
        0,               # indirectsymoff
        0,               # nindirectsyms
        0,               # extreloff
        0,               # nextrel
        0,               # locreloff
        0,               # nlocrel
    ))

    symtab_offs = position(buf)
    write(buf, takebuf_array(symtab))

    strtab_offs = position(buf)
    data = takebuf_array(strtab)
    strtab_size = sizeof(data)
    write(buf, data)

    seek(buf, symtablc_offs)
    pack(buf, symtab_command(
        symtab_offs,
        length(symbols),
        strtab_offs,
        strtab_size
    ))

    seekend(buf)

    # Pad the file
    if (position(buf) < 4096)
        write(buf, zeros(UInt8,4096-position(buf)))
    end

    buf
end

function export_symbols(symbols)
    # Reserve some memory right up until we're ready to dlopen
    x = Ref{Ptr{Void}}()
    x[] = 0
    task = unsafe_load(cglobal(:mach_task_self_,Ptr{Void}))
    @assert task != 0
    size = DyldExport.compute_size(symbols) + 4096
    ret = ccall(:vm_allocate,Cint,(Ptr{Void},Ref{Ptr{Void}},Csize_t,Bool),task,x,size,true)
    @assert ret == 0

    buf = DyldExport.create_image(symbols, x[])
    data = takebuf_array(buf)
    name = tempname()
    open(name,"w") do f
        write(f,data)
    end
    @show x[]
    @assert ccall(:vm_deallocate,Cint,(Ptr{Void},Ptr{Void},Csize_t),
        task,x[],size) == 0
    @show name
    handle = Libdl.dlopen(name,Libdl.RTLD_GLOBAL)
    #rm(name)
    name,handle
end

const NSObjectFileImageFailure = 0
const NSObjectFileImageSuccess = 1
const NSObjectFileImageInappropriateFile = 2
const NSObjectFileImageArch = 3
const NSObjectFileImageFormat = 4
const NSObjectFileImageAccess = 5

function NSCreateObjectFileImageFromMemory(data)
    out = Ref{Ptr{Void}}()
    ret = ccall(:NSCreateObjectFileImageFromMemory,Cint,
        (Ptr{Void},Csize_t,Ref{Ptr{Void}}),pointer(data),sizeof(data),out)
    if ret != NSObjectFileImageSuccess
        error("Failed to open object file")
    end
    out[]
end

# Hack for now
if isdefined(Main,:Cxx)
    include("cxxsupport.jl")
end

end # module


text:
abstract type AbstractBasis{T,I} <: AbstractVector{T} end

struct Basis{T,I,A<:AbstractVector{T}} <: AbstractBasis{T,I}
    basis::A
    rbasis::Dict{T,I}
end

function Basis{I}(basis::AbstractVector) where {I}
    length(basis) <= typemax(I) ||
        throw("index type $I is to small for basis of length $(length(basis))")
    @assert !(eltype(basis) <: Integer)
    return Basis(basis, Dict(b => I(idx) for (idx, b) in pairs(basis)))
end

Base.size(b::Basis) = size(b.basis)

Base.IndexStyle(::Type{<:Basis{T,I,A}}) where {T,I,A} = Base.IndexStyle(A)

Base.@propagate_inbounds Base.getindex(b::Basis, i::Integer) = b.basis[i]
Base.@propagate_inbounds Base.getindex(b::Basis{T}, g::T) where {T} = b.rbasis[g]

Base.in(g, b::Basis) = haskey(b.rbasis, g)

## are these really that useful?
SparseArrays.indtype(::Type{<:Basis{T,I}}) where {T,I} = I
SparseArrays.indtype(b::Basis) = SparseArrays.indtype(typeof(b))


text:
########################################################
function lsststd(; zpwave=3670.7, z=0.5, lum=5e43) # Angstrom, unitless, erg/s
########################################################

    #Define QSO variability-luminosity-redshift relation
    #Morganson et al. (2014) ApJ 784:92
    #input luminosity in erg/s
    #See also Roberto Cid Fernandes, Jr, et al. (1996)
    angstonm = 0.1
    zpwave   = zpwave*angstonm
    sigmavar = 0.079*(1+z)^0.15 * (lum/1e46)^-0.2 * (zpwave/1000.)^-0.44
    return sigmavar
end

########################################################
function setlsstwavelengths(λ = [3670.7, 4826.9, 6223.2, 7546.0, 8690.0, 9710.3]) # Angstrom
########################################################

    lsstwaves() =  λ

end

########################################################
function unitgalaxyvector(T)
########################################################

    G = galaxyvector(T)

    G / norm(G)

end

########################################################
function galaxyvector(T)
########################################################

    G = map(λ -> B(; T=T, λ=λ), lsstwaves())

end

########################################################
function B(; λ=λ, T=T)   # λ is in Angstrom, T in Kelvin
########################################################

    h = 6.6261e-27      # Planck constant    [cgs]
    c = 2.9979245800e10 # Speed of light     [cm/s]
    k = 1.380648813e-16 # Boltzmann constant [cgs]

    wavcm = λ*1e-8      # convert Ang to cm
    nu    = c/wavcm     # convert wavelength to frequency [Hz]

    #--frequency units--
    # emittence [erg/sec/cm^2/Hz]

    In   = 2.0*pi*h*nu^3.0/c^2 * (1.0/(exp(h*nu/(k*T)) - 1.0))

    ImJy = In / 1.0e-26

    return ImJy # Flux in absolute units of milli jansky
                # Doesn't not contain observational effects
end


text:
using VortexDistributions, FourierGPE, Plots, LaTeXStrings
using SpecialFunctions

## Initialize simulation
L = (18.0,18.0)
N = (256,256)
sim = Sim(L,N)
@unpack_Sim sim

## set simulation parameters
μ = 12.0
tf = 10.0
t = LinRange(ti,tf,Nt)

# Time dependent potential function (here trivial t dep)
import FourierGPE.V
V(x,y,t) = 0.5*(x^2 + y^2)
ψ0(x,y,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,0.0)/μ,0.0)+im*0.0)

## make initial state
x,y = X
ψi = ψ0.(x,y',μ,g)
ϕi = kspace(ψi,sim)
@pack_Sim! sim

## evolve
@time sol = runsim(sim)

## plot ground state
ϕg = sol[end]
ψg = xspace(ϕg,sim)
healinglength(x,y,μ,g) = 1/sqrt(g*abs2(ψ0(x,y,μ,g)))
R(μ,g) = sqrt(2*μ)
showpsi(x,y,ψg)
Rtf = R(μ,g)
n0 = maximum(abs2.(ψg))

## spectra of ground state
ξ = healinglength(0,0,μ,g)
kR = 2*pi/Rtf
kξ = 2*pi

kmin = 0.1*kR
kmax = kξ

Np = 300
k = log10range(kmin,kmax,Np)

@time Ek = kespectrum(k,ψg,X,K)
plot(k*ξ,Ek,scale=:log10)

# analytic
f(x) = (sin(x)-x*cos(x))^2/x^3
Ektf(k,R) = π*n0*R*f.(k*R)

# linear scale
plot(k*ξ,Ek)
plot!(k*ξ,Ektf.(k,Rtf), c=:red)

## spectrum
kR = 2*pi/Rtf
kξ = 2*pi

kmin = 0.1*kR
kmax = 1.3kξ
Np = 200
k = log10range(kmin,kmax,Np)

Eki = ikespectrum(k,ψg,X,K)
Ekc = ckespectrum(k,ψg,X,K)
Ekq = qpespectrum(k,ψg,X,K)
plot(k*ξ,Eki,label=L"E_k^i(k)")
plot!(k*ξ,Ekc,label=L"E_k^c(k)")
plot!(k*ξ,Ekq,label=L"E_k^q(k)")

plot(k*ξ,Ekq,scale=:log10)
vline!([kR*ξ],ls=:dash)
vline!([1.0],ls=:dash)

## also plot interaction and trap energies without vortex
using QuadGK
# trap energy
int1(k)=quadgk(x-> x^2*sqrt(1-x^2)*besselj0(k*x),0,1)[1]

# interaction energy
int2(k)=quadgk(x-> x*(1-x^2)*besselj0(k*x),0,1)[1]

plot(k,)
## make central vortex
ψd = copy(ψg)
ξv = healinglength(0,0.,μ,g)
# can use default ξ = 1
xp,yp = 0.,d/2
xn,yn = 0.,-d/2
pv = PointVortex(xp,yp,1)
nv = PointVortex(xn,yn,-1)
dipole = ScalarVortex(ξv,[nv;pv])

psi = Torus(ψd,x,y) # set field topology for VortexDistributions
vortex!(psi,dipole) # make dipole
showpsi(x,y,ψd)

## evolve
tf = 1.0
t = LinRange(ti,tf,Nt)
γ = 0.
ψi .= psi.ψ
ϕi = kspace(ψi,sim)

@pack_Sim! sim
@time sol = runsim(sim)

## plot final state
ϕf = sol[end]
ψf = xspace(ϕf,sim)
showpsi(x,y,ψf)

## analyse
kL = 2*pi/L[1]
kξ = 2*pi

kmin = 0.1*kL
kmax = kξ

Np = 200
k = log10range(kmin,kmax,Np)

@time Ek = kespectrum(k,ψf,X,K)
plot(k,Ek,scale=:log10)

# heatmap(log.(abs2.(A |> fftshift) .+ eps.()))

## incompressible spectrum
kL = 2*pi/L[1]
kξ = 2*pi
kd = 2*pi/d

kmin = 0.1*kL
kmax = 1.3kξ
Np = 400
k = log10range(kmin,kmax,Np)

@time Eki = ikespectrum(k,ψf,X,K)

## power-law plot in logspace
n0 = abs2.(ψf[1,1])
E0 = π*n0
p1 = plot(k,Eki/E0,scale=:log10,grid=false,label=L"E_i(k)",legend=:bottomleft)
# plot!(k,200k.^(-3),label=L"k^{-3}")
# plot!(k,2000*k,label=L"k")
# ylims!(2e-3,1e2)
# xlims!(0.02,10)
vline!([kL],ls=:dash,label=L"k_L")
vline!([kd],ls=:dash,label=L"k_d")
vline!([kξ],ls=:dash,label=L"k_\xi")
xlabel!(L"k\xi")
ylabel!(L"E_i(k)")

# test against analytic form
# use the numerical dipole distance
# vort = findvortices(Torus(ψf,x,y)) |> rawData
# d = vort[:,2] |> diff
d = 2
Λ = 0.8249
f(x) = x*(besselk(1,x)*besseli(0,x)-besselk(0,x)*besseli(1,x))
F(x,Λ) = f(x/(2Λ))^2/x
F(x) = F(x,Λ)
Ed(k,d) = 2*F(k)*(1-besselj0(k*d))
plot!(k,Ed.(k,[d]),label=L"{\cal E}_i^a(k)")


text:
export Derivation, NullDerivation, BasicDerivation, ExtensionDerivation,
   CoefficientLiftingDerivation, TowerOfDifferentialFields,
   BaseDerivation, MonomialDerivative, domain, constant_field,
   isbasic, isprimitive, ishyperexponential, isnonlinear, ishypertangent,
   iscompatible, isnormal, isspecial, issimple, isreduced, is_Sirr1_eq_Sirr

abstract type Derivation end

domain(D::Derivation) = D.domain
iscompatible(p::RingElement, D::Derivation) = parent(p)==domain(D)
iscompatible(f::FracElem, D::Derivation) =
    base_ring(parent(f))==D.domain

is_Sirr1_eq_Sirr(D::Derivation) = true
# For the time being it is assumed that S₁ⁱʳʳ==Sⁱʳʳ for all derivations
# to be considered.
# See Bronstein's book, Theorems 5.1.1, 5.1.2, 5.10.1 .

struct NullDerivation <: Derivation
    domain #::Ring does not work ... :(
end

NullDerivation(domain::FracField{T}) where T<:RingElement = NullDerivation(base_ring(domain))

function (D::NullDerivation)(p::RingElement)
    parent(p)==D.domain || error("p not in domain of D")
    zero(p)
end

function (D::NullDerivation)(f::FracElem{T}) where {T<:RingElement}
    base_ring(parent(f))==D.domain || error("f not in domain of D")
    zero(f)
end

function constant_field(D::NullDerivation)
    R = D.domain
    if isa(R, AbstractAlgebra.Field)
        return R
    else
        return FractionField(R)
    end
end

Base.show(io::IO, D::NullDerivation) = print(io, "Null derivation D=0 on ", domain(D))

struct BasicDerivation{T<:RingElement} <: Derivation
    domain::PolyRing{T}
end

BasicDerivation(domain::FracField{P}) where P<:PolyElem = BasicDerivation(base_ring(domain))

function (D::BasicDerivation{T})(p::PolyElem{T}) where T<:RingElement
    iscompatible(p, D) || error("p not in domain of D")
    derivative(p)
end

function (D::BasicDerivation{T})(f::FracElem{P}) where {T<:RingElement, P<:PolyElem{T}}
    iscompatible(f, D) || error("f not in domain of D")
    derivative(f)
end

BaseDerivation(D::BasicDerivation) = NullDerivation(base_ring(D.domain))
MonomialDerivative(D::BasicDerivation) = one(D.domain)
function constant_field(D::BasicDerivation)
    R = base_ring(D.domain)
    if isa(R, AbstractAlgebra.Field)
        return R
    else
        return FractionField(R)
    end
end

Base.show(io::IO, D::BasicDerivation) = print(io, "Basic derivation D=d/d", gen(domain(D)), " on ", domain(D))

struct ExtensionDerivation{T<:RingElement} <: Derivation
    domain::PolyRing{T}
    D::Derivation
    H::PolyElem{T}
    function ExtensionDerivation(domain::PolyRing{R}, D::Derivation, H::PolyElem{R}) where R<:RingElement
        base_ring(domain)==D.domain || error("base ring of domain must be domain of D")
        new{R}(domain, D, H)
    end

    function ExtensionDerivation(domain::PolyRing{F}, D::Derivation, H::PolyElem{F}) where
        {R<:RingElement, F<:FracElem{R}}
        base_ring(base_ring(domain))==D.domain || error("base ring of domain must be domain of D")
        new{F}(domain, D, H)
    end
end

function ExtensionDerivation(domain::FracField{P}, D::Derivation, H::P) where {T<:RingElement, P<:PolyElem{T}}
    ExtensionDerivation(base_ring(domain), D, H)
end

function CoefficientLiftingDerivation(domain::PolyRing{T}, D::Derivation) where T<:RingElement
    ExtensionDerivation(domain, D, zero(domain))
end

function CoefficientLiftingDerivation(domain::FracField{T}, D::Derivation) where T<:RingElement
    ExtensionDerivation(base_ring(domain), D, zero(base_ring(domain)))
end

function (D::ExtensionDerivation{T})(p::PolyElem{T}) where T<:RingElement
    iscompatible(p, D) || error("p not in domain of D")
    if iszero(D.H)
        return map_coefficients(c->D.D(c), p)
    else
        return map_coefficients(c->D.D(c), p) + D.H*derivative(p)
    end
end

function (D::ExtensionDerivation{T})(f::FracElem{P}) where {T<:RingElement, P<:PolyElem{T}}
    iscompatible(f, D) || error("f not in domain of D")
    a = numerator(f)
    b = denominator(f)
    if isone(b)
        if iszero(D.H)
            return map_coefficients(c->D.D(c), a) + zero(f)
        else
            return map_coefficients(c->D.D(c), a) + D.H*derivative(a) + zero(f)
        end
    else
        if iszero(D.H)
            da = map_coefficients(c->D.D(c), a)
            db = map_coefficients(c->D.D(c), b)
        else
            da = map_coefficients(c->D.D(c), a) + D.H*derivative(a)
            db = map_coefficients(c->D.D(c), b) + D.H*derivative(b)
        end
        return (b*da - a*db)//b^2
    end
end

BaseDerivation(D::ExtensionDerivation) = D.D
MonomialDerivative(D::ExtensionDerivation) = D.H
constant_field(D::ExtensionDerivation) = constant_field(D.D)

Base.show(io::IO, D::ExtensionDerivation) = print(io, "Extension by D",
    gen(domain(D))," = ", MonomialDerivative(D),
    " of ", BaseDerivation(D), " on ", domain(D))

struct AlgebraicExtensionDerivation{T<:FieldElement, P<:PolyElem{T}} <: Derivation
    domain::AbstractAlgebra.ResField{P}
    D::Derivation
    dy::AbstractAlgebra.ResFieldElem{P}
    function AlgebraicExtensionDerivation(domain::AbstractAlgebra.ResField{P}, D::Derivation) where {T<:FieldElement, P<:PolyElem{T}}
        base_ring(base_ring(base_ring(domain)))==D.domain || error("base ring of domain must be domain of D")
        p = modulus(domain)
        y = domain(gen(base_ring(domain)))
        dy = -map_coefficients(derivative, p)(y)//derivative(p)(y)
        new{T,P}(domain, D, dy)
    end
end

function (D::AlgebraicExtensionDerivation)(f::AbstractAlgebra.ResFieldElem{P}) where {T<:FieldElement, P<:PolyElem{T}}
    iscompatible(f, D) || error("f not in domain of D")
    map_coefficients(derivative, data(f)) + derivative(data(f))*D.dy
end

BaseDerivation(D::AlgebraicExtensionDerivation) = D.D
#Note: implementation of constant_field requires further thought ...
#constant_field(D::AlgebraicExtensionDerivation) = constant_field(D.D)

Base.show(io::IO, D::AlgebraicExtensionDerivation) = print(io, "Algebraic extension of ",
    BaseDerivation(D), " on ", domain(D))

AbstractAlgebra.degree(D::Derivation) =
    degree(MonomialDerivative(D)) # \delta(t), see Def. 3.4.1
AbstractAlgebra.leading_coefficient(D::Derivation) =
leading_coefficient(MonomialDerivative(D)) # \lambda(t), see Def. 3.4.1
Base.iszero(D::Derivation) = false
Base.iszero(D::NullDerivation) = true
isbasic(D::Derivation) = false
isbasic(D::BasicDerivation) = true
isbasic(D::ExtensionDerivation) =
    iszero(BaseDerivation(D)) && isone(MonomialDerivative(D))
isprimitive(D::Derivation) = degree(D)==0 # see Def. 5.1.1
ishyperexponential(D::Derivation) = # see Def, 5.1.1
    degree(D)==1 && iszero(constant_coefficient(MonomialDerivative(D)))
isnonlinear(D::Derivation) = degree(D)>=2 # see Def. 3.4.1

function ishypertangent(D::Derivation)
    # see Def. 5.10.1
    isnonlinear(D) || return false
    t = gen(domain(D))
    q, r = divrem(MonomialDerivative(D), t^2+1)
    iszero(r) && degree(q)<=0
end

isnormal(p::PolyElem, D::Derivation) =
    # see Def. 3.4.2
    iscompatible(p, D) && degree(gcd(p, D(p)))==0

isspecial(p::PolyElem, D::Derivation) =
    # see Def. 3.4.2
    iscompatible(p, D) && iszero(rem(D(p), p)) # gcd(p, D(p))==p

issimple(f::FracElem{P}, D::Derivation) where P<:PolyElem =
    #see Def. 3.5.2.
    iscompatible(f, D) && isnormal(denominator(f), D)

isreduced(f::FracElem{P}, D::Derivation) where P<:PolyElem =
    #see Def. 3.5.2.
    iscompatible(f, D) && isspecial(denominator(f), D)

function isconstant(x::T, D::NullDerivation) where T<:RingElement
    @assert iscompatible(x, D)
    true
end

function isconstant(x::F, D::NullDerivation) where {P<:PolyElem, F<:FracElem{P}}
    #this version for disambiguation
    @assert iscompatible(x, D)
    true
end

function isconstant(x::T, D::Derivation) where T<:RingElement
    @assert iscompatible(x, D)
    false
end

function isconstant(x::P, D::BasicDerivation) where P<:PolyElem
    @assert iscompatible(x, D)
    degree(x)<=0
end

function isconstant(x::P, D::ExtensionDerivation) where P<:PolyElem
    @assert iscompatible(x, D)
    if degree(x)>0
        return false
    else
        return isconstant(constant_coefficient(x), BaseDerivation(D))
    end
end

function isconstant(x::F, D::Derivation) where {P<:PolyElem, F<:FracElem{P}}
    @assert iscompatible(x, D)
    isone(denominator(x)) && isconstant(numerator(x), D)
end

function constantize(x::T, D::NullDerivation) where T<:RingElement
    @assert iscompatible(x, D)
    x
end

function constantize(x::F, D::NullDerivation) where {P<:PolyElem, F<:FracElem{P}}
    #this version for disambiguation
    @assert iscompatible(x, D)
    x
end

function constantize(x::T, D::Derivation) where T<:RingElement
    @assert iscompatible(x, D)
    error("not a constant")
end

function constantize(x::P, D::BasicDerivation) where P<:PolyElem
    @assert iscompatible(x, D)
    degree(x)<=0 || error("not a constant")
    constant_coefficient(x)
end

function constantize(x::P, D::ExtensionDerivation) where P<:PolyElem
    @assert iscompatible(x, D)
    degree(x)<=0 || error("not a constant")
    constantize(constant_coefficient(x), BaseDerivation(D))
end

function constantize(x::F, D::Derivation) where {P<:PolyElem, F<:FracElem{P}}
    @assert iscompatible(x, D)
    isone(denominator(x)) || error("not a constant")
    constantize(numerator(x), D)
end

function constant_roots(f::PolyElem{T}, D::Derivation; useQQBar::Bool=false) where T<:FieldElement
    @assert iscompatible(f, D)
    p = map_coefficients(c->constantize(c, BaseDerivation(D)), constant_factors(f))
    if useQQBar
        return roots(p, QQBar)
    else
        return roots(p)
    end
end

function constant_roots(f::PolyElem{T}, D::Derivation; useQQBar::Bool=false) where
    {T<:AbstractAlgebra.ResFieldElem}
    @assert iscompatible(f, D)
    p = map_coefficients(c->constantize(c, BaseDerivation(D)), constant_factors(f))
    pp = map_coefficients(c->real(c), p*map_coefficients(c->conj(c), p))
    g = gcd(pp, derivative(pp))
    if useQQBar
        return roots(g, QQBar)
    else
        return roots(g)
    end
end

"""
    SplitFactor(p, D) -> (pₙ, pₛ)

Splitting factorization.

Given a field `k`, a derivation `D` on `k[t]` and `p` in `k[t]`, return
`pₙ`, `pₛ` in `k[t]` such that `p=pₙ*pₛ`, `pₛ` is special, and each squarefree
factor of `pₙ` is normal.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 3.5, p. 100.
"""
function SplitFactor(p::PolyElem{T}, D::Derivation) where T<:FieldElement
    iscompatible(p, D) || error("polynomial p must be in the domain of derivation D")
    S = divexact(gcd(p, D(p)), gcd(p, derivative(p)))
    if degree(S)==0
        return(p, one(p))
    end
    (qn, qs) = SplitFactor(divexact(p, S), D)
    qn, S*qs
end

"""
    SplitSquarefreeFactor(p, D) -> (Ns, Ss)

Splitting squarefree factorization.

Given a field `k`, a derivation `D` on `k[t]` and `p` in `k[t]`, return
`Ns=[N₁,...,Nₘ]`, `Ss=[S₁,...,Sₘ]` with  `Nᵢ`, `Sᵢ` in `k[t]` such that
`p=(N₁*N₂²*...*Nₘᵐ)*(S₁*S₂²*...*Sₘᵐ)` is a splitting factorization of `p`
and the `Nᵢ` and `Sᵢ` are squarefree and coprime.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 3.5, p. 102.
"""
function SplitSquarefreeFactor(p::PolyElem{T}, D::Derivation) where T<:FieldElement
    iscompatible(p, D) || error("polynomial p must be in the domain of derivation D")
    ps = Squarefree(p)
    Ss = [gcd(ps[i], D(ps[i])) for i=1:length(ps)]
    Ns = [divexact(ps[i], Ss[i]) for i=1:length(ps)]
    return Ns, Ss
end

"""
CanonicalRepresentation(f, D) -> (fₚ, fₛ, fₙ)

Canonical representation.

Given a field `k`, a derivation `D` on `k[t]` and `f` in `k(t)`, return
fₚ, fₛ, fₙ in `k(t)` such that `f=fₚ+fₛ+fₙ` is the canonical representation
of `f`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 3.5, p. 103.
"""
function CanonicalRepresentation(f::FracElem{P}, D::Derivation) where {T<:FieldElement, P<:PolyElem{T}}
    # See Bronstein's book, Section 3.5, p. 103
    iscompatible(f, D) || error("rational function f must be in the domain of derivation D")
    a = numerator(f)
    d = denominator(f)
    q, r = divrem(a, d)
    dn, ds = SplitFactor(d, D)
    b, c = gcdx(dn, ds, r)
    q, b//ds, c//dn
end


text:
@symbol_func function epsilon_P(cur_reactor::AbstractReactor)
  cur_epsilon = a_P(cur_reactor)

  cur_epsilon /= safe_symbol(cur_reactor, :R_0)

  cur_epsilon
end


text:
using JSON
using JWTs

function main()
    largs = length(ARGS)

    if (largs < 2) || (largs > 3)
        println("Usage: julia jwt_create.jl <claims file> <keys> [key id]")
        println("    claims file: text file with JSON claim to be encoded")
        println("    keys: URL/file listing keys in standard OIDC key listing format")
        println("    key id: key id to use from the keys listed (optional, first key in list if not specified)")
        exit(1)
    end

    keysurl = ARGS[2]
    if !(startswith(keysurl, "http://") || startswith(keysurl, "https://") || startswith(keysurl, "file://")) && isfile(keysurl)
        keysurl = "file://" * keysurl
    end
    @info("using keys from " * keysurl)
    keyset = JWKSet(keysurl)
    refresh!(keyset)
    kid = (largs == 3) ? ARGS[3] : first(keys(keyset.keys))
    @info("using key id " * kid)

    jwt = JWT(; payload=JSON.parse(readchomp(ARGS[1])))
    sign!(jwt, keyset, kid)
    println(string(jwt))
end

main()


text:
erode_ref(img, dims, r) = erode_ref(img, strel_box(img, dims; r))
erode_ref(img, se) = extreme_filter(min, img, se)

@testset "erode" begin
    test_types = [Bool, Int, Float64, Gray{N0f8}, Gray{Float64}]
    for T in test_types
        for N in (1, 2, 3)
            sz = ntuple(_ -> 32, N)
            img = rand(T, sz...)

            out = erode(img)
            @test out == erode(img, strel_box(img, ntuple(identity, N); r=1))
            @test out == erode_ref(img, ntuple(identity, N), 1)

            @test erode(img; dims=(1,), r=2) == erode_ref(img, (1,), 2)

            se = centered(rand(Bool, ntuple(_ -> 3, N)))
            @test erode(img, se) == erode_ref(img, se)
        end
    end

    img = rand(1:5, 7, 7)
    out = similar(img)
    erode!(out, img)
    @test out == erode(img)

    img = rand(RGB, 7, 7)
    @test_throws ArgumentError("color image is not supported") erode(img)

    # ImageMeta
    @testset "ImageMeta" begin
        A = rand(1:5, 7, 7)
        @test arraydata(erode(ImageMeta(A))) == erode(A)
        @test arraydata(erode(ImageMeta(A); dims=1:2)) == erode(A; dims=1:2)
    end
end


text:
using MPI, Test

MPI.Init()

include("topology.jl")

MPI.Finalize()

@testset "MPI Jobs" begin
  # The code below was modified from the MPI.jl file runtests.jl
  #
  # Code coverage command line options; must correspond to src/julia.h
  # and src/ui/repl.c
  JL_LOG_NONE = 0
  JL_LOG_USER = 1
  JL_LOG_ALL = 2
  coverage_opts = Dict{Int, String}(JL_LOG_NONE => "none",
                                    JL_LOG_USER => "user",
                                    JL_LOG_ALL => "all")
  coverage_opt = coverage_opts[Base.JLOptions().code_coverage]
  testdir = dirname(@__FILE__)

  for (n, f) in [(5, "mpi_connect_1d.jl")
                 (2, "mpi_connect_ell.jl")
                 (3, "mpi_connect.jl")
                 (3, "mpi_connect_stacked.jl")
                 (2, "mpi_connect_stacked_3d.jl")
                 (5, "mpi_connect_sphere.jl")]
    cmd =  `mpiexec -n $n $(Base.julia_cmd()) --startup-file=no --project=$(Base.active_project()) --code-coverage=$coverage_opt $(joinpath(testdir, f))`
    @info "Running MPI test..." n f cmd
    @test (run(cmd); true)
  end
end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
using Test
using Elemental

A = Elemental.DistMatrix(Float32)
B = Elemental.DistMatrix(Float32)
X = Elemental.DistMatrix(Float32)

copy!(A, Float32[2 1; 1 2])
copy!(B, Float32[4, 5])

Elemental.leastSquares!(A, B, X)
@test isapprox(Array(X), [1, 2])

Elemental.ridge!(A, B, 0f0, X)
@test isapprox(Array(X), [1, 2])


text:
VERSION >= v"0.6-" && __precompile__(true)

module Temporal

    if VERSION >= v"0.7-"
        using Dates
    else
        using Base.Dates
    end

    using HTTP
    using JSON
    using RecipesBase

    export
        # Foundational
        TS, ts,
        # Merging/Joining
        ojoin, ijoin, ljoin, rjoin, merge, hcat, vcat, head, tail,
        # Missing Data
        nanrows, nancols, dropnan, fillnan, fillnan!, ffill!, bfill!, linterp!,
        # Operations
        numfun, arrfun, operation,
        ones, zeros, trues, falses, isnan, countnz, sign, round,
        sum, mean, maximum, minimum, prod, cumsum, cumprod, diff, lag, nans,
        shift, pct_change,
        # Aggregation
        mondays, tuesdays, wednesdays, thursdays, fridays, saturdays, sundays,
        bow, eow, bom, eom, boq, eoq, boy, eoy, collapse, apply,
        # IO
        tsread, tswrite,
        quandl, quandl_auth, quandl_meta, quandl_search,
        yahoo_get_crumb, yahoo,
        google,
        # Financial
        has_open, has_high, has_low, has_close, has_volume, is_ohlc, is_ohlcv,
        op, hi, lo, cl, vo, ohlc, ohlcv, hlc, hl, hl2, hlc3, ohlc4,
        # Utilities
        namefix, namefix!, autoidx, findalphanum,
        SANITIZE_NAMES, RANGE_DELIMITER, set_sanitize_names_option, set_range_delimiter_option,
        # Modeling
        acf

    include("ts.jl")
    include("show.jl")
    include("utils.jl")
    include("indexing/getindex.jl")
    include("indexing/setindex.jl")
    include("indexing/stringrange.jl")
    include("combine.jl")
    include("collapse.jl")
    include("operations.jl")
    include("models.jl")
    include("slice.jl")
    include("io.jl")
    include("ohlc.jl")
    include("viz.jl")

end


text:
@given("connection open schema session for database: typedb") do context
    typedb_session = TypeDBClient.CoreSession(client, "typedb" , TypeDBClient.typedb.protocol.Session_Type[:SCHEMA] , TypeDBClient.TypeDBOptions())
    context[:session] = typedb_session
    @expect context[:session] !== nothing
end


text:
module system

export Spin, SpinCollection, CavityMode, CavitySpinCollection

using LinearAlgebra

"""
Abstract base class for all systems defined in this library.

Currently there are following concrete systems:

* Spin
* SpinCollection
* CavityMode
* CavitySpinCollection
"""
abstract type System end

"""
A class representing a single spin.

A spin is defined by its position and its detuning to a main
frequency.

# Arguments
* `position`: A vector defining a point in R3.
* `delta`: Detuning.
"""
struct Spin{T1,T2} <: System
    position::T1
    delta::T2
end
Spin(position::Vector; delta::Real=0) = Spin(position, delta)

"""
A class representing a system consisting of many spins.

# Arguments
* `spins`: Vector of spins.
* `polarizations`: Unit vectors defining the directions of the spins.
* `gammas`: Decay rates.
"""
struct SpinCollection{S<:Spin,P<:Vector,G<:Real} <: System
    spins::Vector{S}
    polarizations::Vector{P}
    gammas::Vector{G}
    function SpinCollection{S,P,G}(spins::Vector{S}, polarizations::Vector{P}, gammas::Vector{G}) where {S,P<:Vector{<:Number},G}
        @assert length(polarizations)==length(spins)
        @assert length(gammas)==length(spins)
        new(spins,normalize.(polarizations),gammas)
    end
end
SpinCollection(spins::Vector{S}, polarizations::Vector{P}, gammas::Vector{G}) where {S,P,G} = SpinCollection{S,P,G}(spins, polarizations, gammas)
SpinCollection(spins::Vector{<:Spin}, polarizations::Vector{<:Vector{<:Number}}, gammas::Number) = SpinCollection(spins, polarizations, [gammas for i=1:length(spins)])
SpinCollection(spins::Vector{<:Spin}, polarizations::Vector{<:Number}, args...) = SpinCollection(spins, [polarizations for i=1:length(spins)], args...)
SpinCollection(spins::Vector{<:Spin}, polarizations; gammas=ones(length(spins))) = SpinCollection(spins, polarizations, gammas)

"""
Create a SpinCollection without explicitly creating single spins.

# Arguments
* `positions`: Vector containing the positions of all single spins.
* `polarizations`: Unit vectors defining the directions of the spins.
* `deltas=0`: Detunings.
* `gammas=1`: Decay rates.
"""
function SpinCollection(positions::Vector{<:Vector{<:Real}}, args...; deltas::Vector=zeros(length(positions)), kwargs...)
    SpinCollection([Spin(positions[i]; delta=deltas[i]) for i=1:length(positions)], args...; kwargs...)
end

"""
A class representing a single mode in a cavity.

# Arguments
* `cutoff`: Number of included Fock states.
* `delta=0` Detuning.
* `eta=0`: Pump strength.
* `kappa=0`: Decay rate.
"""
struct CavityMode{C<:Int,T1<:Number,T2<:Number,T3<:Number} <: System
    cutoff::C
    delta::T1
    eta::T2
    kappa::T3
end
CavityMode(cutoff::Int; delta::Number=0, eta::Number=0, kappa::Number=0) = CavityMode(cutoff,delta,eta,kappa)

"""
A class representing a cavity coupled to many spins.

# Arguments
* `cavity`: A CavityMode.
* `spincollection`: A SpinCollection.
* `g`: A vector defing the coupling strengths between the i-th spin and
    the cavity mode. Alternatively a single number can be given for
    identical coupling for all spins.
"""
struct CavitySpinCollection{C<:CavityMode,S<:SpinCollection,G<:Number} <: System
    cavity::C
    spincollection::S
    g::Vector{G}
    function CavitySpinCollection{C,S,G}(cavity::C, spincollection::S, g::Vector{G}) where {C,S,G}
        @assert length(g) == length(spincollection.spins)
        new(cavity, spincollection, g)
    end
end
CavitySpinCollection(cavity::C, spincollection::S, g::Vector{G}) where {C,S,G} = CavitySpinCollection{C,S,G}(cavity, spincollection, g)
CavitySpinCollection(cavity::CavityMode, spincollection::SpinCollection, g::Real) = CavitySpinCollection(cavity, spincollection, [g for i=1:length(spincollection.spins)])

end # module