text:
# ---
# title: 665. Non-decreasing Array
# id: problem665
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/non-decreasing-array/description/>
# hidden: true
# ---
#
# Given an array `nums` with `n` integers, your task is to check if it could
# become non-decreasing by modifying **at most** `1` element.
#
# We define an array is non-decreasing if `nums[i] <= nums``[i + 1]` holds for
# every `i` (0-based) such that `(0 <= i <= n - 2)`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [4,2,3]
#     Output: true
#     Explanation: You could modify the first 4 to 1 to get a non-decreasing array.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [4,2,1]
#     Output: false
#     Explanation: You can't get a non-decreasing array by modify at most one element.
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 10 ^ 4`
#   * `- 10 ^ 5 <= nums[i] <= 10 ^ 5`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#####################################################################
# BIM
#=
Building Information Modeling is more than just shapes.
Each BIM element is connected to other BIM elements. A wall is connected to its
windows and doors. A floor is connected to its walls, and stairs. Stairs connect
different floors, etc.

This graph of objects needs to be build programmatically and that is one problem.
The other problem is portability, as we want this to work in backends that are
not BIM tools.

A third problem is related to level of detail. Even in the same backend, we might
need to represent BIM elements at different levels of detail.

To solve these problems, we will use a protocol to realize BIM objects. To avoid
confusion, we need to ensure a series of invariants, namely, that entities are
portable above a certain level and non-portable below that level. As an example,
a Slab and a SlabFamily are portable entities.

We will start with the simplest of BIM elements, namely, levels, slabs, walls,
windows and doors.
=#

abstract type BIMElement <: Proxy end
abstract type Measure <: BIMElement end
BIMElements = Vector{<:BIMElement}

# Level
@defproxy(level, Measure, height::Real=0, elements::BIMElements=BIMElement[])
levels_cache = Dict{Real,Level}()
maybe_replace(level::Level) = get!(levels_cache, level.height, level)

convert(::Type{Level}, h::Real) = level(h)

current_levels() = values(level_cache)
default_level = OptionParameter{Level}(level())
default_level_to_level_height = Parameter{Real}(3)
upper_level(lvl=default_level(), height=default_level_to_level_height()) = level(lvl.height + height)
Base.:(==)(l1::Level, l2::Level) = l1.height == l2.height

#default implementation
realize(b::Backend, s::Level) = s.height

export all_levels, default_level, default_level_to_level_height, upper_level

#=
@defproxy(polygonal_mass, BIMShape, points::Locs, height::Real)
@defproxy(rectangular_mass, BIMShape, center::Loc, width::Real, len::Real, height::Real)

@defproxy(column, BIMShape, center::Loc, bottom_level::Any, top_level::Any, family::Any)
=#

abstract type BIMShape <: Shape3D end

#=

We need to provide defaults for a lot of things. For example, we want to specify
a wall that goes through a path without having to specify the kind of wall or its
thickness and height.

This means that, apart from the wall's path, all other wall features will come
from defaults. The base height will be determined by the current level and the
wall height by the current level-to-level height. Finally, the wall thickness,
constituent parts, thermal characteristics, and so on will come from the wall
defaults.  In the case of wall, we will assume that current_wall_defaults() is
a parameter that contains a set of wall parameters.  As an example of use, we
might have:

current_wall_defaults(wall_defaults(thickness=10))

Another option is the definition of different defaults:

thick_wall_defaults = wall_defaults(thickness=10)
thin_wall_defaults = wall_defaults(thickness=5)

which then can be made current:

current_wall_defaults(thin_wall_defaults)

In most cases, the defaults are not just one value, but a bunch of them. For a
beam, we might have:

standard_beam = beam_defaults(width=10, height=20)
current_beam_defaults(standard_beam)

Another useful feature is the ability to adapt defaults. For example:

current_beam_defaults(beam_with(standard_beam, height=20))

Finally, defaults can be created for anything. For example, in a building, we
might want to define a bunch of parameters that are relevant. The syntax is as
follows:

@defaults(building,
    width::Real=20,
    length::Real=30,
    height::Real=50)

In order to access these defaults, we can use the following:

current_building_defaults().width

In some cases, defaults are supported by the backend itself. For example, in
Revit, a wall can be specified using a family. In order to realize the wall
defaults in the current backend, we need to map from the wall parameters to the
corresponding BIM family parameters. This mapping must be described in a
different structure.

For example, a beam element might have a section with a given width and height
but, in Revit, a beam element such as "...\\Structural Framing\\Wood\\M_Timber.rfa"
has, as parameters, the dimensions b and d.  This means that we need a map, such
as Dict(:width => "b", :height => "d")))). So, for a Revit family, we might use:

RevitFamily(
    "C:\\ProgramData\\Autodesk\\RVT 2017\\Libraries\\US Metric\\Structural Framing\\Wood\\M_Timber.rfa",
    Dict(:width => "b", :height => "d"))

To make things more interesting, some families might require instantiation on
different levels. For example, a circular window family in Revit needs to be
loaded using RVTLoadFamily, which requires the name of the family, then it needs
to be instantiated with RVTFamilyElement, which requires the radius of the window,
and finally needs to be inserted on the wall, using RVTInsertWindow, which requires
the opening angle. This might be different for different families, so we need a
flexible way of using the parameters. One hipotesis is to specify those different
moments as different dictionaries.

RevitFamily(
  "C:\\ProgramData\\Autodesk\\RVT 2019\\Libraries\\US Metric\\Windows\\CIRCULAR WINDOW.rfa",
  ("radius_window" => :radius),
  ("angle_window" => :opening_angle))

If no parameters are needed on a particular phase, we might use an empty dictionary
to describe that phase. Given the typical Revit families, defaults parameters seem
to be the best approach here. This means that the previous beam family might be
equivalent to

RevitFamily(
  "C:\\ProgramData\\Autodesk\\RVT 2017\\Libraries\\US Metric\\Structural Framing\\Wood\\M_Timber.rfa",
  ("b" => :width, "d" => :height),
  ())

However, the same beam might have a different mapping in a different backend.
This means that we need another mapping to support different backends. One
possibility is to use something similar to:

backend_family(
    revit => RevitFamily(
        "C:\\ProgramData\\Autodesk\\RVT 2017\\Libraries\\US Metric\\Structural Framing\\Wood\\M_Timber.rfa",
        ("b" => :width, "d" => :height)),
    archicad => ArchiCADFamily(
        "BeamElement",
        ("size_x" => :width, "size_y" => :height)),
    autocad => AutoCADFamily())

Then, we need an operation that instantiates a family. This can be done on two different
levels: (1) from a backend-specific family (e.g., RevitFamily), for example:

beam_family = RevitFamily(
    "C:\\ProgramData\\Autodesk\\RVT 2017\\Libraries\\US Metric\\Structural Framing\\Wood\\M_Timber.rfa",
    ("b" => :width, "d" => :height))

current_beam_defaults(beam_family_instance(beam_family, width=10, height=20)

or from a generic backend family, for example:

beam_family = backend_family(
    revit => RevitFamily(
        "C:\\ProgramData\\Autodesk\\RVT 2017\\Libraries\\US Metric\\Structural Framing\\Wood\\M_Timber.rfa",
        ("b" => :width, "d" => :height)),
    archicad => ArchiCADFamily(
        "BeamElement"
        ("size_x" => :width, "size_y" => :height)),
    autocad => AutoCADFamily())

current_beam_defaults(beam_family_instance(beam_family, width=10, height=20)

In this last case, the generic family will use the top_backend value to identify
which family to use.

Another important feature is the use of a delegation-based implementation for
family instances. This means that we might do

current_beam_defaults(beam_family_instance(current_beam_defaults(), width=20)

to instantiate a family that uses, by default, the same parameter values used by
another family instance.

Finally, because we might need to know which families are available (e.g., for
metaprogramming purposes), we need a database of families for each BIM element.

A relatively simple approach is to have a registry for families:

family_registry(WallFamily) -> Vector{WallFamily}
register_family(family)
=#

export Family, FamilyInstance, family, family_ref

abstract type Family <: Proxy end
abstract type FamilyInstance <: Family end

family(f::Family) = f
family(f::FamilyInstance) = f.family

#=
Caching of families is more complex than caching of shapes because families
are created independently of the current backend and then they are lazily
instantiated, possibly multiples times as the current backend is being changed.
=#

family_ref(s::BIMShape) =
   family_ref(backend(s), s.family)
family_ref(b::Backend, f::Family) =
  get!(f.ref, b) do
    realize(b, f)
  end

# CAD tools that support Layers might benefit from this typical implementation:
struct LayerFamily <: Family
  name::String
  color::RGB
  ref::Parameter{Any}
end

export LayerFamily, layer_family
layer_family(name, color::RGB=rgb(1,1,1)) =
  LayerFamily(name, color, Parameter{Any}(nothing))

backend_get_family_ref(b::Backend, f::Family, af::LayerFamily) =
  b_layer(b, af.name, true, af.color)

macro deffamily(name, parent, fields...)
  name_str = string(name)
  abstract_name = esc(Symbol(string))
  struct_name = esc(Symbol(string(map(uppercasefirst,split(name_str,'_'))...)))
  # We always add a name field
  fields = [:(name::String=$(name_str)), fields...]
  field_names = map(field -> field.args[1].args[1], fields)
  field_types = map(field -> field.args[1].args[2], fields)
  field_inits = map(field -> field.args[2], fields)
  field_renames = map(esc ∘ Symbol ∘ uppercasefirst ∘ string, field_names)
  field_replacements = Dict(zip(field_names, field_renames))
  struct_fields = map((name,typ) -> :($(name) :: $(typ)), field_names, field_types)
#  opt_params = map((name,typ,init) -> :($(name) :: $(typ) = $(init)), field_renames, field_types, field_inits)
#  key_params = map((name,typ,rename) -> :($(name) :: $(typ) = $(rename)), field_names, field_types, field_renames)
#  mk_param(name,typ) = Expr(:kw, Expr(:(::), name, typ))
  mk_param(name,typ,init) = Expr(:kw, Expr(:(::), name, typ), init)
  opt_params = map(mk_param, field_renames, field_types, map(init -> replace_in(init, field_replacements), field_inits))
  key_params = map(mk_param, field_names, field_types, field_renames)
  instance_params = map(mk_param, field_names, field_types, map(name -> :(family.$(name)), field_names))
  constructor_name = esc(name)
  instance_name = esc(Symbol(name_str, "_element")) #"_instance")) beam_family_element or beam_family_instance?
  default_name = esc(Symbol("default_", name_str))
  predicate_name = esc(Symbol("is_", name_str))
  selector_names = map(field_name -> esc(Symbol(name_str, "_", string(field_name))), field_names)
  with_name = esc(Symbol("with_", name_str))
  quote
    export $(constructor_name), $(instance_name), $(default_name), $(predicate_name), $(struct_name), $(with_name)
    struct $struct_name <: $parent
      $(struct_fields...)
      based_on::Union{Family, Nothing}
      implemented_as::IdDict{<:Backend, <:Family}
      ref::IdDict{<:Backend, Any}
      data::BackendParameter
    end
    $(constructor_name)($(opt_params...);
                        $(key_params...),
                        based_on=nothing,
                        implemented_as=IdDict{Backend, Family}(),
                        data=BackendParameter()) =
      $(struct_name)($(field_names...), based_on, implemented_as, IdDict{Backend, Any}(), data)
    $(instance_name)(family:: Family, implemented_as=copy(family.implemented_as); $(instance_params...)) =
      $(struct_name)($(field_names...), family, implemented_as, IdDict{Backend, Any}(), copy(family.data))
    $(default_name) = Parameter{$struct_name}($(constructor_name)())
    $(predicate_name)(v::$(struct_name)) = true
    $(predicate_name)(v::Any) = false
    $(with_name)(f::Function; family :: $(struct_name) = $(default_name)(), $(instance_params...)) =
      with(f, $(default_name), $(instance_name)(family; $([:($(esc(p))=$(esc(p))) for p in field_names]...)))
#    $(map((selector_name, field_name) -> :($(selector_name)(v::$(struct_name)) = v.$(field_name)),
#          selector_names, field_names)...)
    KhepriBase.meta_program(v::$(struct_name)) =
        Expr(:call, $(Expr(:quote, name)), $(map(field_name -> :(meta_program(v.$(field_name))), field_names)...))
    KhepriBase.meta_program(v::Parameter{$struct_name}) =
        Expr(:call, $(Expr(:quote, default_name)))
  end
end

export set_family
const set_family = set_on!

#family_ref(b::Backend, m::Family) = ref(b, m).value
#family_ref(b::Backend, s::BIMShape) = family_ref(b, s.family)

#=

We also need to register families so that we know, at runtime, which ones are available.

We need BIMElement
avaliable_families(element::BIMElement) => Array of families
default_family_parameter(element::BIMElement) => Parameter

ALSO

Give families names so that the interface can present them

We need to detect using traceability the Khepri primitive that is called. So we cannot exclude the
Khepri module.

=#

# When dispatching a BIM operation to a backend, we also need to dispatch the family

backend_family(b::Backend, family::Family) =
  get(family.implemented_as, b) do
    isnothing(family.based_on) ? # this is not a family_element (nor a derivation of a family_element)
      error("Family $(family) is missing the implementation for backend $(b)") :
      backend_family(b, family.based_on)
  end

copy_struct(s::T) where T = T([getfield(s, k) for k ∈ fieldnames(T)]...)

# Backends will install their own families on top of the default families, e.g.,
# set_backend_family(default_beam_family(), revit, revit_beam_family)
set_backend_family(family::Family, backend::Backend, backend_family::Family) =
  begin
    family.implemented_as[backend]=backend_family
    delete!(family.ref, backend) #force recreation
  end

realize(b::Backend, f::Family) =
  backend_get_family_ref(b, f, backend_family(b, f))

backend_get_family_ref(b::Backend, f::Family, bf) = bf
export backend_family, set_backend_family

#=
We can now define specific families for slabs, beams, etc., the corresponding
specific building elements and even, when possible, default implementations.
=#

@deffamily(slab_family, Family,
  thickness::Real=0.2,
  coating_thickness::Real=0.0,
  bottom_material::Material=material_concrete,
  top_material::Material=material_concrete,
  side_material::Material=material_concrete)

slab_family_elevation(b::Backend, family::SlabFamily) =
  family.coating_thickness - family.thickness
slab_family_thickness(b::Backend, family::SlabFamily) =
  family.coating_thickness + family.thickness

@defproxy(slab, BIMShape, region::Region=rectangular_path(),
          level::Level=default_level(), family::SlabFamily=default_slab_family())

# Default implementation: dispatch on the slab elements
realize(b::Backend, s::Slab) =
  b_slab(b, s.region, s.level, s.family)

#=
export add_slab_opening
add_slab_opening(s::Slab=required(), contour::ClosedPath=circular_path()) =
    let b = backend(s)
        push!(s.openings, contour)
        if realized(s)
            set_ref!(s, realize_slab_openings(b, s, ref(s), [contour]))
        end
        s
    end

realize_slab_openings(b::Backend, s::Slab, s_ref, openings) =
    let s_base_height = s.level.height,
        s_thickness = slab_family_thickness(b, s.family)
        for opening in openings
            op_path = translate(opening, vz(s_base_height-1.1*s_thickness))
            op_ref = ensure_ref(b, backend_slab(b, op_path, s_thickness*1.2))
            s_ref = ensure_ref(b, subtract_ref(b, s_ref, op_ref))
        end
        s_ref
    end
=#

# Roof

@deffamily(roof_family, Family,
    thickness::Real=0.2,
    coating_thickness::Real=0.0,
    bottom_material::Material=material_concrete,
    top_material::Material=material_concrete,
    side_material::Material=material_concrete)

slab_family_elevation(b::Backend, family::RoofFamily) = 0
slab_family_thickness(b::Backend, family::RoofFamily) =
  family.coating_thickness + family.thickness

@defproxy(roof, BIMShape, region::Region=rectangular_path(),
          level::Level=default_level(), family::RoofFamily=default_roof_family())
realize(b::Backend, s::Roof) =
  b_roof(b, s.region, s.level, s.family)

# Panel

@deffamily(panel_family, Family,
  thickness::Real=0.02,
  material::Material=material_glass)

@defproxy(panel, BIMShape, region::Region=rectangular_path(), family::PanelFamily=default_panel_family())

realize(b::Backend, s::Panel) =
  b_panel(b, s.region, s.family)

#=

A wall contains doors and windows

=#

# Wall

@deffamily(wall_family, Family,
    thickness::Real=0.2,
    left_coating_thickness::Real=0.0,
    right_coating_thickness::Real=0.0,
    right_material::Material=material_plaster,
    left_material::Material=material_plaster,
    side_material::Material=material_plaster)

@defproxy(wall, BIMShape, path::Path=rectangular_path(),
          bottom_level::Level=default_level(),
          top_level::Level=upper_level(convert(Level, bottom_level)),
          family::WallFamily=default_wall_family(),
          offset::Real=is_closed_path(path) ? 1/2 : 0, # offset is relative to the thickness
          doors::Shapes=Shape[], windows::Shapes=Shape[])

# To deal with incremental creation, it is necessary to use transactions.
export with_wall
with_wall(f, args...) =
  with_transaction() do
    let w = wall(args...)
      f(w)
      w
    end
  end

# The protocol starts by identifying the approach to use. It can be either
# based on Boolean operations or on the construction of polygonal elements.
# Then, for each approach, the appropriate implementation is selected.
export HasBooleanOps, has_boolean_ops

struct HasBooleanOps{T} end
# By default, we rely on boolean operations
# By default, we DON'T rely on boolean operations
has_boolean_ops(::Type{<:Backend}) = HasBooleanOps{false}()

realize(b::B, w::Wall) where B<:Backend =
  realize(has_boolean_ops(B), b, w)

realize(::HasBooleanOps{true}, b::Backend, w::Wall) =
  with_material_as_layer(b, w.family) do
    realize_wall_openings(b, w, realize_wall_no_openings(b, w), [w.doors..., w.windows...])
  end

realize_wall_no_openings(b::Backend, w::Wall) =
  let w_base_height = w.bottom_level.height,
      w_height = w.top_level.height - w_base_height,
      w_path = translate(w.path, vz(w_base_height)),
      r_thickness = r_thickness(w),
      l_thickness = l_thickness(w)
    ensure_ref(b, backend_wall(b, w_path, w_height, l_thickness, r_thickness, w.family))
  end

realize_wall_openings(b::Backend, w::Wall, w_ref, openings) =
  let w_base_height = w.bottom_level.height,
      w_height = w.top_level.height - w_base_height,
      w_path = translate(w.path, vz(w_base_height)),
      r_thickness = r_thickness(w),
      l_thickness = l_thickness(w)
    for opening in openings
      w_ref = realize_wall_opening(b, w_ref, w_path, l_thickness, r_thickness, opening, w.family)
      ref(b, opening)
    end
    w_ref
  end

realize_wall_opening(b::Backend, w_ref, w_path, l_thickness, r_thickness, op, family) =
  let op_base_height = op.loc.y,
      op_height = op.family.height,
      op_path = translate(subpath(w_path, op.loc.x, op.loc.x + op.family.width), vz(op_base_height)),
      op_ref = ensure_ref(b, backend_wall(b, op_path, op_height, l_thickness+0.1, r_thickness+0.1, family))
    ensure_ref(b, subtract_ref(b, w_ref, op_ref))
  end

# For backends that do not support boolean operations, we use a different approach

realize(::HasBooleanOps{false}, b::Backend, w::Wall) =
  let w_base_height = w.bottom_level.height,
      w_height = w.top_level.height - w_base_height,
      r_thickness = r_thickness(w),
      l_thickness = l_thickness(w),
      w_path = translate(w.path, vz(w_base_height)),
      w_paths = subpaths(w_path),
      r_w_paths = subpaths(offset(w_path, -r_thickness)),
      l_w_paths = subpaths(offset(w_path, l_thickness)),
      openings = [w.doors..., w.windows...],
      prevlength = 0,
      matright = w.family.right_material,
      matleft = w.family.left_material,
      matside = w.family.side_material,
      refs = []
    for (w_seg_path, r_w_path, l_w_path) in zip(w_paths, r_w_paths, l_w_paths)
      let currlength = prevlength + path_length(w_seg_path),
          c_r_w_path = closed_path_for_height(r_w_path, w_height),
          c_l_w_path = closed_path_for_height(l_w_path, w_height)
        append!(refs, materialize_path(b, c_r_w_path, c_l_w_path, matside))
        openings = filter(openings) do op
          if prevlength <= op.loc.x < currlength ||
             prevlength <= op.loc.x + op.family.width <= currlength # contained (at least, partially)
            let op_height = op.family.height,
                op_at_start = op.loc.x <= prevlength,
                op_at_end = op.loc.x + op.family.width >= currlength,
                op_path = subpath(w_path,
                                  max(prevlength, op.loc.x),
                                  min(currlength, op.loc.x + op.family.width)),
                r_op_path = offset(op_path, -r_thickness),
                l_op_path = offset(op_path,  l_thickness),
                fixed_r_op_path =
                  open_polygonal_path([path_start(op_at_start ? r_w_path : r_op_path),
                                       path_end(op_at_end ? r_w_path : r_op_path)]),
                fixed_l_op_path =
                  open_polygonal_path([path_start(op_at_start ? l_w_path : l_op_path),
                                       path_end(op_at_end ? l_w_path : l_op_path)]),
                c_r_op_path = closed_path_for_height(translate(fixed_r_op_path, vz(op.loc.y)), op_height),
                c_l_op_path = closed_path_for_height(translate(fixed_l_op_path, vz(op.loc.y)), op_height)
              append!(refs, materialize_path(b, reverse(c_r_op_path), reverse(c_l_op_path), matside))
              c_r_w_path, c_l_w_path = subtract_paths(b, c_r_w_path, c_l_w_path, c_r_op_path, c_l_op_path)
              # preserve if not totally contained
              ! (op.loc.x >= prevlength && op.loc.x + op.family.width <= currlength)
            end
          else
            true
          end
        end
        prevlength = currlength
        push!(refs, materialize_path(b, reverse(c_l_w_path), matleft))
        push!(refs, materialize_path(b, c_r_w_path, matright))
      end
    end
    refs
  end

closed_offsetted_path(path, v) =
  let ps = path_vertices(path)
    closed_polygonal_path([ps..., reverse(map(p -> p+v, ps))...])
  end

closed_path_for_height(path, h) =
  closed_offsetted_path(path, vz(h))

subtract_paths(b::Backend, c_r_w_path, c_l_w_path, c_r_op_path, c_l_op_path) =
  region(c_r_w_path, c_r_op_path), region(c_l_w_path, c_l_op_path)
  #=
  let idxs = closest_vertices_indexes(path_vertices(c_r_w_path), path_vertices(c_r_op_path))
    closed_polygonal_path(
      inject_polygon_vertices_at_indexes(path_vertices(c_r_w_path), path_vertices(c_r_op_path), idxs)),
    closed_polygonal_path(
      inject_polygon_vertices_at_indexes(path_vertices(c_l_w_path), path_vertices(c_l_op_path), idxs))
  end=#

#=
Walls can be joined. That is very important because the wall needs to have
uniform thickness along the entire path.
=#
export join_walls
join_walls(wall1, wall2) =
  if wall1.bottom_level != wall2.bottom_level
    error("Walls with different bottom levels")
  elseif wall1.top_level != wall2.top_level
    error("Walls with different top levels")
  elseif wall1.family != wall2.family
    error("Walls with different families")
  elseif wall1.offset != wall2.offset
    error("Walls with different offsets")
  else
    let w = wall(join_paths(wall1.path, wall2.path),
                 wall1.bottom_level, wall1.top_level,
                 wall1.family, wall1.offset),
        len = path_length(wall1.path)
      for (es,l) in ((wall1.doors, 0), (wall2.doors, len))
        for e in es
          add_door(w, e.loc+vx(l), e.family)
        end
      end
      for (es,l) in ((wall1.windows, 0), (wall2.windows, len))
        for e in es
          add_window(w, e.loc+vx(l), e.family)
        end
      end
      for w in (wall1, wall2)
        delete_shapes(w.doors)
        delete_shapes(w.windows)
        delete_shape(w)
      end
      w
    end
  end

join_walls(walls...) =
  reduce(join_walls, walls)

# Right and Left considering observer looking along with curve direction
# a non-closed wall should have a wall offset of zero and a r_thickness and a l_thickness of 1/2*thickness
# a closed wall should have a wall offset of 1/2 and a r_thickness of zero and a l_thickness of 1*thickness
# Thus, we have:
r_thickness(offset, thickness) = (1/2 - offset)*thickness
l_thickness(offset, thickness) = (1/2 + offset)*thickness

r_thickness(w::Wall) = r_thickness(w.offset, w.family.thickness + w.family.right_coating_thickness)
l_thickness(w::Wall) = l_thickness(w.offset, w.family.thickness + w.family.left_coating_thickness)

# Door

@deffamily(door_family, Family,
  width::Real=1.0,
  height::Real=2.0,
  thickness::Real=0.05,
  right_material::Material=material_wood,
  left_material::Material=material_wood,
  side_material::Material=material_wood)

@defproxy(door, BIMShape, wall::Wall=required(), loc::Loc=u0(), flip_x::Bool=false, flip_y::Bool=false, family::DoorFamily=default_door_family())

# Window

@deffamily(window_family, Family,
  width::Real=1.0,
  height::Real=1.0,
  thickness::Real=0.05,
  right_material::Material=material_glass,
  left_material::Material=material_glass,
  side_material::Material=material_glass)

@defproxy(window, BIMShape, wall::Wall=required(), loc::Loc=u0(), flip_x::Bool=false, flip_y::Bool=false, family::WindowFamily=default_window_family())

realize(b::Backend, s::Union{Door, Window}) =
  let base_height = s.wall.bottom_level.height + s.loc.y,
      height = s.family.height,
      subpath = translate(subpath(s.wall.path, s.loc.x, s.loc.x + s.family.width), vz(base_height)),
      r_thickness = r_thickness(s.wall),
      l_thickness = l_thickness(s.wall),
      thickness = s.family.thickness
    b_wall(b, subpath, height, (l_thickness - r_thickness + thickness)/2, (r_thickness - l_thickness + thickness)/2,
           s.family)
  end
##

export add_door
add_door(w::Wall=required(), loc::Loc=u0(), family::DoorFamily=default_door_family()) =
  let d = door(w, loc, family=family)
    push!(w.doors, d)
    delete_shape(w)
    force_realize(w)
    w
  end

export add_window
add_window(w::Wall=required(), loc::Loc=u0(), family::WindowFamily=default_window_family()) =
  let d = window(w, loc, family=family)
    push!(w.windows, d)
    delete_shape(w)
    force_realize(w)
    w
  end

#=
A curtain wall is a special kind of wall that is made of a frame with windows.
=#

@deffamily(curtain_wall_frame_family, Family,
  width::Real=0.1,
  depth::Real=0.1,
  depth_offset::Real=0.25,
  material::Material=material_metal)

@deffamily(curtain_wall_family, Family,
  max_panel_dx::Real=1,
  max_panel_dy::Real=2,
  panel::PanelFamily=panel_family(thickness=0.05),
  boundary_frame::CurtainWallFrameFamily=
    curtain_wall_frame_family(width=0.1,depth=0.1,depth_offset=0.25),
  mullion_frame::CurtainWallFrameFamily=
    curtain_wall_frame_family(width=0.08,depth=0.09,depth_offset=0.2),
  transom_frame::CurtainWallFrameFamily=
    curtain_wall_frame_family(width=0.06,depth=0.1,depth_offset=0.11))

@defproxy(curtain_wall, BIMShape,
          path::Path=rectangular_path(),
          bottom_level::Level=default_level(),
          top_level::Level=upper_level(bottom_level),
          family::CurtainWallFamily=default_curtain_wall_family(),
          offset::Real=0.0)
curtain_wall(p0::Loc, p1::Loc;
     bottom_level::Level=default_level(),
     top_level::Level=upper_level(bottom_level),
     family::CurtainWallFamily=default_curtain_wall_family(),
     offset::Real=0.0) =
  curtain_wall([p0, p1], bottom_level=bottom_level, top_level=top_level,
         family=family, offset=offset)

realize(b::Backend, s::CurtainWall) =
  with_material_as_layer(b, s.family) do
    let th = s.family.panel.thickness,
        bfw = s.family.boundary_frame.width,
        bfd = s.family.boundary_frame.depth,
        bfdo = s.family.boundary_frame.depth_offset,
        mfw = s.family.mullion_frame.width,
        mfd = s.family.mullion_frame.depth,
        mdfo = s.family.mullion_frame.depth_offset,
        tfw = s.family.transom_frame.width,
        tfd = s.family.transom_frame.depth,
        tfdo = s.family.transom_frame.depth_offset,
        path = curtain_wall_path(b, s, s.family.panel),
        path_length = path_length(path),
        bottom = level_height(b, s.bottom_level),
        top = level_height(b, s.top_level),
        height = top - bottom,
        x_panels = ceil(Int, path_length/s.family.max_panel_dx),
        y_panels = ceil(Int, height/s.family.max_panel_dy),
        refs = []
      push!(refs, backend_curtain_wall(b, s, subpath(path, bfw, path_length-bfw), bottom+bfw, height-2*bfw, th/2, th/2, :panel))
      push!(refs, backend_curtain_wall(b, s, path, bottom, bfw, l_thickness(bfdo, bfd), r_thickness(bfdo, bfd), :boundary_frame))
      push!(refs, backend_curtain_wall(b, s, path, top-bfw, bfw, l_thickness(bfdo, bfd), r_thickness(bfdo, bfd), :boundary_frame))
      push!(refs, backend_curtain_wall(b, s, subpath(path, 0, bfw), bottom+bfw, height-2*bfw, l_thickness(bfdo, bfd), r_thickness(bfdo, bfd), :boundary_frame))
      push!(refs, backend_curtain_wall(b, s, subpath(path, path_length-bfw, path_length), bottom+bfw, height-2*bfw, l_thickness(bfdo, bfd), r_thickness(bfdo, bfd), :boundary_frame))
      for i in 1:y_panels-1
        l = height/y_panels*i
        sub = subpath(path, bfw, path_length-bfw)
        push!(refs, backend_curtain_wall(b, s, sub, bottom+l-tfw/2, tfw, l_thickness(tfdo, tfd), r_thickness(tfdo, tfd), :transom_frame))
      end
      for i in 1:x_panels-1
        l = path_length/x_panels*i
        push!(refs, backend_curtain_wall(b, s, subpath(path, l-mfw/2, l+mfw/2), bottom+bfw, height-2*bfw, l_thickness(mdfo, mfd), r_thickness(mdfo, mfd), :mullion_frame))
      end
      [ensure_ref(b,r) for r in refs]
    end
  end

# By default, curtain wall panels are planar
curtain_wall_path(b::Backend, s::CurtainWall, panel_family::Family) =
  s.path

curtain_wall_path(b::Backend, s::CurtainWall, panel_family::PanelFamily) =
  let path_length = path_length(s.path),
      x_panels = ceil(Int, path_length/s.family.max_panel_dx),
      pts = map(t->in_world(location_at_length(s.path, t)),
                division(0, path_length, x_panels))
    polygonal_path(pts)
  end

#
# We need to redefine the default method (maybe add an option to the macro to avoid defining the meta_program)
# This needs to be fixed for windows
#=
meta_program(w::Wall) =
    if isempty(w.doors)
        Expr(:call, :wall,
             meta_program(w.path),
             meta_program(w.bottom_level),
             meta_program(w.top_level),
             meta_program(w.family))
    else
        let door = w.doors[1]
            Expr(:call, :add_door,
                 meta_program(wall(w.path, w.bottom_level, w.top_level, w.family, w.doors[2:end], w.windows)),
                 meta_program(door.loc),
                 meta_program(door.family))
        end
    end
=#

# Beam
# Beams are mainly horizontal elements. By default, a beam is aligned along its top axis
@deffamily(beam_family, Family,
  profile::ClosedPath=top_aligned_rectangular_profile(1, 2),
  material::Material=material_metal)

@defproxy(beam, BIMShape, cb::Loc=u0(), h::Real=1, angle::Real=0, family::BeamFamily=default_beam_family())
beam(cb::Loc, ct::Loc, Angle::Real=0, Family::BeamFamily=default_beam_family(); angle::Real=Angle, family::BeamFamily=Family) =
    let (c, h) = position_and_height(cb, ct)
      beam(c, h, angle, family)
    end

realize(b::Backend, s::Beam) =
  b_beam(b, s.cb, s.h, s.angle, s.family)

# Column
# Columns are mainly vertical elements. A column has its center axis aligned with a line defined by two points

@deffamily(column_family, Family,
  profile::ClosedPath=rectangular_profile(0.2, 0.2),
  material::Material=material_concrete)

@defproxy(free_column, BIMShape, cb::Loc=u0(), h::Real=1, angle::Real=0, family::ColumnFamily=default_column_family())
free_column(cb::Loc, ct::Loc, Angle::Real=0, Family::ColumnFamily=default_column_family(); angle::Real=Angle, family::ColumnFamily=Family) =
  let (c, h) = position_and_height(cb, ct)
    free_column(c, h, angle, family)
  end

realize(b::Backend, s::FreeColumn) =
  b_free_column(b, s.cb, s.h, s.angle, s.family)

@defproxy(column, BIMShape, cb::Loc=u0(), angle::Real=0,
  bottom_level::Level=default_level(), top_level::Level=upper_level(bottom_level),
  family::ColumnFamily=default_column_family())

realize(b::Backend, s::Column) =
  b_column(b, s.cb, s.angle, s.bottom_level, s.top_level, s.family)

# Tables and chairs

@deffamily(table_family, Family,
    length::Real=1.6,
    width::Real=0.9,
    height::Real=0.75,
    top_thickness::Real=0.05,
    leg_thickness::Real=0.05,
    material::Material=material_wood)

@deffamily(chair_family, Family,
    length::Real=0.4,
    width::Real=0.4,
    height::Real=1.0,
    seat_height::Real=0.5,
    thickness::Real=0.05,
    material::Material=material_wood)

@deffamily(table_chair_family, Family,
    table_family::TableFamily=default_table_family(),
    chair_family::ChairFamily=default_chair_family(),
    chairs_top::Int=1,
    chairs_bottom::Int=1,
    chairs_right::Int=2,
    chairs_left::Int=2,
    spacing::Real=0.7)

@defproxy(table, BIMShape, loc::Loc=u0(), level::Level=default_level(), family::TableFamily=default_table_family())
table(loc::Loc, angle::Real, level=default_level(), family::TableFamily=default_table_family()) =
  table(loc_from_o_phi(loc, angle), level, family)
realize(b::Backend, s::Table) =
  let tf = s.family
    b_table(b, add_z(s.loc, s.level.height),
            tf.length, tf.width, tf.height,
            tf.top_thickness, tf.leg_thickness,
            material_ref(b, tf.material))
  end

@defproxy(chair, BIMShape, loc::Loc=u0(), level::Level=default_level(), family::ChairFamily=default_chair_family())
chair(loc::Loc, angle::Real, level::Level=default_level(), family::ChairFamily=default_chair_family()) =
  chair(loc_from_o_phi(loc, angle), level, family)
realize(b::Backend, s::Chair) =
  let cf = s.family
    b_chair(b, add_z(s.loc, s.level.height),
            cf.length, cf.width, cf.height,
            cf.seat_height, cf.thickness,
            material_ref(b, cf.material))
  end

@defproxy(table_and_chairs, BIMShape, loc::Loc=u0(), level::Level=default_level(), family::TableChairFamily=default_table_chair_family())
table_and_chairs(loc::Loc, angle::Real, level::Level=default_level(), family::TableChairFamily=default_table_chair_family()) =
  table_and_chairs(loc_from_o_phi(loc, angle), level, family)
realize(b::Backend, s::TableAndChairs) =
  let f = s.family,
      tf = f.table_family,
      cf = f.chair_family,
      tmat = material_ref(b, tf.material),
      cmat = material_ref(b, cf.material)
    b_table_and_chairs(b,
      add_z(s.loc, s.level.height),
      p->b_table(b, p, tf.length, tf.width, tf.height, tf.top_thickness, tf.leg_thickness, tmat),
      p->b_chair(b, p, cf.length, cf.width, cf.height, cf.seat_height, cf.thickness, cmat),
      tf.length,
      tf.width,
      f.chairs_top,
      f.chairs_bottom,
      f.chairs_right,
      f.chairs_left,
      f.spacing)
  end

# Lights

@defproxy(pointlight, BIMShape, loc::Loc=z(3), color::RGB=rgb(1,1,1), range::Real=10, intensity::Real=4, level::Level=default_level())

realize(b::Backend, s::Pointlight) =
  b_pointlight(b, add_z(s.loc, s.level.height), s.color, s.range, s.intensity)

@defproxy(spotlight, BIMShape, loc::Loc=z(3), dir::Vec=vz(-1), hotspot::Real=pi/4, falloff::Real=pi/3)

realize(b::Backend, s::Spotlight) =
  b_spotlight(b, s.loc, s.dir, s.hotspot, s.falloff)

@defproxy(ieslight, BIMShape, file::String=required(), loc::Loc=z(3), dir::Vec=vz(-1), alpha::Real=0, beta::Real=0, gamma::Real=0)

realize(b::Backend, s::Ieslight) =
  b_ieslight(b, s.file, s.loc, s.dir, s.alpha, s.beta, s.gamma)

#################################
# Node support
Base.@kwdef struct TrussNodeSupport
    ux::Bool=false
    uy::Bool=false
    uz::Bool=false
    rx::Bool=false
    ry::Bool=false
    rz::Bool=false
end

const truss_node_support = TrussNodeSupport

@deffamily(truss_node_family, Family,
    radius::Real=0.03,
    support::Any=false,
    material::Material=material_metal)

@deffamily(truss_bar_family, Family,
    radius::Real=0.02,
    inner_radius::Real=0,
    material::Material=material_metal)

truss_bar_family_cross_section_area(f::TrussBarFamily) =
  error("This should be computed by the backend family") #truss_bar_family_cross_section_area(back)

# We need a few families by default
export free_truss_node_family, fixed_truss_node_family, truss_node_support

free_truss_node_family =
  truss_node_family_element(default_truss_node_family(),
                            support=truss_node_support())
fixed_truss_node_family =
  truss_node_family_element(default_truss_node_family(),
                            support=truss_node_support(ux=true, uy=true, uz=true))

@defproxy(truss_node, BIMShape, p::Loc=u0(), family::TrussNodeFamily=default_truss_node_family())
@defproxy(truss_bar, BIMShape, p0::Loc=u0(), p1::Loc=u0(), angle::Real=0, family::TrussBarFamily=default_truss_bar_family())

realize(b::Backend, s::TrussNode) =
  truss_node_is_supported(s) ?
    [b_truss_node(b, s.p, s.family), b_truss_node_support(b, s.p, s.family)] :
    b_truss_node(b, s.p, s.family)

realize(b::Backend, s::TrussBar) =
  b_truss_bar(b, s.p0, s.p1, s.family)

export truss_nodes, truss_bars
truss_nodes(ps, family=default_truss_node_family()) =
  [truss_node(p, family) for p in ps]
truss_bars(ps, qs, family=default_truss_bar_family()) =
  [truss_bar(p, q, 0, family) for (p, q) in zip(ps, qs)]

export truss_node_is_supported
truss_node_is_supported(n) =
  let s = n.family.support
    s != false && (s.ux || s.uy || s.uz || s.rx || s.ry || s.rz)
  end

truss_bar_cross_section_area(s::TrussBar) =
  truss_bar_family_cross_section_area(family_ref(top_backend(), s.family))
truss_bar_volume(s::TrussBar) =
  truss_bar_cross_section_area(s)*distance(s.p0, s.p1)

# Should we merge coincident nodes?
export merge_coincident_truss_nodes,
       coincident_truss_nodes_distance,
       merge_coincident_truss_bars,
       maybe_merged_node,
       maybe_merged_bar

const merge_coincident_truss_nodes = Parameter(true)
const coincident_truss_nodes_distance = Parameter(1e-6)

maybe_merged_node(b::Backend, s::TrussNode) =
  # We are allowed to replace a node that we just created with one that already exists
  let epsilon = coincident_truss_nodes_distance(),
      p = s.p
    for n in b.truss_nodes
      if distance(n.p, p) < epsilon
        merge_coincident_truss_nodes() || error("Coincident nodes $(s) and $(n) at $(p)")
        return n
      end
    end
    b.realized(false)
    push!(b.truss_nodes, s)
    s
  end

# Should we merge coincident bars?
const merge_coincident_truss_bars = Parameter(true)

maybe_merged_bar(b::Backend, s::TrussBar) =
  # We are allowed to replace a bar that we just created with one that already exists
  let epsilon = coincident_truss_nodes_distance(),
      p0 = s.p0,
      p1 = s.p1
    for b in b.truss_bars
      if (distance(b.p0, p0) < epsilon && distance(b.p1, p1) < epsilon) ||
         (distance(b.p0, p1) < epsilon && distance(b.p1, p0) < epsilon)
        merge_coincident_truss_bars() || error("Coincident bars $(s) and $(n) at $(p0) and $(p1)")
        return b
      end
    end
    b.realized(false)
    push!(b.truss_bars, s)
    s
  end

# Many analysis tools prefer a simplified node-and-bar representation.
# To merge nodes and bars we need a different structure. Maybe we should merge
# this with the BIM information
export TrussNodeData, truss_node_data, TrussBarData, truss_bar_data

struct TrussNodeData
    id::Int
    loc::Loc
    family::Any
    load::Any
end

truss_node_data(id::Int, loc::Loc, family::Any, load::Vec) =
  TrussNodeData(id, loc, family, load)

#
struct TrussBarData
    id::Int
    node1::TrussNodeData
    node2::TrussNodeData
    rotation::Real
    family::Any
end

truss_bar_data(id::Int, node0::TrussNodeData, node1::TrussNodeData, rotation::Real, family::Any) =
  TrussBarData(id, node0, node1, rotation, family)

export process_nodes, process_bars
process_nodes(nodes, load=vz(0), loads_points=Dict()) =
  let point_loads = Dict()
    for k in keys(loads_points)
      for p in loads_points[k]
        point_loads[p] = k
      end
    end
    [truss_node_data(i, in_world(node.p), node.family, load + get(point_loads, node.p, vz(0)))
     for (i, node) in enumerate(nodes)]
  end
process_bars(bars, processed_nodes) =
  let epsilon = coincident_truss_nodes_distance(),
      node_data_near(loc) =
        begin
          for nd in processed_nodes
            distance(loc, nd.loc) < epsilon && return nd
          end
          error("Bar without corresponding node at location $(loc)!")
        end
    [truss_bar_data(
      i,
      node_data_near(in_world(bar.p0)),
      node_data_near(in_world(bar.p1)),
      bar.angle,
      bar.family)
     for (i, bar) in enumerate(bars)]
  end

# Analysis
@defcb truss_analysis(load::Vec=vz(-1e5), self_weight::Bool=false, point_loads::Dict=Dict())
@defcb truss_bars_volume()
@defcb node_displacement_function(res::Any)

export view_truss_deformation
view_truss_deformation(
  results::Any=nothing,
  visualizer::Backend=autocad;
  factor::Real=100) =
  let disp = node_displacement_function(results),
      b = top_backend()
    with(current_backend, visualizer) do
      for node in b.truss_node_data
        d = disp(node)*factor
        p = node.loc
        truss_node(p+d, family=node.family)
      end
      for bar in b.truss_bar_data
        let (node1, node2) = (bar.node1, bar.node2),
            (p1, p2) = (node1.loc, node2.loc),
            (d1, d2) = (disp(node1)*factor, disp(node2)*factor)
          truss_bar(p1+d1, p2+d2, family=bar.family)
        end
      end
    end
  end

export show_truss_deformation
show_truss_deformation(
    results::Any=nothing,
    visualizer::Backend=autocad;
    node_radius::Real=0.08, bar_radius::Real=0.02, factor::Real=100,
    deformation_name::String="Deformation",
    deformation_color::RGB=rgb(1, 0, 0),
    no_deformation_name::String="No deformation",
    no_deformation_color::RGB=rgb(0, 1, 0)) =
  let disp = node_displacement_function(results),
      b = top_backend()
    with(current_backend, visualizer) do
      delete_all_shapes()
      with(current_layer, create_layer(no_deformation_name, true, no_deformation_color)) do
        for node in b.truss_node_data
          p = node.loc
          sphere(p, node_radius)
        end
        for bar in b.truss_bar_data
          let (node1, node2) = (bar.node1, bar.node2),
              (p1, p2) = (node1.loc, node2.loc)
            cylinder(p1, bar_radius, p2)
          end
        end
      end
      with(current_layer, create_layer(deformation_name, true, deformation_color)) do
        for node in b.truss_node_data
          d = disp(node)*factor
          p = node.loc
          sphere(p+d, node_radius)
        end
        for bar in b.truss_bar_data
          let (node1, node2) = (bar.node1, bar.node2),
              (p1, p2) = (node1.loc, node2.loc),
              (d1, d2) = (disp(node1)*factor, disp(node2)*factor)
            cylinder(p1+d1, bar_radius, p2+d2)
          end
        end
      end
    end
  end

export max_displacement
max_displacement(results, b::Backend=top_backend()) =
  let disp = node_displacement_function(results)
    maximum(map(norm∘disp, b.truss_node_data))
  end

@defcb lighting_analysis()

###################################
# BIM
@defcb all_levels()
@defcb all_walls()
@defcb all_walls_at_level(level)
@defcbs realize_beam_profile(s::BIMShape, profile::Path, cb::Loc, length::Real)
#@defcbs slab(profile, holes, thickness, family)
#@defcbs wall(path, height, l_thickness, r_thickness, family)
#@defcbs panel(bot::Locs, top::Locs, family)

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

#=
The typical family change is:

with(default_XPTO_family, XPTO_family_element(default_XPTO_family(), param=value)) do
  ...
end

We will simplify this pattern by allowing the following syntax:

with(XPTO_family, param=value) do
  ...
end
=#

family_profile(b::Backend, family) =
  family.profile

ext: jl


text:
@params mutable struct SIMPResult{T}
    topology::AbstractVector{T}
    objval::T
    convstate
    nsubproblems::Int
end
Base.show(::IO, ::MIME{Symbol("text/plain")}, ::SIMPResult) = println("TopOpt SIMP result")

function NewSIMPResult(::Type{T}, optimizer, ncells) where {T}
    return SIMPResult(fill(T(NaN), ncells), T(NaN), Nonconvex.ConvergenceState(), 0)
end

"""
The vanilla SIMP algorithm, see [Bendsoe1989](@cite).
"""
@params mutable struct SIMP{T} <: AbstractSIMP
    optimizer
    penalty
    prev_penalty
    solver
    result::SIMPResult{T}
    tracing::Bool
end
Base.show(::IO, ::MIME{Symbol("text/plain")}, ::SIMP) = println("TopOpt SIMP algorithm")

function SIMP(optimizer, solver, p::T; tracing=true) where T
    penalty = getpenalty(solver)
    prev_penalty = deepcopy(penalty)
    setpenalty!(penalty, p)
    ncells = getncells(solver.problem)
    result = NewSIMPResult(T, optimizer, ncells)
    return SIMP(optimizer, penalty, prev_penalty, solver, result, tracing)
end

Utilities.getpenalty(s::AbstractSIMP) = s.penalty
function Utilities.setpenalty!(s::AbstractSIMP, p::Number)
    penalty = s.penalty
    s.prev_penalty = deepcopy(penalty)
    setpenalty!(penalty, p)
    setpenalty!(s.solver, p)
    return s
end

function (s::SIMP{T, TO})(x0 = s.solver.vars) where {T, TO <: Optimizer}
    setpenalty!(s.solver, s.penalty.p)
    mma_results = s.optimizer(x0)
    update_result!(s, mma_results)
    return s.result
end

function (s::SIMP{T, TO})(workspace::Nonconvex.Workspace) where {T, TO <: Optimizer}
    mma_results = s.optimizer(workspace)
    update_result!(s, mma_results)
    return s.result
end

function get_topologies(problem, trace::TopOptTrace)
    @unpack black, white, varind = problem
    @unpack x_hist = trace
    nel = length(black)
    topologies = Vector{Float64}[]
    topology = zeros(T, nel)
    for i in 1:length(x_hist)
        update_topology!(topology, black, white, x_hist[i], varind)
        push!(topologies, copy(topology))
    end
    return topologies
end

function update_result!(s::SIMP{T}, mma_results) where T
    # Postprocessing
    @unpack result, optimizer = s
    @unpack problem = s.solver
    @unpack black, white, varind = problem
    nel = getncells(problem)

    update_topology!(result.topology, black, white, mma_results.minimizer, varind)
    result.objval = mma_results.minimum
    if hasproperty(mma_results, :convstate)
        result.convstate = deepcopy(mma_results.convstate)
    end
    return result
end

function update_topology!(topology, black, white, x, varind)
    for i in 1:length(black)
        if black[i]
            topology[i] = 1
        elseif white[i]
            topology[i] = 0
        else
            topology[i] = x[varind[i]]
        end
    end
    return topology
end


text:
using Dates
using Test

using MarketData

using TimeSeries

struct _TestType
end

@testset "combine" begin

@testset "collapse operations" begin
    @testset "collapse squishes correctly" begin
        @test values(collapse(cl, week, first))[2]    == 97.75
        @test timestamp(collapse(cl, week, first))[2] == Date(2000,1,10)

        @test values(collapse(cl, week, first, last))[2]    == 100.44
        @test timestamp(collapse(cl, week, first, last))[2] == Date(2000,1,10)

        @test values(collapse(cl, month, first))[2]    == 100.25
        @test timestamp(collapse(cl, month, first))[2] == Date(2000,2,1)

        @test values(collapse(ohlc, week, first))[2, :] == [102.0, 102.25, 94.75, 97.75]
        @test timestamp(collapse(ohlc, week, first))[2] == Date(2000,1,10)

        @test values(collapse(ohlc, week, first, last))[2, :] == [100.0, 102.25, 99.38, 100.44]
        @test timestamp(collapse(ohlc, week, first, last))[2] == Date(2000,1,10)

        @test values(collapse(ohlc, month, first))[2, :] == [104.0, 105.0, 100.0, 100.25]
        @test timestamp(collapse(ohlc, month, first))[2] == Date(2000,2,1)
    end

    # https://github.com/JuliaStats/TimeSeries.jl/pull/397
    @testset "Array of String" begin
        A = string.(Char.(rand(97:97+25, size(cl))))
        ts = timestamp(cl)
        ta = TimeArray(ts, A)

        let
            ta = collapse(ta, week, first)

            @test values(ta)[2]    == A[6]
            @test timestamp(ta)[2] == ts[6]
        end
    end
end

@testset "merge works correctly" begin
    cl1  = cl[1:3]
    op1  = cl[2:4]
    aapl = tail(AAPL)
    ba   = tail(BA)

    @testset "takes colnames kwarg correctly" begin
        @test colnames(merge(cl, ohlc[:High, :Low], colnames = [:a, :b, :c])) == [:a, :b, :c]
        @test colnames(merge(cl, op, colnames = [:a, :b]))                    == [:a, :b]
        @test_throws ArgumentError merge(cl, op, colnames = [:a])
        @test_throws ArgumentError merge(cl, op, colnames = [:a, :b, :c])

        for mode ∈ [:inner, :left, :right, :outer]
            @test colnames(merge(cl, ohlc[:High, :Low], method = mode, colnames = [:a, :b, :c])) == [:a, :b, :c]
            @test colnames(merge(cl, op, method = mode, colnames = [:a, :b])) == [:a, :b]
            @test_throws ArgumentError merge(cl, op, method = mode, colnames = [:a])
            @test_throws ArgumentError merge(cl, op, method = mode, colnames = [:a, :b, :c])
        end

        # issue #475
        @test colnames(merge(cl, cl, cl, colnames = [:a, :b, :c])) == [:a, :b, :c]
    end

    @testset "returns correct alignment with Dates and values" begin
        @test values(merge(cl, op))     == values(merge(cl, op, method = :inner))
        @test values(merge(cl,op))[2,1] == values(cl)[2,1]
        @test values(merge(cl,op))[2,2] == values(op)[2,1]
    end

    @testset "aligns with disparate sized objects" begin
        @test values(merge(cl, op[2:5]))[1,1]  == values(cl)[2,1]
        @test values(merge(cl, op[2:5]))[1,2]  == values(op)[2,1]
        @test timestamp(merge(cl, op[2:5]))[1] == Date(2000,1,4)
        @test length(merge(cl, op[2:5]))       == 4

        @test length(merge(cl1, op1, method = :inner))      == 2
        @test values(merge(cl1, op1, method = :inner))[2,1] == values(cl1)[3,1]
        @test values(merge(cl1, op1, method = :inner))[2,2] == values(op1)[2,1]

        @test length(merge(cl1, op1, method = :left))      == 3
        @test values(merge(cl1, op1, method = :left))[2,1] == values(cl1)[2,1]
        @test values(merge(cl1, op1, method = :left))[2,2] == values(op1)[1,1]
        @test isnan(values(merge(cl1, op1, method = :left))[1,2])

        @test length(merge(cl1, op1, method = :right))      == 3
        @test values(merge(cl1, op1, method = :right))[2,1] == values(cl1)[3,1]
        @test values(merge(cl1, op1, method = :right))[2,2] == values(op1)[2,1]
        @test isnan(values(merge(cl1, op1, method = :right))[3,1])

        @test length(merge(cl1, op1, method = :outer))      == 4
        @test values(merge(cl1, op1, method = :outer))[2,1] == values(cl1)[2,1]
        @test values(merge(cl1, op1, method = :outer))[2,2] == values(op1)[1,1]
        @test isnan(values(merge(cl1, op1, method = :outer))[1,2])
        @test isnan(values(merge(cl1, op1, method = :outer))[4,1])
    end

    @testset "column names match the correct values" begin
        @test colnames(merge(cl, op[2:5])) == [:Close, :Open]
        @test colnames(merge(op[2:5], cl)) == [:Open, :Close]

        @test colnames(merge(cl, op[2:5], method = :inner)) == [:Close, :Open]
        @test colnames(merge(op[2:5], cl, method = :inner)) == [:Open, :Close]

        @test colnames(merge(cl, op[2:5], method = :left))  == [:Close, :Open]
        @test colnames(merge(op[2:5], cl, method = :left))  == [:Open, :Close]

        @test colnames(merge(cl, op[2:5], method = :right)) == [:Close, :Open]
        @test colnames(merge(op[2:5], cl, method = :right)) == [:Open, :Close]

        @test colnames(merge(cl, op[2:5], method = :outer)) == [:Close, :Open]
        @test colnames(merge(op[2:5], cl, method = :outer)) == [:Open, :Close]
    end

    @testset "unknown method" begin
        @test_throws ArgumentError merge(cl, op, method = :unknown)
    end

    @testset "custom missing values" begin
        ts1 = TimeArray([Date(2018, 1, 1), Date(2018, 1, 2)], [1, 2])
        ts2 = TimeArray([Date(2018, 1, 2), Date(2018, 1, 3)], [3, 4])

        m1 = merge(ts1, ts2, method =  :left, padvalue = 0)
        @test timestamp(m1) == [Date(2018, 1, 1), Date(2018, 1, 2)]
        @test values(m1)    == [1 0; 2 3]

        m2 = merge(ts1, ts2, method = :right, padvalue = 0)
        @test timestamp(m2) == [Date(2018, 1, 2), Date(2018, 1, 3)]
        @test values(m2)    == [2 3; 0 4]

        m3 = merge(ts1, ts2, method = :outer, padvalue = 0)
        @test timestamp(m3) == [Date(2018, 1, 1), Date(2018, 1, 2), Date(2018, 1, 3)]
        @test values(m3)    == [1 0; 2 3; 0 4]
    end

    @testset "vararg input" begin
        ta1 = merge(op, ohlc, cl)
        ta2 = merge(merge(op, ohlc), cl)

        @test timestamp(ta1) == timestamp(ta2)
        @test values(ta1)    == values(ta2)

        # test keyword arguments passing
        ta1 = merge(op, ohlc, cl, method = :outer)
        ta2 = merge(merge(op, ohlc, method = :outer), cl, method = :outer)

        @test timestamp(ta1) == timestamp(ta2)
        @test values(ta1)    == values(ta2)
    end
end

@testset "hcat" begin
    let ta = [cl op cl]
        @test meta(ta)             == meta(cl)
        @test length(colnames(ta)) == 3
        @test colnames(ta)[1]      != colnames(ta)[3]
        @test colnames(ta)[1]      == :Close
        @test colnames(ta)[2]      == :Open
        @test values(ta)           == [values(cl) values(op) values(cl)]
    end

    let ta = [cl ohlc]
        @test meta(ta)             == meta(ta)
        @test length(colnames(ta)) == 5
        @test colnames(ta)[1]      == :Close
        @test colnames(ta)[1]      != colnames(ta)[end]
        @test values(ta)           == [values(cl) values(ohlc)]
    end

    @test_throws DimensionMismatch [cl[1:3] cl]
end

@testset "vcat works correctly" begin
    @testset "concatenates time series correctly in 1D" begin
        a = TimeArray([Date(2015, 10, 01), Date(2015, 11, 01)], [15, 16], [:Number])
        b = TimeArray([Date(2015, 12, 01)], [17], [:Number])
        c = vcat(a, b)

        @test length(c)   == (length(a) + length(b))
        @test colnames(c) == colnames(a)
        @test colnames(c) == colnames(b)
        @test values(c)   == [15, 16, 17]
    end

    @testset "concatenates time series correctly in 2D" begin
        a = TimeArray([Date(2015, 09, 01), Date(2015, 10, 01), Date(2015, 11, 01)], [[15 16]; [17 18]; [19 20]], [:Number1, :Number2])
        b = TimeArray([Date(2015, 12, 01)], [18 18], [:Number1, :Number2])
        c = vcat(a, b)

        @test length(c)   == length(a) + length(b)
        @test colnames(c) == colnames(a)
        @test colnames(c) == colnames(b)
        @test values(c)   == [[15 16]; [17 18]; [19 20]; [18 18]]
    end

    @testset "rejects when column names do not match" begin
        a = TimeArray([Date(2015, 10, 01), Date(2015, 11, 01)], [15, 16], [:foo])
        b = TimeArray([Date(2015, 12, 01)], [17], [:bar])

        @test_throws ArgumentError vcat(a, b)
    end

    @testset "rejects when metas do not match" begin
        a = TimeArray([Date(2015, 10, 01), Date(2015, 11, 01)], [15, 16], [:A], :FirstMeta)
        b = TimeArray([Date(2015, 12, 01)], [17], [:A], :SecondMeta)

        @test_throws ArgumentError vcat(a, b)
    end

    @testset "duplicated timestamps order" begin
        a = TimeArray([Date(2015, 10, 1), Date(2015, 10, 2), Date(2015, 11, 1)], [15, 16, 17])
        b = TimeArray([Date(2015, 10, 2), Date(2015, 11, 1)], [18, 19])

        ts = [Date(2015, 10, 1),
              Date(2015, 10, 2), Date(2015, 10, 2),
              Date(2015, 11, 1), Date(2015, 11, 1)]
        ta = vcat(a, b)
        @test timestamp(ta) == ts
        @test values(ta)    == [15, 16, 18, 17, 19]
    end

    @testset "still works when dates are mixed" begin
        a = TimeArray([Date(2015, 10, 01), Date(2015, 12, 01)], [15, 17])
        b = TimeArray([Date(2015, 11, 01)], [16])
        c = vcat(a, b)

        @test length(c)   == length(a) + length(b)
        @test colnames(c) == colnames(a)
        @test colnames(c) == colnames(b)
        @test values(c)   == [15, 16, 17]
        @test issorted(timestamp(c))
    end
end

@testset "map works correctly" begin
    @testset "works on both time stamps and 1D values" begin
        a = TimeArray([Date(2015, 10, 01), Date(2015, 11, 01)], [15, 16], [:A], :Something)
        b = map((timestamp, values) -> (timestamp + Dates.Year(1), values .- 1), a)

        @test length(b)                   == length(a)
        @test colnames(b)                 == colnames(a)
        @test Dates.year(timestamp(b)[1]) == Dates.year(timestamp(a)[1]) + 1
        @test values(b)[1]                == values(a)[1] - 1
        @test meta(b)                     == meta(a)
    end

    @testset "works on both time stamps and 2D values" begin
        a = TimeArray([Date(2015, 09, 01), Date(2015, 10, 01), Date(2015, 11, 01)], [[15 16]; [17 18]; [19 20]], [:A, :B])
        b = map((timestamp, values) -> (timestamp + Dates.Year(1), [values[1] + 2, values[2] - 1]), a)

        @test length(b)                   == length(a)
        @test colnames(b)                 == colnames(a)
        @test Dates.year(timestamp(b)[1]) == Dates.year(timestamp(a)[1]) + 1
        @test values(b)[1, 1]             == values(a)[1, 1] + 2
        @test values(b)[1, 2]             == values(a)[1, 2] - 1
    end

    @testset "works with order of elements that varies after modifications" begin
        a = TimeArray([Date(2015, 10, 01), Date(2015, 12, 01)], [15, 16], [:A])
        b = map((timestamp, values) -> (timestamp + Dates.Year((timestamp >= Date(2015, 11, 01)) ? -1 : 1), values), a)

        @test length(b) == length(a)
        @test issorted(timestamp(b))
    end

    @testset "map callable object" begin
        (::_TestType)(ts, x) = (ts, x + 42)

        ta = map(_TestType(), cl)
        @test timestamp(ta) == timestamp(cl)
        @test values(ta)    == values(cl) .+ 42
        @test meta(ta)      == meta(cl)
    end
end

end  # @testset "combine"


text:
module BiochemicalLibrary

using ModelingToolkit
using OrderedCollections

export biochemical_library

const _R = 8.314
const _T = 310.0

@parameters t
D = Differential(t)

@variables E[1:2](t) F[1:2](t)

# Chemical species (:Ce)
@parameters k R T
@variables q(t)
Ce_dict = Dict(
    :description => "Chemical species",
    :numports => 1,
    :parameters => OrderedDict(
        k => "Biochemical Constant; exp(mu_0/RT)/V",
        R => "Universal Gas Constant",
        T => "Temperature"
    ),
    :states => OrderedDict(
        q => "Molar Quantity"
    ),
    :equations => [
        0 ~ R * T * log(k * q) - E[1],
        D(q) ~ F[1]
    ],
    :defaults => Dict(k => 1.0, R => _R, T => _T)
)

# Normalised chemical species (:ce)
@parameters k
@variables q(t)
ce_dict = Dict(
    :description => "Chemical species (normalised)",
    :numports => 1,
    :parameters => OrderedDict(
        k => "Biochemical Constant; exp(mu_0/RT)/V"
    ),
    :states => OrderedDict(
        q => "Molar Quantity"
    ),
    :equations => [
        0 ~ log(k * q) - E[1],
        D(q) ~ F[1]
    ],
    :defaults => Dict(k => 1.0, q => 0.0)
)

# Chemical reaction (:Re)
@parameters r R T
Re_dict = Dict(
    :description => "Biochemical reaction",
    :numports => 2,
    :parameters => OrderedDict(
        r => "Reaction Rate",
        R => "Universal Gas Constant",
        T => "Temperature"
    ),
    :equations => [
        0 ~ F[1] + F[2],
        0 ~ F[1] - r * (exp(E[1] / R / T) - exp(E[2] / R / T))
    ],
    :defaults => Dict(r => 1.0, R => _R, T => _T)
)

# Normalised chemical reaction (:re)
@parameters r
re_dict = Dict(
    :description => "Biochemical reaction (normalised)",
    :numports => 2,
    :parameters => OrderedDict(
        r => "Reaction Rate"
    ),
    :equations => [
        0 ~ F[1] + F[2],
        0 ~ F[1] - r * (exp(E[1]) - exp(E[2]))
    ],
    :defaults => Dict(r => 1.0)
)

const biochemical_library = Dict(
    :Ce => Ce_dict,
    :ce => ce_dict,
    :Re => Re_dict,
    :re => re_dict
)

end

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
using LaserTypes
using LaserTypes: pochhammer, immutable_cache
using CSV
using UnPack
using Test
using Base.Threads

ω = 0.057
c = 137.03599908330935
λ = 2π*c/ω
w₀ = 75*λ

normalization(m, p, E₀, Nₚₘ) = factorial(p)/pochhammer(abs(m)+1,p)*E₀*Nₚₘ

csv_folder = joinpath(@__DIR__, "Laguerre-Gauss/fields_csv/")
ntrds = nthreads()

@testset "Thread safety" begin
    for (roots, dirs, files) in walkdir(csv_folder)
        for file in files
            i, p, m = tryparse.(Int, (split(file,"_"))[2:4])
            ξx = cos(i*π/4) + 0im; ξy = im*sin(i*π/4);
            s = setup_laser(LaguerreGaussLaser, :atomic; λ, w₀, p, m, ξx, ξy, profile = ConstantProfile())
            @unpack E₀, Nₚₘ = immutable_cache(s)

            current_file = CSV.File(joinpath(roots, file), header = false)
            num_of_rows = length(current_file)
            test_bool = Atomic{Bool}(true)

            current_row = zeros(Int, ntrds)
            @threads for id in 1:ntrds
                while current_row[id] + id ≤ num_of_rows
                    x = current_file[current_row[id] + id][1]
                    y = current_file[current_row[id] + id][2]
                    z = current_file[current_row[id] + id][3]

                    EX = tryparse(ComplexF64, replace(replace(current_file[current_row[id] + id][4],"*I"=>"im"),"*^"=>"e"))
                    EY = tryparse(ComplexF64, replace(replace(current_file[current_row[id] + id][5],"*I"=>"im"),"*^"=>"e"))
                    EZ = tryparse(ComplexF64, replace(replace(current_file[current_row[id] + id][6],"*I"=>"im"),"*^"=>"e"))

                    BX = tryparse(ComplexF64, replace(replace(current_file[current_row[id] + id][7],"*I"=>"im"),"*^"=>"e"))
                    BY = tryparse(ComplexF64, replace(replace(current_file[current_row[id] + id][8],"*I"=>"im"),"*^"=>"e"))
                    BZ = tryparse(ComplexF64, replace(replace(current_file[current_row[id] + id][9],"*I"=>"im"),"*^"=>"e"))

                    juliaE = E([x,y,z], s)
                    juliaB = B([x,y,z], s)

                    atomic_and!(test_bool, juliaE[1] ≈ normalization(s.m, s.p, E₀, Nₚₘ)*EX)
                    atomic_and!(test_bool, juliaE[2] ≈ normalization(s.m, s.p, E₀, Nₚₘ)*EY)
                    atomic_and!(test_bool, juliaE[3] ≈ normalization(s.m, s.p, E₀, Nₚₘ)*EZ)

                    atomic_and!(test_bool, juliaB[1] ≈ normalization(s.m, s.p, E₀, Nₚₘ)*BX)
                    atomic_and!(test_bool, juliaB[2] ≈ normalization(s.m, s.p, E₀, Nₚₘ)*BY)
                    atomic_and!(test_bool, juliaB[3] ≈ normalization(s.m, s.p, E₀, Nₚₘ)*BZ)
                    current_row[id] += ntrds
                end
            end
            @test test_bool.value
        end
    end
end


text:
##################################################################################
# This file is part of ModelBaseEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

export export_model

_check_name(name) = Base.isidentifier(name) ? true : throw(ArgumentError("Model name must be a valid Julia identifier."))

"""
    export_model(model, name, file::IO)
    export_model(model, name, path::String)

Export the model into a module file. The `name` parameter is used for the name
of the module as well as the module file. The module file is created in the
directory specified by the optional third argument.
"""
function export_model(m::Model, name::AbstractString, path::AbstractString=".")
    if !endswith(path, ".jl")
        path = joinpath(path, name * ".jl")
    end
    open(path, "w") do fd
        export_model(m, name, IOContext(fd, :compact => false, :limit => false))
    end
    return nothing
end

function export_model(m::Model, name::AbstractString, fio::IO)
    _check_name(name)
    println(fio, "module ", name)
    println(fio)
    println(fio, "using ModelBaseEcon")
    println(fio)
    println(fio, "const model = Model()")
    println(fio)

    function _print_modified_options(opts, default, prefix)
        for (ok, ov) in pairs(opts)
            dv = getoption(default, ok, :not_a_default)
            if ov isa Options && dv isa Options
                _print_modified_options(ov, dv, prefix * "$ok.")
            elseif dv == :not_a_default || dv != ov
                println(fio, prefix, ok, " = ", ov)
            end
        end
    end

    println(fio, "# options")
    _print_modified_options(m.options, defaultoptions, "model.")
    println(fio)

    println(fio, "# flags")
    for fld in fieldnames(ModelFlags)
        fval = getfield(m.flags, fld)
        if fval != getfield(ModelFlags(), fld)
            println(fio, "model.", fld, " = ", fval)
        end
    end
    println(fio)

    if !isempty(parameters(m))
        println(fio, "@parameters model begin")
        for (n, p) in m.parameters
            println(fio, "    ", n, " = ", p)
        end
        println(fio, "end # parameters")
        println(fio)
    end

    allvars = m.allvars
    if !isempty(allvars)
        println(fio, "@variables model begin")
        has_exog = false
        has_shocks = false
        for v in allvars
            if isexog(v)
                has_exog = true
            elseif isshock(v)
                has_shocks = true
            else
                println(fio, "    ", v)
            end
        end
        println(fio, "end # variables")
        println(fio)
        if has_exog
            println(fio, "@exogenous model begin")
            for v in allvars
                if isexog(v)
                    doc = ifelse(isempty(v.doc), "", v.doc * " ")
                    println(fio, "    ", doc, v.name)
                end
            end
            println(fio, "end # exogenous")
            println(fio)
        end
        if has_shocks
            println(fio, "@shocks model begin")
            for v in allvars
                if isshock(v)
                    println(fio, "    ", v)
                end
            end
            println(fio, "end # shocks")
            println(fio)
        end
    end

    if !isempty(m.autoexogenize)
        println(fio, "@autoexogenize model begin")
        for (k, v) in pairs(m.autoexogenize)
            println(fio, "    ", k, " = ", v)
        end
        println(fio, "end # autoexogenize")
        println(fio)
    end

    alleqns = m.alleqns
    if !isempty(alleqns)
        println(fio, "@equations model begin")
        for eqn in alleqns
            str = sprint(print, eqn, context=fio, sizehint=0)
            str = replace(str, r"(\s*\".*\"\s*)" => s"\1\\n    ")
            println(fio, "    ", unescape_string(str))
        end
        println(fio, "end # equations")
        println(fio)
    end

    println(fio, "@initialize model")

    sd = sstate(m)
    for cons in sd.constraints
        println(fio)
        println(fio, "@steadystate model ", cons)
    end

    println(fio)
    println(fio, "end # module ", name)
    println(fio)

    return nothing
end


text:
using GLVisualize, GeometryTypes, Colors, GLAbstraction

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0,1,20))
end

description = """
Displaying changing unicode text with GLVisualize.
"""

# GLAbstraction.const_lift is defined as
# const_lift(F, args...) = Reactive.map(F, map(Signal, args)...)
# and allows to also use constant arguments without manually wrapping
# them into signals
π_signal = const_lift(*, timesignal, pi)

# Reactive.map takes a signal and transforms it by applying F into a new signal
# it's similar to Base.map, just over the time dimension (like registering a
# callback in a more traditional event system)
s = map(π_signal) do π
    """The quick brown fox jumped over
    some lazy text sample.
    He wasn't really into numbers, but it's
    really important to try out number rendering:
    This number goes from 0 to $(@sprintf("%0.6f", π)) in no time!
    And then back to 0 again... Wow!
    This is real crazy stuff,
    but it gets even more ludicrous:
    ∮ E⋅da = Q,  n → ∞, ∑ f(i) = ∏ g(i),
    ∀x∈ℝ: ⌈x⌉ = −⌊−x⌋, α ∧ ¬β = ¬(¬α ∨ β),
    ℕ ⊆ ℕ₀ ⊂ ℤ ⊂ ℚ ⊂ ℝ ⊂ ℂ,
    ⊥ < a ≠ b ≡ c ≤ d ≪ ⊤ ⇒ (A ⇔ B),
    2H₂ + O₂ ⇌ 2H₂O, R = 4.7 kΩ, ⌀
    I can't even...
    """
end

# _view and visualize it!
_view(visualize(s), window, camera = :orthographic_pixel)

# text won't come out centered by default, so lets center it!
center!(window, :orthographic_pixel)

if !isdefined(:runtests)
    renderloop(window)
end


text:
using ArrayInterface
using EllipsisNotation
using Test

include("basic.jl")
include("more_generic.jl")


text:
# fizzbuzz
for i=1:100
    if mod(i,15)==0     println("FizzBuzz")
    elseif mod(i,5)==0  println("Buzz")
    elseif mod(i,3)==0  println("Fizz")
    else println(i)
    end
end

max_stars_repo_path: Module 1/Chapter01/fizzbuzz.jl

max_issues_repo_path: Module 1/Chapter01/fizzbuzz.jl

max_forks_repo_path: Module 1/Chapter01/fizzbuzz.jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule isoband_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("isoband")
JLLWrappers.@generate_main_file("isoband", UUID("9a68df92-36a6-505f-a73e-abb412b6bfb4"))
end  # module isoband_jll


text:
include(joinpath(@__DIR__, "../../src/FLS.jl"))

"""
Cartesian Product
"""

pairschema = @acset LabeledGraph begin
    V = 2
    E = 2
    vlabel = [:s, :s2]
    elabel = [:p1, :p2]
    src = [2,2]
    tgt = [1,1]
end

paird= @acset LabeledGraph begin
    V = 2
    vlabel = [:s, :s]
end

pairs = FLS(:pairs, pairschema, [
    Cone(lconed, :s2, [1=>:p1,2=>:p2])], [])


text:
"""
ULMFiT - Binary Sentiment Analyzer

    BinSentimentClassifier()

This is a binary sentiment classifier developed after
fine-tuning the ULMFiT language model on IMDB movie reviews dataset.

# Usage:

julia> sc = BinSentimentClassifier()

julia> doc = StringDocument("this classifier is great")

julia> sc(doc)
"positive"

"""

struct BinSentimentClassifier
    vocab::Vector
    rnn_layers::Flux.Chain
    linear_layers::Flux.Chain
end

function BinSentimentClassifier()
    BSON.@load datadep"ULMFiT Sentiment Classifier/umlfit_sentiment_en.bson" weights
    vocab_sz, em_sz = size(weights[1])
    hid_lstm_sz = 1150
    out_lstm_sz = em_sz
    clsfr_hid_sz = 50
    clsfr_out_sz = 2
    vocab = (string.(readdlm("vocabs/sc_vocab.csv", ',')))[:, 1]
    sc = SentimentClassifier(
        vocab,
        Chain(
            DroppedEmbeddings(vocab_sz, em_sz),
            LSTM(em_sz, hid_lstm_sz),
            LSTM(hid_lstm_sz, hid_lstm_sz),
            LSTM(hid_lstm_sz, out_lstm_sz)
        ),
        Chain(
            PooledDense(out_lstm_sz, clsfr_hid_sz),
            BatchNorm(clsfr_hid_sz, relu),
            Dense(clsfr_hid_sz, clsfr_out_sz, sigmoid),
            BatchNorm(clsfr_out_sz),
            softmax
        )
    )
    Flux.loadparams!(sc, weights)
    sc = mapleaves(Tracker.data, sc)
    Flux.testmode!(sc)
    return sc
end

Flux.@treelike BinSentimentClassifier

function (sc::BinSentimentClassifier)(x::TokenDocument)
    remove_case!(x)
    idxs = map(w -> indices([w], sc.vocab, "_unk_"), tokens(x))
    h = sc.rnn_layers.(idxs)
    h = sc.linear_layers(h)
    Flux.reset!(sc.rnn_layers)
    return argmax(h)[1] == 1 ? "positive" : "negative"
end


text:
using Distributed

# Number of matrices in the sum
nmats = 2
addprocs(nmats)

@everywhere begin
  using LinearAlgebra
  using KrylovKit

  BLAS.set_num_threads(1)

  include(joinpath("src", "DistributedSums.jl"))
  include(joinpath("src", "MatrixMaps.jl"))

  # Matrix dimension
  n = 100

  # Starting vector that we will apply
  # the matrix map to
  v0 = randn(n)
end

# Make a distributed sum of MatrixMaps, which is
# just a Matrix with special multiplication
# notation that can work with KrylovKit
A = DistributedSum(_ -> MatrixMap(randn(n, n)), nmats)

Av0 = A(v0)

@show norm(Av0 - sum(fetch.(A))(v0))

# Solve eigenvector equation
λ⃗, v⃗ = eigsolve(A, v0, 1, :LM)
λ = λ⃗[1]
v = v⃗[1]

@show norm(A(v) - λ * v)

max_stars_repo_name: mtfishman/ITensorParallel


text:
using FinancialSeries

tt = readdlm(Pkg.dir("FinancialSeries/data/tickdata.csv"), '\t')

t  = convert(Array{Float64}, tt[2:end, 2:end])
cn = convert(Array{ASCIIString}, tt[1,2:end])
ct = convert(Array{ASCIIString}, [c for c in cn])
dt = convert(Vector{ASCIIString}, tt[2:end,1])
da = DateTime{ISOCalendar, UTC}[parsedatetime1(d) for d in dt];
td = TickData(da, t, ct, Stock(Ticker("Foo")))


text:
#=
Hackerrank problem - Geometric Distribution
https://www.hackerrank.com/challenges/s10-geometric-distribution-1/problem
https://www.hackerrank.com/challenges/s10-geometric-distribution-2/problem

Date: 2020-07-30
Author: Eric Nam
=#

using Printf

function dist_geomteric(p, n)
    (1. - p) ^ (n - 1) * p
end

# Part 1
function day4_1()
    # STDIN
    nu, denom = parse.(Int, split(readline(stdin)))
    nth = parse(Int, readline(stdin))

    result = dist_geomteric(nu / denom, nth)
    println("$(@sprintf("%.3f", result))")
end

# day4_1()

# Part 2
function day4_2()
    # STDIN
    nu, denom = parse.(Int, split(readline(stdin)))
    nth = parse(Int, readline(stdin))

    result = sum(x -> dist_geomteric(nu / denom, x), 1:5)
    println("$(@sprintf("%.3f", result))")
end

day4_2()


text:
using BinaryBuilder

name = "mongo-c-driver-builder"
version = v"1.14.0"

# Collection of sources required to build mongo-c-driver-builder
sources = [
    "https://github.com/mongodb/mongo-c-driver/releases/download/1.14.0/mongo-c-driver-1.14.0.tar.gz" =>
    "ebe9694f7fa6477e594f19507877bbaa0b72747682541cf0cf9a6c29187e97e8",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd mongo-c-driver-1.14.0/
mkdir cmake-build
cd cmake-build/
cmake -DENABLE_AUTOMATIC_INIT_AND_CLEANUP=OFF -DCMAKE_BUILD_TYPE=Release -DENABLE_SSL=AUTO -DENABLE_STATIC=ON -DENABLE_MONGOC=ON -DENABLE_BSON=ON -DENABLE_ZLIB=BUNDLED -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain ..
make -j${nproc}
make install
exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    MacOS(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libmongoc", :libmongoc),
    LibraryProduct(prefix, "libbson", :libbson)
]

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

]

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

max_stars_repo_name: felipenoris/mongo-c-driver-builder

max_issues_repo_name: felipenoris/mongo-c-driver-builder

max_forks_repo_name: felipenoris/mongo-c-driver-builder


text:
# module CommonTools
function filesdirs(path)
    cd = readdir(path)
    filter!(x->x[1] != '.',cd)
    dirs = isdir.(cd);
    return cd[!dirs] , cd[dirs]
end
export filesdirs

function includeallafiles(expr)
    myfullname = expr;
    localname = basename(myfullname);
    path = dirname(myfullname);
    (files,dirs) = filesdirs(path);
    for f in files
        if f == localname;
            continue;
        end;
        include(joinpath(path,f))
    end
end
export includeallafiles

ishidden(fname::String) = match(r"\\\.",fname) != nothing || match(r"^\.",fname) != nothing
ishidden(fnames::Vector{String}) = map(x->ishidden(x),fnames)
export ishidden

function listdirs_rec(path::String)
    V = Vector{String}();
    for (root, dirs, files) in walkdir(path)
        for dir in dirs
            if ishidden(root) || ishidden(dir) continue; end
            push!(V,joinpath(root, dir)) # path to directories
        end
    end
    V
end
export listdirs_rec

"return (path,file,ext)"
function fileparts(fullname)
    path = dirname(fullname)
    file = basename(fullname)
    m = match(r"([^\.]+)(\.(.+))?",file);
    f(x) = x == nothing ? "" : x;
    file = f(m.captures[1])
    ext = f(m.captures[3])
    return (path,file,ext)
end
export fileparts

saveitem(fname,item) = open(fname,"w") do f
    serialize(f,item)
end
export saveitem

loaditem(fname) = open(fname,"r") do f
    deserialize(f)
end
export loaditem

u8(x) = convert(Matrix{UInt8},round.(x))
export u8

using FileIO
cache = Dict{String,Any}()
function cached_load(fname)
    global cache
    if fname in keys(cache)
        return cache[fname]
    else
        return cache[fname] = load(fname)
    end
end
export cached_load

""" forward abspath """
function fabspath(path...)
    return replace(abspath(path...),"\\", "/")
end
export fabspath


text:
module Runtests

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    # Desired:
    #   using Test
    #
    # In order that Test need not to be defined in the user environment, it is included via ModiaMath:
    using ModiaMath.Test
end

@testset "Test ModiaMath" begin
    include(joinpath("result", "_includes.jl"))
    include(joinpath("variables", "_includes.jl"))
    include(joinpath("frames", "_includes.jl"))
    include(joinpath("simulation", "_includes.jl"))
end

end


text:
using MCParticles
using Test

@testset "MCParticles.jl constructor: DynamicParticle" begin
    # Constructor tests
    f(x,t) = -sum(x) + t
    dyn_p1_1 = DynamicParticle([1.0,1.0], log(1.0), f, 1.0, f([1.0,1.0],1.0))
    dyn_p1_2 = Particle([1.0,1.0], 1.0, f, 1.0)
    @test pequal(dyn_p1_1,dyn_p1_2)

    g(x) = -sum(x)
    dyn_p2_1 = DynamicParticle([1.0 1.0], log(1.0), g, 1, g([1.0 1.0]))
    dyn_p2_2 = Particle([1.0 1.0], 1.0, g, 1)
    @test pequal(dyn_p2_1,dyn_p2_2)
end

@testset "MCParticles.jl constructor: StaticParticle" begin
    # Constructor tests
    stat_p1_1 = StaticParticle([1.0,1.0], log(1.0))
    stat_p1_2 = Particle([1.0,1.0], 1.0)
    @test pequal(stat_p1_1, stat_p1_2)
end

@testset "MCParticles.jl helper functions" begin
    # MCParticles.normalise
    lw = [ -8, -4, -2, 0, 2, 4, 8] .* log(10.0)
    lsumw = log( sum(10.0 .^ [ -8, -4, -2, 0, 2, 4, 8]) )
    nlw = lw .- lsumw
    normalise_nlw = MCParticles.normalise(lw)
    @test normalise_nlw ≈ nlw
end

@testset "MCParticles.jl accessor functions" begin
    #dynamic particles
    f(x,t) = -sum(x) + t
    dyn_p1 = Particle([1.0,1.0], 1.0, f, 1.0)
    dyn_p2 = Particle([1.0,2.0], 0.5, f, 1.0)
    arr_dyn_p = [dyn_p1, dyn_p1, dyn_p2, dyn_p2]
    arr_w = [1.0, 1.0, 0.5, 0.5]
    @test weights(arr_dyn_p, logscale = false, normalscale = false) ≈ arr_w
    @test weights(arr_dyn_p, logscale = true, normalscale = false) ≈ log.(arr_w)
    @test weights(arr_dyn_p, logscale = false, normalscale = true) ≈ arr_w/3.0
    # static particles
    stat_p1 = Particle([1.0,1.0], 1.0)
    stat_p2 = Particle([1.0,1.0], 0.5)
    arr_stat_p = [stat_p1, stat_p1, stat_p2, stat_p2]
    @test weights(arr_stat_p, logscale = false, normalscale = false) ≈ arr_w
    @test weights(arr_stat_p, logscale = true, normalscale = false) ≈ log.(arr_w)
    @test weights(arr_stat_p, logscale = false, normalscale = true) ≈ arr_w/3.0
end

@testset "MCParticles.jl mutate functions" begin
    #DynamicParticle
    f(x,t) = -sum(x) + t
    dyn_p1 = Particle([1.0,1.0], 1.0, f, 1.0)
    dyn_p2 = Particle([1.0,2.0], 1.0, f, 2.0)
    dyn_p1 = mutate(dyn_p1, [1.0,2.0], 2.0)
    @test pequal(dyn_p1, dyn_p2)
    #StaticParticle
    stat_p1 = Particle([1.0,1.0], 1.0)
    stat_p2 = Particle([1.0,2.0], 1.0)
    stat_p1 = mutate(stat_p1, [1.0,2.0])
    @test pequal(stat_p1, stat_p2)
end

@testset "MCParticles.jl reweight functions" begin
    #DynamicParticle
    # f(x,t) = -sum(x) + t
    # dyn_p1 = Particle([1.0,1.0], 1.0, f, 1.0)
    # dyn_p2 = Particle([1.0,2.0], 1.0, f, 2.0)
    # dyn_p1 = mutate(dyn_p1, [1.0,2.0], 2.0)
    # @test pequal(dyn_p1, dyn_p2)
    #StaticParticle
    stat_p1 = Particle([1.0,1.0], 1.0)
    stat_p2 = Particle([1.0,1.0], 2.0)
    stat_p1 = reweight(stat_p1, 2.0)
    @test pequal(stat_p1, stat_p2)
end


text:
module test_loltools_leaguev4

using Test
using LOLTools.LeagueV4

include("helpers.jl")

struct LeagueController <: ApplicationController
    conn::Conn
end

function challengerleagues(c::LeagueController)
    # @info :c c.params.entries
    render(JSON, (leagueId="", tier="CHALLENGER", entries=[(summonerId = "", summonerName = "", leaguePoints = 0, rank = "I", wins = 0, losses = 0, veteran = false, inactive = false, freshBlood = true, hotStreak = false)], queue="RANKED_SOLO_5x5", name="Nautilus's Outriders"))
end

routes() do
    get("/lol/league/v4/challengerleagues/by-queue/:queue", LeagueController, challengerleagues)
end

api_key = ""
platform = "kr"
LeagueV4.challengerleagues(api_key, platform, "RANKED_SOLO_5x5"; action=mock_action)

end # module test_loltools_leaguev4


text:
include("acs.jl")
include("web.jl")
include("aps.jl")


text:
<extend src="layout.jl">
    <block name="content">

Angel

    </block>
</extend>

max_stars_repo_name: misterdev/angel

max_issues_repo_name: misterdev/angel

max_forks_repo_name: misterdev/angel


text:
"""
    FEMMHeatDiffModule

Module for operations on interiors of domains to construct system
matrices and system vectors for linear heat conduction/diffusion.
"""
module FEMMHeatDiffModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: AbstractFESet, nodesperelem, manifdim, gradN!
import FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobianvolume
import FinEtools.CSysModule: CSys, updatecsmat!
import FinEtools.FieldModule: ndofs, gatherdofnums!, gatherfixedvalues_asvec!, gathervalues_asvec!, gathervalues_asmat!
import FinEtools.NodalFieldModule: NodalField
import FinEtools.ElementalFieldModule: ElementalField
import FinEtools.AssemblyModule: AbstractSysvecAssembler, AbstractSysmatAssembler, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!, makevector!, SysvecAssembler
import FinEtools.ForceIntensityModule: ForceIntensity
import FinEtools.FEMMBaseModule: AbstractFEMM, inspectintegpoints
import FinEtools.MatrixUtilityModule: add_gkgt_ut_only!, complete_lt!, locjac!, add_nnt_ut_only!, mulCAtB!, mulCAB!
import LinearAlgebra: norm, dot
import FinEtoolsHeatDiff.MatHeatDiffModule: MatHeatDiff, tangentmoduli!, update!

"""
    FEMMHeatDiff{S<:AbstractFESet, F<:Function, M<:MatHeatDiff} <: AbstractFEMM

Type for heat diffusion finite element modeling machine.
"""
mutable struct FEMMHeatDiff{S<:AbstractFESet, F<:Function, M<:MatHeatDiff} <: AbstractFEMM
    integdomain::IntegDomain{S, F} # geometry data
    mcsys::CSys # updater of the material orientation matrix
    material::M # material object
end

"""
    FEMMHeatDiff(integdomain::IntegDomain{S, F}, material::M) where {S<:AbstractFESet, F<:Function, M<:MatHeatDiff}

Construct with the default orientation matrix (identity).
"""
function FEMMHeatDiff(integdomain::IntegDomain{S, F}, material::M) where {S<:AbstractFESet, F<:Function, M<:MatHeatDiff}
    return FEMMHeatDiff(integdomain, CSys(manifdim(integdomain.fes)), material)
end

function  _buffers1(self::FEMMHeatDiff, geom::NodalField{FFlt}, temp::NodalField{FFlt})
    # Constants
    fes = self.integdomain.fes
    nfes = count(fes); # number of finite elements in the set
    ndn = ndofs(temp); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);   # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    Kedim = ndn*nne;      # dimension of the element matrix
    ecoords = fill(zero(FFlt), nne, ndofs(geom)); # array of Element coordinates
    elmat = fill(zero(FFlt), Kedim, Kedim); # buffer
    elvec = fill(zero(FFlt), Kedim); # buffer
    elvecfix = fill(zero(FFlt), Kedim); # buffer
    dofnums = fill(zero(FInt), Kedim); # buffer
    loc = fill(zero(FFlt), 1, sdim); # buffer
    J = fill(zero(FFlt), sdim, mdim); # buffer
    RmTJ = fill(zero(FFlt), mdim, mdim); # buffer
    gradN = fill(zero(FFlt), nne, mdim); # buffer
    kappa_bar = fill(zero(FFlt), mdim, mdim); # buffer
    kappa_bargradNT = fill(zero(FFlt), mdim, nne); # buffer
    return ecoords, dofnums, loc, J, RmTJ, gradN, kappa_bar, kappa_bargradNT, elmat, elvec, elvecfix
end

function  _buffers2(self::FEMMHeatDiff, geom::NodalField{FFlt}, temp::NodalField{FFlt})
    # Constants
    fes = self.integdomain.fes
    nfes = count(fes); # number of finite elements in the set
    ndn = ndofs(temp); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);   # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    Kedim = ndn*nne;      # dimension of the element matrix
    ecoords = fill(zero(FFlt), nne, ndofs(geom)); # array of Element coordinates
    elmat = fill(zero(FFlt), Kedim, Kedim); # buffer
    elvec = fill(zero(FFlt), Kedim); # buffer
    elvecfix = fill(zero(FFlt), Kedim); # buffer
    dofnums = fill(zero(FInt), Kedim); # buffer
    loc = fill(zero(FFlt), 1, sdim); # buffer
    J = fill(zero(FFlt), sdim, mdim); # buffer
    # RmTJ = fill(zero(FFlt), mdim, mdim); # buffer
    # gradN = fill(zero(FFlt), nne, mdim); # buffer
    # kappa_bar = fill(zero(FFlt), mdim, mdim); # buffer
    # kappa_bargradNT = fill(zero(FFlt), mdim, nne); # buffer
    return ecoords, dofnums, loc, J, elmat, elvec, elvecfix
end

"""
    conductivity(self::FEMMHeatDiff,  assembler::A, geom::NodalField{FFlt},  temp::NodalField{FFlt}) where {A<:AbstractSysmatAssembler}

Compute the conductivity matrix.

# Arguments
- `self` = model machine,
- `assembler` = matrix assembler
- `geom` = geometry field,
- `temp` = temperature field
"""
function conductivity(self::FEMMHeatDiff,  assembler::A, geom::NodalField{FFlt},  temp::NodalField{FFlt}) where {A<:AbstractSysmatAssembler}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    # Prepare assembler and buffers
    ecoords, dofnums, loc, J, RmTJ, gradN, kappa_bar, kappa_bargradNT, elmat = _buffers1(self, geom, temp)
    # Thermal conductivity matrix is in local  material coordinates.
    kappa_bar = tangentmoduli!(self.material, kappa_bar)
    startassembly!(assembler, size(elmat,1), size(elmat,2), count(fes), temp.nfreedofs, temp.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom, ecoords, fes.conn[i]);
        fill!(elmat,  0.0); # Initialize element matrix
        for j=1:npts # Loop over quadrature points
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            updatecsmat!(self.mcsys, loc, J, fes.label[i]);
            mulCAtB!(RmTJ,  self.mcsys.csmat,  J); # local Jacobian matrix
            gradN!(fes, gradN, gradNparams[j], RmTJ);
            # Add the product gradN*kappa_bar*gradNT*(Jac*w[j])
            add_gkgt_ut_only!(elmat, gradN, (Jac*w[j]), kappa_bar, kappa_bargradNT)
        end # Loop over quadrature points
        complete_lt!(elmat)
        gatherdofnums!(temp, dofnums, fes.conn[i]);# retrieve degrees of freedom
        assemble!(assembler, elmat, dofnums, dofnums);# assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function conductivity(self::FEMMHeatDiff, geom::NodalField{FFlt},  temp::NodalField{FFlt})
    assembler = SysmatAssemblerSparseSymm();
    return conductivity(self, assembler, geom, temp);
end

"""
    nzebcloadsconductivity(self::FEMMHeatDiff, assembler::A,  geom::NodalField{FFlt},  temp::NodalField{FFlt}) where {A<:AbstractSysvecAssembler}

Compute load vector for nonzero EBC of prescribed temperature.

# Arguments
- `self` = model machine,
- `assembler` = matrix assembler
- `geom` = geometry field,
- `temp` = temperature field
"""
function nzebcloadsconductivity(self::FEMMHeatDiff, assembler::A,  geom::NodalField{FFlt},  temp::NodalField{FFlt}) where {A<:AbstractSysvecAssembler}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    # Prepare assembler and buffers
    ecoords, dofnums, loc, J, RmTJ, gradN, kappa_bar, kappa_bargradNT, elmat, elvec, elvecfix = _buffers1(self, geom, temp)
    # Thermal conductivity matrix is in local  material coordinates.
    kappa_bar = tangentmoduli!(self.material, kappa_bar)
    startassembly!(assembler,  temp.nfreedofs);
    # Now loop over all finite elements in the set
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom, ecoords, fes.conn[i]);
        gatherfixedvalues_asvec!(temp, elvecfix, fes.conn[i]);# retrieve element coordinates
        if norm(elvecfix, Inf) != 0. # Is the load nonzero?
            fill!(elmat,  0.0);
            for j=1:npts # Loop over quadrature points
                locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
                Jac = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
                updatecsmat!(self.mcsys, loc, J, fes.label[i]);
                mulCAtB!(RmTJ,  self.mcsys.csmat,  J); # local Jacobian matrix
                gradN!(fes, gradN, gradNparams[j], RmTJ);
                # Add the product gradN*kappa_bar*gradNT*(Jac*w[j])
                add_gkgt_ut_only!(elmat, gradN, (Jac*w[j]), kappa_bar, kappa_bargradNT)
            end # Loop over quadrature points
            complete_lt!(elmat)
            mulCAB!(elvec, elmat, elvecfix) # compute  the load vector
            gatherdofnums!(temp, dofnums, fes.conn[i]); # retrieve degrees of freedom
            assemble!(assembler,  -elvec,  dofnums); # assemble element load vector
        end
    end
    return makevector!(assembler);
end

function nzebcloadsconductivity(self::FEMMHeatDiff,  geom::NodalField{FFlt},   temp::NodalField{FFlt})
    assembler = SysvecAssembler()
    return  nzebcloadsconductivity(self, assembler, geom, temp);
end

"""
    energy(self::FEMMHeatDiff, geom::NodalField{FFlt},  temp::NodalField{FFlt})

Compute the "energy" integral over the interior domain.

The "energy" density is the dot product of the gradient of temperature
and the heat flux.

# Arguments
- `self` = model machine,
- `geom` = geometry field,
- `temp` = temperature field
"""
function energy(self::FEMMHeatDiff, geom::NodalField{FFlt},  temp::NodalField{FFlt})
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    # Prepare assembler and buffers
    ecoords, dofnums, loc, J, RmTJ, gradN, kappa_bar, kappa_bargradNT, elmat, elvec, elvecfix = _buffers1(self, geom, temp)
    # Thermal conductivity matrix is in local  material coordinates.
    kappa_bar = tangentmoduli!(self.material, kappa_bar)
    gradT = fill(0.0, 1, size(gradN, 2))
    fluxT = reshape(deepcopy(gradT), length(gradT), 1)
    energy = 0.0
    # Now loop over all finite elements in the set
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom, ecoords, fes.conn[i]);
        gathervalues_asvec!(temp, elvec, fes.conn[i]);# retrieve element coordinates
        for j=1:npts # Loop over quadrature points
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            updatecsmat!(self.mcsys, loc, J, fes.label[i]);
            mulCAtB!(RmTJ,  self.mcsys.csmat,  J); # local Jacobian matrix
            gradN!(fes, gradN, gradNparams[j], RmTJ);
            mulCAtB!(gradT, reshape(elvec, length(elvec), 1), gradN)
           	mulCAB!(fluxT, kappa_bar, reshape(gradT, length(gradT), 1))
            energy += dot(vec(gradT), vec(fluxT)) * (Jac*w[j])
        end # Loop over quadrature points
    end
    return energy;
end

"""
    inspectintegpoints(self::FEMMHeatDiff, geom::NodalField{FFlt}, u::NodalField{T}, temp::NodalField{FFlt}, felist::FIntVec, inspector::F, idat, quantity=:heatflux; context...) where {T<:Number, F<:Function}

Inspect integration point quantities.

# Arguments
- `geom` - reference geometry field
- `u` - displacement field (ignored)
- `temp` - temperature field
- `felist` - indexes of the finite elements that are to be inspected:
    The fes to be included are: `fes[felist]`.
- `context`    - structure: see the update!() method of the material.
- `inspector` - function with the signature
        `idat = inspector(idat, j, conn, x, out, loc);`
   where
    `idat` - a structure or an array that the inspector may
           use to maintain some state,  for instance heat flux, `j` is the
          element number, `conn` is the element connectivity, `out` is the
          output of the `update!()` method,  `loc` is the location of the
          integration point in the *reference* configuration.
# Output
The updated inspector data is returned.
"""
function inspectintegpoints(self::FEMMHeatDiff, geom::NodalField{FFlt}, u::NodalField{T}, temp::NodalField{FFlt}, felist::FIntVec, inspector::F, idat, quantity=:heatflux; context...) where {T<:Number, F<:Function}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    ecoords, dofnums, loc, J, RmTJ, gradN, kappa_bar, kappa_bargradNT, elmat, elvec, elvecfix = _buffers1(self, geom, temp)
    # Thermal conductivity matrix is in local  material coordinates.
    kappa_bar = tangentmoduli!(self.material, kappa_bar)
    # Sort out  the output requirements
    outputcsys = self.mcsys; # default: report the vector quantities in the material coord system
    for apair in pairs(context)
        sy, val = apair
        if sy == :outputcsys
            outputcsys = val
        end
    end
    t= 0.0
    dt = 0.0
    Te = fill(zero(FFlt), nodesperelem(fes)) # nodal temperatures -- buffer
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    qpgradT = fill(zero(FFlt), 1, sdim); # Temperature gradient -- buffer
    qpflux = fill(zero(FFlt), sdim); # thermal strain -- buffer
    out1 = fill(zero(FFlt), sdim); # output -- buffer
    out =  fill(zero(FFlt), sdim);# output -- buffer
    # Loop over  all the elements and all the quadrature points within them
    for ilist = 1:length(felist) # Loop over elements
        i = felist[ilist];
        gathervalues_asmat!(geom, ecoords, fes.conn[i]);
        gathervalues_asvec!(temp, Te, fes.conn[i]);# retrieve element temperatures
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            updatecsmat!(self.mcsys, loc, J, fes.label[i]);
            mulCAtB!(RmTJ,  self.mcsys.csmat,  J); # local Jacobian matrix
            gradN!(fes, gradN, gradNparams[j], RmTJ);
            # Quadrature point quantities
            mulCAB!(qpgradT, reshape(Te, 1, :), gradN); # temperature gradient in material coordinates
            out = update!(self.material, qpflux, out, vec(qpgradT), 0.0, 0.0, loc, fes.label[i], quantity)
            if (quantity == :heatflux)   # Transform heat flux vector,  if that is "out"
                mulCAB!(out1, transpose(self.mcsys.csmat), out);# To global coord sys
                mulCAB!(out, outputcsys.csmat, out1);# To output coord sys
            end
            # Call the inspector
            idat = inspector(idat, i, fes.conn[i], ecoords, out, loc);
        end # Loop over quadrature points
    end # Loop over elements
    return idat; # return the updated inspector data
end

"""
    capacity(self::FEMMHeatDiff,  assembler::A, geom::NodalField{FFlt},  temp::NodalField{FFlt}) where {A<:AbstractSysmatAssembler}

Compute the capacity matrix.

# Arguments
- `self` = model machine,
- `assembler` = matrix assembler
- `geom` = geometry field,
- `temp` = temperature field
"""
function capacity(self::FEMMHeatDiff,  assembler::A, geom::NodalField{FFlt},  temp::NodalField{FFlt}) where {A<:AbstractSysmatAssembler}
	fes = self.integdomain.fes
	ecoords, dofnums, loc, J, elmat, elvec, elvecfix = _buffers2(self, geom, temp)
	# Precompute basis f. values + basis f. gradients wrt parametric coor
	npts, Ns, gradNparams, w, pc  =  integrationdata(self.integdomain);
	# Material
	specific_heat  =  self.material.specific_heat;
	startassembly!(assembler, size(elmat,1), size(elmat,2), count(fes),
		temp.nfreedofs, temp.nfreedofs);
	for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom, ecoords, fes.conn[i]);
		fill!(elmat, 0.0); # Initialize element matrix
		for j = 1:npts # Loop over quadrature points
			locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
			Jac = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
			ffactor = Jac*specific_heat*w[j]
			add_nnt_ut_only!(elmat, Ns[j], ffactor)
		end # Loop over quadrature points
		complete_lt!(elmat)
		gatherdofnums!(temp, dofnums, fes.conn[i]);# retrieve degrees of freedom
		assemble!(assembler, elmat, dofnums, dofnums);# assemble symmetric matrix
	end # Loop over elements
	return makematrix!(assembler);
end

function capacity(self::FEMMHeatDiff, geom::NodalField{FFlt},  temp::NodalField{FFlt})
    assembler = SysmatAssemblerSparseSymm();
    return capacity(self, assembler, geom, temp);
end

end

max_stars_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_issues_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_forks_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff


text:
include("utils.jl")

include("structures.jl")
include("sto3g.jl")
include("basfun.jl")
include("BasisSet.jl")
include("overlap.jl")
include("kinetic.jl")
include("nuclear.jl")

include("TESTS.jl")

test_utils()
test_pgbf()
test_cgbf()
test_kinetic()

test_fgamma()
test_one()
test_na2()

include("two_el.jl")
test_two_terms()
test_coul1()
test_vrr()
test_hrr()

test_geo_basis()

include("rhf.jl")
include("test_molecules.jl")

#test_h2()
#test_lih()
#test_h2o()
#test_ch4()
#test_c6h6()