text:
""" Lay out and draw wiring diagrams using Graphviz.
"""
module GraphvizWiringDiagrams
export to_graphviz, graphviz_layout

import JSON
using LightGraphs, MetaGraphs
using LinearAlgebra: normalize
using StaticArrays

import ...Theories: HomExpr
using ...WiringDiagrams, ...WiringDiagrams.WiringDiagramSerialization
import ..Graphviz
import ..Graphviz: to_graphviz
using ..Graphviz: parse_graphviz, run_graphviz
using ..WiringDiagramLayouts: BoxLayout, PortLayout, WirePoint,
  LayoutOrientation, LeftToRight, RightToLeft, TopToBottom, BottomToTop,
  is_horizontal, is_vertical, box_label, wire_label, port_sign, svector

# Drawing
#########

# Default Graphviz font. Reference: http://www.graphviz.org/doc/fontfaq.txt
const default_font = "Serif"

# Default graph, node, edge, and cell attributes.
const default_graph_attrs = Graphviz.Attributes(
  :fontname => default_font,
)
const default_node_attrs = Graphviz.Attributes(
  :fontname => default_font,
  :shape => "none",
  :width => "0",
  :height => "0",
  :margin => "0",
)
const default_edge_attrs = Graphviz.Attributes(
  :arrowsize => "0.5",
  :fontname => default_font,
)
const default_cell_attrs = Graphviz.Attributes(
  :border => "1",
  :cellpadding => "4",
)

struct GraphvizBox
  stmts::Vector{Graphviz.Statement} # Usually Graphviz.Node
  input_ports::Vector{Graphviz.NodeID}
  output_ports::Vector{Graphviz.NodeID}
end

""" Draw a wiring diagram using Graphviz.

The input `f` can also be a morphism expression, which is converted into a
wiring diagram.

# Arguments
- `graph_name="G"`: name of Graphviz digraph
- `orientation=TopToBottom`: orientation of layout.
  One of `LeftToRight`, `RightToLeft`, `TopToBottom`, or `BottomToTop`.
- `node_labels=true`: whether to label the nodes
- `labels=false`: whether to label the edges
- `label_attr=:label`: what kind of edge label to use (if `labels` is true).
  One of `:label`, `:xlabel`, `:headlabel`, or `:taillabel`.
- `port_size="24"`: minimum size of ports on box, in points
- `junction_size="0.05"`: size of junction nodes, in inches
- `outer_ports=true`: whether to display the outer box's input and output ports.
  If disabled, no incoming or outgoing wires will be shown either!
- `anchor_outer_ports=true`: whether to enforce ordering of the outer box's
  input and output, i.e., ordering of the incoming and outgoing wires
- `graph_attrs=default_graph_attrs`: top-level graph attributes
- `node_attrs=default_node_attrs`: top-level node attributes
- `edge_attrs=default_edge_attrs`: top-level edge attributes
- `cell_attrs=default_cell_attrs`: main cell attributes in node HTML-like label
"""
function to_graphviz(f::WiringDiagram;
    graph_name::String="G", orientation::LayoutOrientation=TopToBottom,
    node_labels::Bool=true, labels::Bool=false, label_attr::Symbol=:label,
    port_size::String="24", junction_size::String="0.05",
    outer_ports::Bool=true, anchor_outer_ports::Bool=true,
    graph_attrs::AbstractDict=Dict(), node_attrs::AbstractDict=Dict(),
    edge_attrs::AbstractDict=Dict(), cell_attrs::AbstractDict=Dict())::Graphviz.Graph
  @assert label_attr in (:label, :xlabel, :headlabel, :taillabel)

  # State variables.
  stmts = Graphviz.Statement[]
  port_map = Dict{Port,Graphviz.NodeID}()
  update_port_map! = (v::Int, kind::PortKind, node_ids) -> begin
    for (i, node_id) in enumerate(node_ids)
      port_map[Port(v,kind,i)] = node_id
    end
  end

  # Invisible nodes for incoming and outgoing wires.
  if outer_ports
    gv_box = graphviz_outer_box(f;
      anchor=anchor_outer_ports, orientation=orientation)
    append!(stmts, gv_box.stmts)
    update_port_map!(input_id(f), OutputPort, gv_box.input_ports)
    update_port_map!(output_id(f), InputPort, gv_box.output_ports)
  end

  # Visible nodes for boxes.
  cell_attrs = merge(default_cell_attrs, Graphviz.as_attributes(cell_attrs))
  for v in box_ids(f)
    gv_box = graphviz_box(box(f,v), box_id([v]),
      orientation=orientation, labels=node_labels, port_size=port_size,
      junction_size=junction_size, cell_attrs=cell_attrs)
    append!(stmts, gv_box.stmts)
    update_port_map!(v, InputPort, gv_box.input_ports)
    update_port_map!(v, OutputPort, gv_box.output_ports)
  end

  # Edges.
  for (i, wire) in enumerate(wires(f))
    source, target = wire.source, wire.target
    if !(haskey(port_map, source) && haskey(port_map, target))
      continue
    end
    # Use the port value to label the wire. We take the source port.
    # In most wiring diagrams, the source and target ports should yield the
    # same label, but that is not guaranteed. An advantage of choosing the
    # source port over the target port is that it will carry the
    # "more specific" type when implicit conversions are allowed.
    port = port_value(f, source)
    attrs = Graphviz.Attributes(
      :id => wire_id(Int[], i),
      :comment => edge_label(port),
    )
    if labels
      attrs[label_attr] = edge_label(port)
    end
    edge = Graphviz.Edge(port_map[source], port_map[target]; attrs...)
    push!(stmts, edge)
  end

  # Graph.
  Graphviz.Digraph(graph_name, stmts;
    graph_attrs=merge(default_graph_attrs, Graphviz.as_attributes(graph_attrs),
      Graphviz.Attributes(
        :rankdir => rank_dir(orientation)
      )),
    node_attrs=merge(default_node_attrs, Graphviz.as_attributes(node_attrs)),
    edge_attrs=merge(default_edge_attrs, Graphviz.as_attributes(edge_attrs)))
end

function to_graphviz(f::HomExpr; kw...)::Graphviz.Graph
  to_graphviz(to_wiring_diagram(f); kw...)
end

""" Graphviz box for a generic box.
"""
function graphviz_box(box::AbstractBox, node_id;
    orientation::LayoutOrientation=TopToBottom,
    labels::Bool=true, port_size::String="0",
    cell_attrs::AbstractDict=Dict(), kw...)
  # Main node.
  nin, nout = length(input_ports(box)), length(output_ports(box))
  text_label = labels ? node_label(box.value) : ""
  html_label = node_html_label(nin, nout, text_label;
    orientation=orientation, port_size=port_size, attrs=cell_attrs)
  # Note: The `id` attribute is included in the Graphviz output but is not used
  # internally by Graphviz. It is for use by downstream applications.
  # Reference: http://www.graphviz.org/doc/info/attrs.html#d:id
  node = Graphviz.Node(node_id,
    id = node_id,
    comment = node_label(box.value),
    label = html_label,
  )

  # Input and output ports.
  graphviz_port = (kind::PortKind, port::Int) -> begin
    Graphviz.NodeID(node_id, port_name(kind, port), port_anchor(kind, orientation))
  end
  inputs = [ graphviz_port(InputPort, i) for i in 1:nin ]
  outputs = [ graphviz_port(OutputPort, i) for i in 1:nout ]

  GraphvizBox([node], inputs, outputs)
end

""" Graphviz box for a junction.
"""
function graphviz_box(junction::Junction, node_id;
    junction_size::String="0", kw...)
  graphviz_junction(junction, node_id;
    comment = "junction",
    label = "",
    shape = "circle",
    style = "filled",
    fillcolor = "black",
    width = junction_size,
    height = junction_size,
  )
end

function graphviz_junction(junction::Junction, node_id; kw...)
  node = Graphviz.Node(node_id; id=node_id, kw...)
  nin, nout = length(input_ports(junction)), length(output_ports(junction))
  GraphvizBox([node],
    repeat([Graphviz.NodeID(node_id)], nin),
    repeat([Graphviz.NodeID(node_id)], nout))
end

""" Create "HTML-like" node label for a box.
"""
function node_html_label(nin::Int, nout::Int, text_label::String;
    orientation::LayoutOrientation=TopToBottom, port_size::String="0",
    attrs::AbstractDict=Dict())::Graphviz.Html
  if is_vertical(orientation)
    input_label = ports_horizontal_html_label(InputPort, nin, port_size)
    output_label = ports_horizontal_html_label(OutputPort, nout, port_size)
    Graphviz.Html("""

$(orientation == TopToBottom ? input_label : output_label)

$(escape_html(text_label))

$(orientation == BottomToTop ? input_label : output_label)

""")
  else
    input_label = ports_vertical_html_label(InputPort, nin, port_size)
    output_label = ports_vertical_html_label(OutputPort, nout, port_size)
    Graphviz.Html("""

      $(orientation == LeftToRight ? input_label : output_label)
      $(escape_html(text_label))
      $(orientation == RightToLeft ? input_label : output_label)

""")
  end
end

""" Create horizontal "HTML-like" label for the input or output ports of a box.
"""
function ports_horizontal_html_label(kind::PortKind, nports::Int,
    port_size::String="0")::Graphviz.Html
  cols = if nports > 0
    join("""	"""
         for i in 1:nports)
  else
    """	"""
  end
  Graphviz.Html("""

$cols

""")
end

""" Create vertical "HTML-like" label for the input or output ports of a box.
"""
function ports_vertical_html_label(kind::PortKind, nports::Int,
    port_size::String="0")::Graphviz.Html
  rows = if nports > 0
    join("""

"""
         for i in 1:nports)
  else
    """

"""
  end
  Graphviz.Html("""

$rows
""")
end

""" Graphviz box for the outer box of a wiring diagram.
"""
function graphviz_outer_box(f::WiringDiagram;
    anchor::Bool=true, orientation::LayoutOrientation=TopToBottom)
  # Subgraphs containing invisible nodes.
  stmts = Graphviz.Statement[]
  ninputs, noutputs = length(input_ports(f)), length(output_ports(f))
  if ninputs > 0
    push!(stmts, graphviz_outer_ports(input_id(f), InputPort, ninputs;
      anchor=anchor, orientation=orientation))
  end
  if noutputs > 0
    push!(stmts, graphviz_outer_ports(output_id(f), OutputPort, noutputs;
      anchor=anchor, orientation=orientation))
  end

  # Input and output ports.
  graphviz_port = (port::Port) -> Graphviz.NodeID(
    port_node_name(port.box, port.port),
    port_anchor(port.kind, orientation)
  )
  inputs = [ graphviz_port(Port(input_id(f), OutputPort, i)) for i in 1:ninputs ]
  outputs = [ graphviz_port(Port(output_id(f), InputPort, i)) for i in 1:noutputs ]

  GraphvizBox(stmts, inputs, outputs)
end

""" Create invisible nodes for the input or output ports of an outer box.
"""
function graphviz_outer_ports(v::Int, kind::PortKind, nports::Int;
    anchor::Bool=true, orientation::LayoutOrientation=TopToBottom)::Graphviz.Subgraph
  @assert nports > 0
  port_width = "$(round(24/72,digits=3))" # port width in inches
  nodes = [ port_node_name(v, i) for i in 1:nports ]
  stmts = Graphviz.Statement[
    Graphviz.Node(nodes[i], id=port_name(kind, i)) for i in 1:nports
  ]
  if anchor
    push!(stmts, Graphviz.Edge(nodes))
  end
  Graphviz.Subgraph(
    stmts,
    graph_attrs=Graphviz.Attributes(
      :rank => kind == InputPort ? "source" : "sink",
      :rankdir => is_vertical(orientation) ? "LR" : "TB"
    ),
    node_attrs=Graphviz.Attributes(
      :style => "invis",
      :shape => "none",
      :label => "",
      :width => is_vertical(orientation) ? port_width : "0",
      :height => is_horizontal(orientation) ? port_width : "0",
    ),
    edge_attrs=Graphviz.Attributes(
      :style => "invis",
    ),
  )
end
port_node_name(v::Int, port::Int) = string(box_id([v]), "p", port)

""" Graphviz rank direction (`rankdir`) for layout orientation.

Reference: https://www.graphviz.org/doc/info/attrs.html#k:rankdir
"""
rank_dir(orientation::LayoutOrientation) = rank_dirs[orientation]

const rank_dirs = Dict{LayoutOrientation,String}(
  TopToBottom => "TB",
  BottomToTop => "BT",
  LeftToRight => "LR",
  RightToLeft => "RL",
)

""" Graphviz anchor for port.
"""
function port_anchor(kind::PortKind, orientation::LayoutOrientation)
  port_anchors[(kind, orientation)]
end

const port_anchors = Dict{Tuple{PortKind,LayoutOrientation},String}(
  (InputPort, TopToBottom) => "n",
  (OutputPort, TopToBottom) => "s",
  (InputPort, BottomToTop) => "s",
  (OutputPort, BottomToTop) => "n",
  (InputPort, LeftToRight) => "w",
  (OutputPort, LeftToRight) => "e",
  (InputPort, RightToLeft) => "e",
  (OutputPort, RightToLeft) => "w",
)

""" Create a label for the main content of a box.
"""
node_label(box_value) = box_label(MIME("text/plain"), box_value)

""" Create a label for an edge.
"""
edge_label(port_value) = wire_label(MIME("text/plain"), port_value)

""" Encode attributes for Graphviz HTML-like labels.
"""
function html_attributes(attrs::AbstractDict)::String
  join([ "$(uppercase(string(k)))=\"$v\""
         for (k,v) in Graphviz.as_attributes(attrs) ], " ")
end

""" Escape special HTML characters: &, <, >, ", '

Adapted from [HttpCommon.jl](https://github.com/JuliaWeb/HttpCommon.jl).
"""
function escape_html(s::AbstractString)
  reduce(replace, [
    "&" => "&",
    "\"" => """,
    "'" => "'",
    "<" => "<",
    ">" => ">",
  ], init=s)
end

# Layout
########

""" Lay out wiring diagram using Graphviz.

Note: At this time, only the positions and sizes of the boxes, and the positions
of the outer ports, are used. The positions of the box ports and the splines for
the wires are ignored.
"""
function graphviz_layout(diagram::WiringDiagram; kw...)
  graph = to_graphviz(diagram; kw...)
  doc = JSON.parse(run_graphviz(graph, format="json0"))
  graphviz_layout(diagram, parse_graphviz(doc, multigraph=true))
end

function graphviz_layout(diagram::WiringDiagram, graph::MetaDiGraph)
  # Graphviz uses the standard Cartesian coordinate system (with origin in
  # bottom left corner), while our layout system uses a coordinate system with
  # centered origin and positive y-axis pointing downwards.
  orientation = inverse_rank_dirs[get_prop(graph, :rankdir)]
  bounds, pad = get_prop(graph, :bounds), get_prop(graph, :pad)
  diagram_size = bounds + 2*pad
  transform_point(p) = SVector(1,-1) .* (p - diagram_size/2)

  # Assumes vertices are in same order as created by `to_graphviz`:
  # 1. Input ports of outer box
  nin, nout = length(input_ports(diagram)), length(output_ports(diagram))
  main_dir = svector(orientation, 1.0, 0.0)
  inputs = map(enumerate(input_ports(diagram))) do (i, value)
    attrs = props(graph, i)
    pos = transform_point(attrs[:position])
    PortLayout(; value=value, position=pos, normal=-main_dir)
  end

  # 2. Output ports of outer box
  outputs = map(enumerate(output_ports(diagram))) do (i, value)
    attrs = props(graph, nin + i)
    pos = transform_point(attrs[:position])
    PortLayout(; value=value, position=pos, normal=main_dir)
  end

  # 3. Inner boxes
  # TODO: Use port positions from Graphviz layout. Obtain these from either
  # the HTML label or, more likely, an incident edge.
  layout = WiringDiagram(BoxLayout(size=diagram_size), inputs, outputs)
  node_map = Dict{Int,Int}()
  for (i, box) in enumerate(boxes(diagram))
    node = nin + nout + i
    attrs = props(graph, node)
    shape = Symbol(get(attrs, :shape, :ellipse))
    if shape == :none; shape = :rectangle end # XXX: Assume HTML label.
    @assert shape == :rectangle "Only the rectangle shape is implemented so far"
    pos, size = transform_point(attrs[:position]), attrs[:size]
    node_map[node] = add_box!(layout, Box(
      BoxLayout(; value=box.value, shape=shape, position=pos, size=size),
      layout_linear_ports(InputPort, input_ports(box), size, orientation),
      layout_linear_ports(OutputPort, output_ports(box), size, orientation)
    ))
  end

  # Add wires using spline points from Graphviz edge layout.
  function map_port(node::Int, portname::Union{String,Nothing})
    if node <= nin
      Port(input_id(layout), OutputPort, node)
    elseif node <= nin + nout
      Port(output_id(layout), InputPort, node - nin)
    else
      # XXX: Parsing the port names is hacky. The alternative is parsing the
      # ordered ports in the HTML node label and matching the port names.
      kind = startswith(portname, "out") ? OutputPort : InputPort
      port = parse(Int, portname[findfirst(r"[0-9]+", portname)])
      Port(node_map[node], kind, port)
    end
  end
  for edge in edges(graph)
    for attrs in get_prop(graph, edge, :edges)
      # Chop off start and end points to obtain intermediate points.
      points = attrs[:spline][3:end-2]
      wire_layout = map(Iterators.partition(points, 3)) do (c1, p, c2)
        WirePoint(transform_point(p),
                  normalize(transform_point(c2) - transform_point(c1)))
      end
      add_wire!(layout, Wire(wire_layout,
        map_port(src(edge), get(attrs, :tailport, nothing)),
        map_port(dst(edge), get(attrs, :headport, nothing))
      ))
    end
  end
  layout
end

""" Lay out ports linearly, equispaced along a rectangular box.

FIXME: Should this be in `WiringDiagramLayouts` as an alternative layout method?
"""
function layout_linear_ports(port_kind::PortKind, port_values::Vector,
                             box_size::StaticVector{2}, orientation::LayoutOrientation)
  n = length(port_values)
  sign = port_sign(port_kind, orientation)
  normal = svector(orientation, sign, 0.0)
  map(port_values, range(-1,1,length=n+2)[2:end-1]) do value, coeff
    pos = svector(orientation, sign, coeff) .* box_size/2
    PortLayout(; value=value, position=pos, normal=normal)
  end
end

const inverse_rank_dirs = Dict{String,LayoutOrientation}(
  "TB" => TopToBottom,
  "BT" => BottomToTop,
  "LR" => LeftToRight,
  "RL" => RightToLeft,
)

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# commented-out definitions are implemented in C

#abstract type Any <: Any end
#abstract type Type{T} end

#abstract type Vararg{T} end

#mutable struct Symbol
#    #opaque
#end

#mutable struct TypeName
#    name::Symbol
#end

#mutable struct DataType <: Type
#    name::TypeName
#    super::Type
#    parameters::Tuple
#    names::Tuple
#    types::Tuple
#    ctor
#    instance
#    size::Int32
#    abstract::Bool
#    mutable::Bool
#    pointerfree::Bool
#end

#struct Union <: Type
#    a
#    b
#end

#mutable struct TypeVar
#    name::Symbol
#    lb::Type
#    ub::Type
#end

#struct UnionAll
#    var::TypeVar
#    body
#end

#struct Nothing
#end
#const nothing = Nothing()

#abstract type AbstractArray{T,N} end
#abstract type DenseArray{T,N} <: AbstractArray{T,N} end

#mutable struct Array{T,N} <: DenseArray{T,N}
#end

#mutable struct Module
#    name::Symbol
#end

#mutable struct Method
#end

#mutable struct MethodInstance
#end

#mutable struct CodeInstance
#end

#mutable struct CodeInfo
#end

#mutable struct TypeMapLevel
#end

#mutable struct TypeMapEntry
#end

#abstract type Ref{T} end
#primitive type Ptr{T} <: Ref{T} {32|64} end

# types for the front end

#mutable struct Expr
#    head::Symbol
#    args::Array{Any,1}
#end

#struct LineNumberNode
#    line::Int
#    file::Union{Symbol,Nothing}
#end

#struct LineInfoNode
#    method::Any
#    file::Symbol
#    line::Int
#    inlined_at::Int
#end

#struct GotoNode
#    label::Int
#end

#struct GotoIfNot
#    cond::Any
#    dest::Int
#end

#struct ReturnNode
#    val::Any
#end

#struct PiNode
#    val
#    typ
#end

#struct PhiNode
#    edges::Vector{Any}
#    values::Vector{Any}
#end

#struct PhiCNode
#    values::Vector{Any}
#end

#struct UpsilonNode
#    val
#end

#struct QuoteNode
#    value
#end

#struct GlobalRef
#    mod::Module
#    name::Symbol
#end

#mutable struct Task
#    parent::Task
#    storage::Any
#    state::Symbol
#    donenotify::Any
#    result::Any
#    exception::Any
#    backtrace::Any
#    logstate::Any
#    code::Any
#end

export
    # key types
    Any, DataType, Vararg, NTuple,
    Tuple, Type, UnionAll, TypeVar, Union, Nothing, Cvoid,
    AbstractArray, DenseArray, NamedTuple,
    # special objects
    Function, Method,
    Module, Symbol, Task, Array, UndefInitializer, undef, WeakRef, VecElement,
    # numeric types
    Number, Real, Integer, Bool, Ref, Ptr,
    AbstractFloat, Float16, Float32, Float64,
    Signed, Int, Int8, Int16, Int32, Int64, Int128,
    Unsigned, UInt, UInt8, UInt16, UInt32, UInt64, UInt128,
    # string types
    AbstractChar, Char, AbstractString, String, IO,
    # errors
    ErrorException, BoundsError, DivideError, DomainError, Exception,
    InterruptException, InexactError, OutOfMemoryError, ReadOnlyMemoryError,
    OverflowError, StackOverflowError, SegmentationFault, UndefRefError, UndefVarError,
    TypeError, ArgumentError, MethodError, AssertionError, LoadError, InitError,
    UndefKeywordError,
    # AST representation
    Expr, QuoteNode, LineNumberNode, GlobalRef,
    # object model functions
    fieldtype, getfield, setfield!, nfields, throw, tuple, ===, isdefined, eval, ifelse,
    # sizeof    # not exported, to avoid conflicting with Base.sizeof
    # type reflection
    <:, typeof, isa, typeassert,
    # method reflection
    applicable, invoke,
    # constants
    nothing, Main

const getproperty = getfield
const setproperty! = setfield!

abstract type Number end
abstract type Real     <: Number end
abstract type AbstractFloat <: Real end
abstract type Integer  <: Real end
abstract type Signed   <: Integer end
abstract type Unsigned <: Integer end

primitive type Float16 <: AbstractFloat 16 end
primitive type Float32 <: AbstractFloat 32 end
primitive type Float64 <: AbstractFloat 64 end

#primitive type Bool <: Integer 8 end
abstract type AbstractChar end
primitive type Char <: AbstractChar 32 end

primitive type Int8    <: Signed   8 end
#primitive type UInt8   <: Unsigned 8 end
primitive type Int16   <: Signed   16 end
primitive type UInt16  <: Unsigned 16 end
#primitive type Int32   <: Signed   32 end
#primitive type UInt32  <: Unsigned 32 end
#primitive type Int64   <: Signed   64 end
#primitive type UInt64  <: Unsigned 64 end
primitive type Int128  <: Signed   128 end
primitive type UInt128 <: Unsigned 128 end

if Int === Int64
    const UInt = UInt64
else
    const UInt = UInt32
end

function iterate end
function Typeof end
ccall(:jl_toplevel_eval_in, Any, (Any, Any),
      Core, quote
      (f::typeof(Typeof))(x) = ($(_expr(:meta,:nospecialize,:x)); isa(x,Type) ? Type{x} : typeof(x))
      end)

macro nospecialize(x)
    _expr(:meta, :nospecialize, x)
end

Expr(@nospecialize args...) = _expr(args...)

abstract type Exception end
struct ErrorException <: Exception
    msg::AbstractString
end

macro _inline_meta()
    Expr(:meta, :inline)
end

macro _noinline_meta()
    Expr(:meta, :noinline)
end

struct BoundsError <: Exception
    a::Any
    i::Any
    BoundsError() = new()
    BoundsError(@nospecialize(a)) = (@_noinline_meta; new(a))
    BoundsError(@nospecialize(a), i) = (@_noinline_meta; new(a,i))
end
struct DivideError         <: Exception end
struct OutOfMemoryError    <: Exception end
struct ReadOnlyMemoryError <: Exception end
struct SegmentationFault   <: Exception end
struct StackOverflowError  <: Exception end
struct UndefRefError       <: Exception end
struct UndefVarError <: Exception
    var::Symbol
end
struct InterruptException <: Exception end
struct DomainError <: Exception
    val
    msg::AbstractString
    DomainError(@nospecialize(val)) = (@_noinline_meta; new(val, ""))
    DomainError(@nospecialize(val), @nospecialize(msg)) = (@_noinline_meta; new(val, msg))
end
struct TypeError <: Exception
    # `func` is the name of the builtin function that encountered a type error,
    # the name of the type that hit an error in its definition or application, or
    # some other brief description of where the error happened.
    # `context` optionally adds extra detail, e.g. the name of the type parameter
    # that got a bad value.
    func::Symbol
    context::Union{AbstractString,Symbol}
    expected::Type
    got
    TypeError(func, context, @nospecialize(expected::Type), @nospecialize(got)) =
        new(func, context, expected, got)
end
TypeError(where, @nospecialize(expected::Type), @nospecialize(got)) =
    TypeError(Symbol(where), "", expected, got)
struct InexactError <: Exception
    func::Symbol
    T  # Type
    val
    InexactError(f::Symbol, @nospecialize(T), @nospecialize(val)) = (@_noinline_meta; new(f, T, val))
end
struct OverflowError <: Exception
    msg::AbstractString
end

struct ArgumentError <: Exception
    msg::AbstractString
end
struct UndefKeywordError <: Exception
    var::Symbol
end

struct MethodError <: Exception
    f
    args
    world::UInt
    MethodError(@nospecialize(f), @nospecialize(args), world::UInt) = new(f, args, world)
end
const typemax_UInt = ccall(:jl_typemax_uint, Any, (Any,), UInt)
MethodError(@nospecialize(f), @nospecialize(args)) = MethodError(f, args, typemax_UInt)

struct AssertionError <: Exception
    msg::AbstractString
end
AssertionError() = AssertionError("")

abstract type WrappedException <: Exception end

struct LoadError <: WrappedException
    file::AbstractString
    line::Int
    error
end

struct InitError <: WrappedException
    mod::Symbol
    error
end

String(s::String) = s  # no constructor yet

const Cvoid = Nothing
Nothing() = nothing

# This should always be inlined
getptls() = ccall(:jl_get_ptls_states, Ptr{Cvoid}, ())

include(m::Module, fname::String) = ccall(:jl_load_, Any, (Any, Any), m, fname)

eval(m::Module, @nospecialize(e)) = ccall(:jl_toplevel_eval_in, Any, (Any, Any), m, e)

kwfunc(@nospecialize(f)) = ccall(:jl_get_keyword_sorter, Any, (Any,), f)

kwftype(@nospecialize(t)) = typeof(ccall(:jl_get_kwsorter, Any, (Any,), t))

mutable struct Box
    contents::Any
    Box(@nospecialize(x)) = new(x)
    Box() = new()
end

# constructors for built-in types

mutable struct WeakRef
    value
    WeakRef() = WeakRef(nothing)
    WeakRef(@nospecialize(v)) = ccall(:jl_gc_new_weakref_th, Ref{WeakRef},
                                      (Ptr{Cvoid}, Any), getptls(), v)
end

TypeVar(n::Symbol) = _typevar(n, Union{}, Any)
TypeVar(n::Symbol, @nospecialize(ub)) = _typevar(n, Union{}, ub)
TypeVar(n::Symbol, @nospecialize(lb), @nospecialize(ub)) = _typevar(n, lb, ub)

UnionAll(v::TypeVar, @nospecialize(t)) = ccall(:jl_type_unionall, Any, (Any, Any), v, t)

Tuple{}() = ()

struct VecElement{T}
    value::T
    VecElement{T}(value::T) where {T} = new(value) # disable converting constructor in Core
end
VecElement(arg::T) where {T} = VecElement{T}(arg)

_new(typ::Symbol, argty::Symbol) = eval(Core, :($typ(@nospecialize n::$argty) = $(Expr(:new, typ, :n))))
_new(:GotoNode, :Int)
_new(:NewvarNode, :SlotNumber)
_new(:QuoteNode, :Any)
_new(:SSAValue, :Int)
_new(:Argument, :Int)
_new(:ReturnNode, :Any)
eval(Core, :(ReturnNode() = $(Expr(:new, :ReturnNode)))) # unassigned val indicates unreachable
eval(Core, :(GotoIfNot(@nospecialize(cond), dest::Int) = $(Expr(:new, :GotoIfNot, :cond, :dest))))
eval(Core, :(LineNumberNode(l::Int) = $(Expr(:new, :LineNumberNode, :l, nothing))))
eval(Core, :(LineNumberNode(l::Int, @nospecialize(f)) = $(Expr(:new, :LineNumberNode, :l, :f))))
LineNumberNode(l::Int, f::String) = LineNumberNode(l, Symbol(f))
eval(Core, :(GlobalRef(m::Module, s::Symbol) = $(Expr(:new, :GlobalRef, :m, :s))))
eval(Core, :(SlotNumber(n::Int) = $(Expr(:new, :SlotNumber, :n))))
eval(Core, :(TypedSlot(n::Int, @nospecialize(t)) = $(Expr(:new, :TypedSlot, :n, :t))))
eval(Core, :(PhiNode(edges::Array{Any, 1}, values::Array{Any, 1}) = $(Expr(:new, :PhiNode, :edges, :values))))
eval(Core, :(PiNode(val, typ) = $(Expr(:new, :PiNode, :val, :typ))))
eval(Core, :(PhiCNode(values::Array{Any, 1}) = $(Expr(:new, :PhiCNode, :values))))
eval(Core, :(UpsilonNode(val) = $(Expr(:new, :UpsilonNode, :val))))
eval(Core, :(UpsilonNode() = $(Expr(:new, :UpsilonNode))))
eval(Core, :(LineInfoNode(@nospecialize(method), file::Symbol, line::Int, inlined_at::Int) =
             $(Expr(:new, :LineInfoNode, :method, :file, :line, :inlined_at))))
eval(Core, :(CodeInstance(mi::MethodInstance, @nospecialize(rettype), @nospecialize(inferred_const),
                          @nospecialize(inferred), const_flags::Int32,
                          min_world::UInt, max_world::UInt) =
                ccall(:jl_new_codeinst, Ref{CodeInstance}, (Any, Any, Any, Any, Int32, UInt, UInt),
                    mi, rettype, inferred_const, inferred, const_flags, min_world, max_world)))
eval(Core, :(Const(@nospecialize(v)) = $(Expr(:new, :Const, :v, false))))
eval(Core, :(Const(@nospecialize(v), actual::Bool) = $(Expr(:new, :Const, :v, :actual))))
eval(Core, :(PartialStruct(@nospecialize(typ), fields::Array{Any, 1}) = $(Expr(:new, :PartialStruct, :typ, :fields))))

Module(name::Symbol=:anonymous, std_imports::Bool=true) = ccall(:jl_f_new_module, Ref{Module}, (Any, Bool), name, std_imports)

function _Task(@nospecialize(f), reserved_stack::Int, completion_future)
    return ccall(:jl_new_task, Ref{Task}, (Any, Any, Int), f, completion_future, reserved_stack)
end

# simple convert for use by constructors of types in Core
# note that there is no actual conversion defined here,
# so the methods and ccall's in Core aren't permitted to use convert
convert(::Type{Any}, @nospecialize(x)) = x
convert(::Type{T}, x::T) where {T} = x
cconvert(::Type{T}, x) where {T} = convert(T, x)
unsafe_convert(::Type{T}, x::T) where {T} = x

const NTuple{N,T} = Tuple{Vararg{T,N}}

## primitive Array constructors
struct UndefInitializer end
const undef = UndefInitializer()
# type and dimensionality specified, accepting dims as series of Ints
Array{T,1}(::UndefInitializer, m::Int) where {T} =
    ccall(:jl_alloc_array_1d, Array{T,1}, (Any, Int), Array{T,1}, m)
Array{T,2}(::UndefInitializer, m::Int, n::Int) where {T} =
    ccall(:jl_alloc_array_2d, Array{T,2}, (Any, Int, Int), Array{T,2}, m, n)
Array{T,3}(::UndefInitializer, m::Int, n::Int, o::Int) where {T} =
    ccall(:jl_alloc_array_3d, Array{T,3}, (Any, Int, Int, Int), Array{T,3}, m, n, o)
Array{T,N}(::UndefInitializer, d::Vararg{Int,N}) where {T,N} =
    ccall(:jl_new_array, Array{T,N}, (Any, Any), Array{T,N}, d)
# type and dimensionality specified, accepting dims as tuples of Ints
Array{T,1}(::UndefInitializer, d::NTuple{1,Int}) where {T} = Array{T,1}(undef, getfield(d,1))
Array{T,2}(::UndefInitializer, d::NTuple{2,Int}) where {T} = Array{T,2}(undef, getfield(d,1), getfield(d,2))
Array{T,3}(::UndefInitializer, d::NTuple{3,Int}) where {T} = Array{T,3}(undef, getfield(d,1), getfield(d,2), getfield(d,3))
Array{T,N}(::UndefInitializer, d::NTuple{N,Int}) where {T,N} = ccall(:jl_new_array, Array{T,N}, (Any, Any), Array{T,N}, d)
# type but not dimensionality specified
Array{T}(::UndefInitializer, m::Int) where {T} = Array{T,1}(undef, m)
Array{T}(::UndefInitializer, m::Int, n::Int) where {T} = Array{T,2}(undef, m, n)
Array{T}(::UndefInitializer, m::Int, n::Int, o::Int) where {T} = Array{T,3}(undef, m, n, o)
Array{T}(::UndefInitializer, d::NTuple{N,Int}) where {T,N} = Array{T,N}(undef, d)
# empty vector constructor
Array{T,1}() where {T} = Array{T,1}(undef, 0)

(Array{T,N} where T)(x::AbstractArray{S,N}) where {S,N} = Array{S,N}(x)

Array(A::AbstractArray{T,N})    where {T,N}   = Array{T,N}(A)
Array{T}(A::AbstractArray{S,N}) where {T,N,S} = Array{T,N}(A)

AbstractArray{T}(A::AbstractArray{S,N}) where {T,S,N} = AbstractArray{T,N}(A)

# primitive Symbol constructors
eval(Core, :(function Symbol(s::String)
    $(Expr(:meta, :pure))
    return ccall(:jl_symbol_n, Ref{Symbol}, (Ptr{UInt8}, Int),
                 ccall(:jl_string_ptr, Ptr{UInt8}, (Any,), s),
                 sizeof(s))
end))
function Symbol(a::Array{UInt8,1})
    return ccall(:jl_symbol_n, Ref{Symbol}, (Ptr{UInt8}, Int),
                 ccall(:jl_array_ptr, Ptr{UInt8}, (Any,), a),
                 Intrinsics.arraylen(a))
end
Symbol(s::Symbol) = s

# module providing the IR object model
module IR
export CodeInfo, MethodInstance, CodeInstance, GotoNode, GotoIfNot, ReturnNode,
    NewvarNode, SSAValue, Slot, SlotNumber, TypedSlot, Argument,
    PiNode, PhiNode, PhiCNode, UpsilonNode, LineInfoNode,
    Const, PartialStruct

import Core: CodeInfo, MethodInstance, CodeInstance, GotoNode, GotoIfNot, ReturnNode,
    NewvarNode, SSAValue, Slot, SlotNumber, TypedSlot, Argument,
    PiNode, PhiNode, PhiCNode, UpsilonNode, LineInfoNode,
    Const, PartialStruct

end

# docsystem basics
const unescape = Symbol("hygienic-scope")
macro doc(x...)
    docex = atdoc(__source__, __module__, x...)
    isa(docex, Expr) && docex.head === :escape && return docex
    return Expr(:escape, Expr(unescape, docex, typeof(atdoc).name.module))
end
macro __doc__(x)
    return Expr(:escape, Expr(:block, Expr(:meta, :doc), x))
end
atdoc     = (source, mod, str, expr) -> Expr(:escape, expr)
atdoc!(λ) = global atdoc = λ

# macros for big integer syntax
macro int128_str end
macro uint128_str end
macro big_str end

# macro for command syntax
macro cmd end

# simple stand-alone print definitions for debugging
abstract type IO end
struct CoreSTDOUT <: IO end
struct CoreSTDERR <: IO end
const stdout = CoreSTDOUT()
const stderr = CoreSTDERR()
io_pointer(::CoreSTDOUT) = Intrinsics.pointerref(Intrinsics.cglobal(:jl_uv_stdout, Ptr{Cvoid}), 1, 1)
io_pointer(::CoreSTDERR) = Intrinsics.pointerref(Intrinsics.cglobal(:jl_uv_stderr, Ptr{Cvoid}), 1, 1)

unsafe_write(io::IO, x::Ptr{UInt8}, nb::UInt) =
    (ccall(:jl_uv_puts, Cvoid, (Ptr{Cvoid}, Ptr{UInt8}, UInt), io_pointer(io), x, nb); nb)
unsafe_write(io::IO, x::Ptr{UInt8}, nb::Int) =
    (ccall(:jl_uv_puts, Cvoid, (Ptr{Cvoid}, Ptr{UInt8}, Int), io_pointer(io), x, nb); nb)
write(io::IO, x::UInt8) =
    (ccall(:jl_uv_putb, Cvoid, (Ptr{Cvoid}, UInt8), io_pointer(io), x); 1)
function write(io::IO, x::String)
    nb = sizeof(x)
    unsafe_write(io, ccall(:jl_string_ptr, Ptr{UInt8}, (Any,), x), nb)
    return nb
end

show(io::IO, @nospecialize x) = ccall(:jl_static_show, Cvoid, (Ptr{Cvoid}, Any), io_pointer(io), x)
print(io::IO, x::AbstractChar) = ccall(:jl_uv_putc, Cvoid, (Ptr{Cvoid}, Char), io_pointer(io), x)
print(io::IO, x::String) = (write(io, x); nothing)
print(io::IO, @nospecialize x) = show(io, x)
print(io::IO, @nospecialize(x), @nospecialize a...) = (print(io, x); print(io, a...))
println(io::IO) = (write(io, 0x0a); nothing) # 0x0a = '\n'
println(io::IO, @nospecialize x...) = (print(io, x...); println(io))

show(@nospecialize a) = show(stdout, a)
print(@nospecialize a...) = print(stdout, a...)
println(@nospecialize a...) = println(stdout, a...)

struct GeneratedFunctionStub
    gen
    argnames::Array{Any,1}
    spnames::Union{Nothing, Array{Any,1}}
    line::Int
    file::Symbol
    expand_early::Bool
end

# invoke and wrap the results of @generated
function (g::GeneratedFunctionStub)(@nospecialize args...)
    body = g.gen(args...)
    if body isa CodeInfo
        return body
    end
    lam = Expr(:lambda, g.argnames,
               Expr(Symbol("scope-block"),
                    Expr(:block,
                         LineNumberNode(g.line, g.file),
                         Expr(:meta, :push_loc, g.file, Symbol("@generated body")),
                         Expr(:return, body),
                         Expr(:meta, :pop_loc))))
    if g.spnames === nothing
        return lam
    else
        return Expr(Symbol("with-static-parameters"), lam, g.spnames...)
    end
end

NamedTuple() = NamedTuple{(),Tuple{}}(())

NamedTuple{names}(args::Tuple) where {names} = NamedTuple{names,typeof(args)}(args)

using .Intrinsics: sle_int, add_int

eval(Core, :(NamedTuple{names,T}(args::T) where {names, T <: Tuple} =
             $(Expr(:splatnew, :(NamedTuple{names,T}), :args))))

# constructors for built-in types

import .Intrinsics: eq_int, trunc_int, lshr_int, sub_int, shl_int, bitcast, sext_int, zext_int, and_int

throw_inexacterror(f::Symbol, ::Type{T}, val) where {T} = (@_noinline_meta; throw(InexactError(f, T, val)))

function is_top_bit_set(x)
    @_inline_meta
    eq_int(trunc_int(UInt8, lshr_int(x, sub_int(shl_int(sizeof(x), 3), 1))), trunc_int(UInt8, 1))
end

function is_top_bit_set(x::Union{Int8,UInt8})
    @_inline_meta
    eq_int(lshr_int(x, 7), trunc_int(typeof(x), 1))
end

function check_top_bit(::Type{To}, x) where {To}
    @_inline_meta
    is_top_bit_set(x) && throw_inexacterror(:check_top_bit, To, x)
    x
end

function checked_trunc_sint(::Type{To}, x::From) where {To,From}
    @_inline_meta
    y = trunc_int(To, x)
    back = sext_int(From, y)
    eq_int(x, back) || throw_inexacterror(:trunc, To, x)
    y
end

function checked_trunc_uint(::Type{To}, x::From) where {To,From}
    @_inline_meta
    y = trunc_int(To, x)
    back = zext_int(From, y)
    eq_int(x, back) || throw_inexacterror(:trunc, To, x)
    y
end

toInt8(x::Int8)       = x
toInt8(x::Int16)      = checked_trunc_sint(Int8, x)
toInt8(x::Int32)      = checked_trunc_sint(Int8, x)
toInt8(x::Int64)      = checked_trunc_sint(Int8, x)
toInt8(x::Int128)     = checked_trunc_sint(Int8, x)
toInt8(x::UInt8)      = bitcast(Int8, check_top_bit(Int8, x))
toInt8(x::UInt16)     = checked_trunc_sint(Int8, check_top_bit(Int8, x))
toInt8(x::UInt32)     = checked_trunc_sint(Int8, check_top_bit(Int8, x))
toInt8(x::UInt64)     = checked_trunc_sint(Int8, check_top_bit(Int8, x))
toInt8(x::UInt128)    = checked_trunc_sint(Int8, check_top_bit(Int8, x))
toInt8(x::Bool)       = and_int(bitcast(Int8, x), Int8(1))
toInt16(x::Int8)      = sext_int(Int16, x)
toInt16(x::Int16)     = x
toInt16(x::Int32)     = checked_trunc_sint(Int16, x)
toInt16(x::Int64)     = checked_trunc_sint(Int16, x)
toInt16(x::Int128)    = checked_trunc_sint(Int16, x)
toInt16(x::UInt8)     = zext_int(Int16, x)
toInt16(x::UInt16)    = bitcast(Int16, check_top_bit(Int16, x))
toInt16(x::UInt32)    = checked_trunc_sint(Int16, check_top_bit(Int16, x))
toInt16(x::UInt64)    = checked_trunc_sint(Int16, check_top_bit(Int16, x))
toInt16(x::UInt128)   = checked_trunc_sint(Int16, check_top_bit(Int16, x))
toInt16(x::Bool)      = and_int(zext_int(Int16, x), Int16(1))
toInt32(x::Int8)      = sext_int(Int32, x)
toInt32(x::Int16)     = sext_int(Int32, x)
toInt32(x::Int32)     = x
toInt32(x::Int64)     = checked_trunc_sint(Int32, x)
toInt32(x::Int128)    = checked_trunc_sint(Int32, x)
toInt32(x::UInt8)     = zext_int(Int32, x)
toInt32(x::UInt16)    = zext_int(Int32, x)
toInt32(x::UInt32)    = bitcast(Int32, check_top_bit(Int32, x))
toInt32(x::UInt64)    = checked_trunc_sint(Int32, check_top_bit(Int32, x))
toInt32(x::UInt128)   = checked_trunc_sint(Int32, check_top_bit(Int32, x))
toInt32(x::Bool)      = and_int(zext_int(Int32, x), Int32(1))
toInt64(x::Int8)      = sext_int(Int64, x)
toInt64(x::Int16)     = sext_int(Int64, x)
toInt64(x::Int32)     = sext_int(Int64, x)
toInt64(x::Int64)     = x
toInt64(x::Int128)    = checked_trunc_sint(Int64, x)
toInt64(x::UInt8)     = zext_int(Int64, x)
toInt64(x::UInt16)    = zext_int(Int64, x)
toInt64(x::UInt32)    = zext_int(Int64, x)
toInt64(x::UInt64)    = bitcast(Int64, check_top_bit(Int64, x))
toInt64(x::UInt128)   = checked_trunc_sint(Int64, check_top_bit(Int64, x))
toInt64(x::Bool)      = and_int(zext_int(Int64, x), Int64(1))
toInt128(x::Int8)     = sext_int(Int128, x)
toInt128(x::Int16)    = sext_int(Int128, x)
toInt128(x::Int32)    = sext_int(Int128, x)
toInt128(x::Int64)    = sext_int(Int128, x)
toInt128(x::Int128)   = x
toInt128(x::UInt8)    = zext_int(Int128, x)
toInt128(x::UInt16)   = zext_int(Int128, x)
toInt128(x::UInt32)   = zext_int(Int128, x)
toInt128(x::UInt64)   = zext_int(Int128, x)
toInt128(x::UInt128)  = bitcast(Int128, check_top_bit(Int128, x))
toInt128(x::Bool)     = and_int(zext_int(Int128, x), Int128(1))
toUInt8(x::Int8)      = bitcast(UInt8, check_top_bit(UInt8, x))
toUInt8(x::Int16)     = checked_trunc_uint(UInt8, x)
toUInt8(x::Int32)     = checked_trunc_uint(UInt8, x)
toUInt8(x::Int64)     = checked_trunc_uint(UInt8, x)
toUInt8(x::Int128)    = checked_trunc_uint(UInt8, x)
toUInt8(x::UInt8)     = x
toUInt8(x::UInt16)    = checked_trunc_uint(UInt8, x)
toUInt8(x::UInt32)    = checked_trunc_uint(UInt8, x)
toUInt8(x::UInt64)    = checked_trunc_uint(UInt8, x)
toUInt8(x::UInt128)   = checked_trunc_uint(UInt8, x)
toUInt8(x::Bool)      = and_int(bitcast(UInt8, x), UInt8(1))
toUInt16(x::Int8)     = sext_int(UInt16, check_top_bit(UInt16, x))
toUInt16(x::Int16)    = bitcast(UInt16, check_top_bit(UInt16, x))
toUInt16(x::Int32)    = checked_trunc_uint(UInt16, x)
toUInt16(x::Int64)    = checked_trunc_uint(UInt16, x)
toUInt16(x::Int128)   = checked_trunc_uint(UInt16, x)
toUInt16(x::UInt8)    = zext_int(UInt16, x)
toUInt16(x::UInt16)   = x
toUInt16(x::UInt32)   = checked_trunc_uint(UInt16, x)
toUInt16(x::UInt64)   = checked_trunc_uint(UInt16, x)
toUInt16(x::UInt128)  = checked_trunc_uint(UInt16, x)
toUInt16(x::Bool)     = and_int(zext_int(UInt16, x), UInt16(1))
toUInt32(x::Int8)     = sext_int(UInt32, check_top_bit(UInt32, x))
toUInt32(x::Int16)    = sext_int(UInt32, check_top_bit(UInt32, x))
toUInt32(x::Int32)    = bitcast(UInt32, check_top_bit(UInt32, x))
toUInt32(x::Int64)    = checked_trunc_uint(UInt32, x)
toUInt32(x::Int128)   = checked_trunc_uint(UInt32, x)
toUInt32(x::UInt8)    = zext_int(UInt32, x)
toUInt32(x::UInt16)   = zext_int(UInt32, x)
toUInt32(x::UInt32)   = x
toUInt32(x::UInt64)   = checked_trunc_uint(UInt32, x)
toUInt32(x::UInt128)  = checked_trunc_uint(UInt32, x)
toUInt32(x::Bool)     = and_int(zext_int(UInt32, x), UInt32(1))
toUInt64(x::Int8)     = sext_int(UInt64, check_top_bit(UInt64, x))
toUInt64(x::Int16)    = sext_int(UInt64, check_top_bit(UInt64, x))
toUInt64(x::Int32)    = sext_int(UInt64, check_top_bit(UInt64, x))
toUInt64(x::Int64)    = bitcast(UInt64, check_top_bit(UInt64, x))
toUInt64(x::Int128)   = checked_trunc_uint(UInt64, x)
toUInt64(x::UInt8)    = zext_int(UInt64, x)
toUInt64(x::UInt16)   = zext_int(UInt64, x)
toUInt64(x::UInt32)   = zext_int(UInt64, x)
toUInt64(x::UInt64)   = x
toUInt64(x::UInt128)  = checked_trunc_uint(UInt64, x)
toUInt64(x::Bool)     = and_int(zext_int(UInt64, x), UInt64(1))
toUInt128(x::Int8)    = sext_int(UInt128, check_top_bit(UInt128, x))
toUInt128(x::Int16)   = sext_int(UInt128, check_top_bit(UInt128, x))
toUInt128(x::Int32)   = sext_int(UInt128, check_top_bit(UInt128, x))
toUInt128(x::Int64)   = sext_int(UInt128, check_top_bit(UInt128, x))
toUInt128(x::Int128)  = bitcast(UInt128, check_top_bit(UInt128, x))
toUInt128(x::UInt8)   = zext_int(UInt128, x)
toUInt128(x::UInt16)  = zext_int(UInt128, x)
toUInt128(x::UInt32)  = zext_int(UInt128, x)
toUInt128(x::UInt64)  = zext_int(UInt128, x)
toUInt128(x::UInt128) = x
toUInt128(x::Bool)    = and_int(zext_int(UInt128, x), UInt128(1))

# TODO: this is here to work around the 4 method limit in inference (#23210).
const BuiltinInts = Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8, Bool}
Int8(x::BuiltinInts)    = toInt8(x)::Int8
Int16(x::BuiltinInts)   = toInt16(x)::Int16
Int32(x::BuiltinInts)   = toInt32(x)::Int32
Int64(x::BuiltinInts)   = toInt64(x)::Int64
Int128(x::BuiltinInts)  = toInt128(x)::Int128
UInt8(x::BuiltinInts)   = toUInt8(x)::UInt8
UInt16(x::BuiltinInts)  = toUInt16(x)::UInt16
UInt32(x::BuiltinInts)  = toUInt32(x)::UInt32
UInt64(x::BuiltinInts)  = toUInt64(x)::UInt64
UInt128(x::BuiltinInts) = toUInt128(x)::UInt128

(::Type{T})(x::T) where {T<:Number} = x

Int(x::Ptr)  = bitcast(Int, x)
UInt(x::Ptr) = bitcast(UInt, x)
if Int === Int32
Int64(x::Ptr) = Int64(UInt32(x))
UInt64(x::Ptr) = UInt64(UInt32(x))
end
Ptr{T}(x::Union{Int,UInt,Ptr}) where {T} = bitcast(Ptr{T}, x)
Ptr{T}() where {T} = Ptr{T}(0)

Signed(x::UInt8)    = Int8(x)
Unsigned(x::Int8)   = UInt8(x)
Signed(x::UInt16)   = Int16(x)
Unsigned(x::Int16)  = UInt16(x)
Signed(x::UInt32)   = Int32(x)
Unsigned(x::Int32)  = UInt32(x)
Signed(x::UInt64)   = Int64(x)
Unsigned(x::Int64)  = UInt64(x)
Signed(x::UInt128)  = Int128(x)
Unsigned(x::Int128) = UInt128(x)

Signed(x::Union{Float32, Float64, Bool})   = Int(x)
Unsigned(x::Union{Float32, Float64, Bool}) = UInt(x)

Integer(x::Integer) = x
Integer(x::Union{Float32, Float64}) = Int(x)

# Binding for the julia parser, called as
#
#    Core._parse(text, filename, offset, options)
#
# Parse Julia code from the buffer `text`, starting at `offset` and attributing
# it to `filename`. `text` may be a `String` or `svec(ptr::Ptr{UInt8},
# len::Int)` for a raw unmanaged buffer. `options` should be one of `:atom`,
# `:statement` or `:all`, indicating how much the parser will consume.
#
# `_parse` must return an `svec` containing an `Expr` and the new offset as an
# `Int`.
#
# The internal jl_parse which will call into Core._parse if not `nothing`.
_parse = nothing

ccall(:jl_set_istopmod, Cvoid, (Any, Bool), Core, true)

ext: jl

max_stars_repo_name: FedericoStra/julia

max_issues_repo_name: FedericoStra/julia

max_forks_repo_name: FedericoStra/julia


text:
module NetworkDesignModels

    import JuMP, Gurobi, LightGraphs, MathProgBase
    using LinearAlgebra, Parameters
    using IterTools, Combinatorics, ProgressMeter

    include("transitnetwork.jl")
    include("masterproblem.jl")
    include("utils.jl")
    include("subproblem.jl")
    include("subproblem-traveltimes.jl")
    include("subproblem-model.jl")
    include("subproblem-utils.jl")
end


text:
using Flexpart
using Test
using Dates

# Flexpart.NCF_OUTPUT
# @show Flexpart.option2dict("outgrid")
# Flexpart.find_ncf()
# # Flexpart.outclear()
# # Flexpart.outinfo()
# Flexpart.outgrid()
# @show Flexpart.conc(9)
# Flexpart.fields(9)
# outgrid = Flexpart.Outgrid(5.009, 50.353, 1111, 593, 0.001, 0.001, [100.0])

# @testset "Flexpart.jl" begin

    ###################################
    ############ UTILS ################
    ###################################$
    area = [60.1875, -5.0625, 34.875, 20.25]
    area = [60, -5, 35, 20]
    grid1 = 0.22
    grid1 = 0.5
    alons = -180.0:grid1:180.0 |> collect

    outb = Flexpart.outer_vals(alons, (area[2], area[4]))
    inb = Flexpart.inner_vals(alons, (area[2], area[4]))

    outa = Flexpart.outer_area(area, grid1)
    ina = Flexpart.inner_area(area, grid1)

    ###################################
    ###### TEST FLEXPART DIR ##########
    ###################################
    dirpath = "test/fp_dir_test"
    newdir = Flexpart.create(dirpath)
    fpdir = FlexpartDir(dirpath)
    # pathnam = Flexpart.pathnames(newdir)
    fpdir[:output]
    fpdir[:input] = "/home/tcarion/.julia/dev/Flexpart/test/fe_template/fedir/input"
    fpdir[:input]
    write(fpdir)

    ###################################
    ###### TEST FLEXPART INPUT ########
    ###################################
    inputdir = "test/fp_dir_test/input"
    # av = Flexpart.Available(fpdir)
    # new_av = Flexpart.update(av, "/home/tcarion/.julia/dev/Flexpart/test/fe_template/fedir/input")
    # new_av2 = Flexpart.update(av, DateTime(2021, 8, 8):Dates.Hour(1):DateTime(2021, 8, 9)|>collect, "PREF")
    newav = Flexpart.updated_available(fpdir)
    Flexpart.write(fpdir, newav)

    ###################################
    ###### TEST FLEXPART OPTIONS ######
    ###################################
    fpoptions = FlexpartOptions(dirpath)
    fpoptions["COMMAND"][:command][:ldirect] = 9
    area = [50, 4, 52, 6]
    newv = area2outgrid(area)
    set!(fpoptions["OUTGRID"][:outgrid], newv)
    write(fpoptions, pwd())

    fopt = Flexpart.getnamelists(Flexpart.getdir(fpdir, :options))

    comp = Flexpart.compare(
        "/home/tcarion/.julia/dev/Flexpart/test/fp_dir_test/output/COMMAND.namelist",
        "/home/tcarion/.julia/dev/Flexpart/test/fp_dir_test/output/COMMAND.namelist2")

    # dirs = files[isdir.(files)]
    ###################################
    ###### TEST FLEXPART OUTPUTS ######
    ###################################
    output_files = ncf_files(fpdir)
    fpoutput = FlexpartOutput(output_files[1])

    var2d = Flexpart.variables2d(fpoutput)
    vars = Flexpart.variables(fpoutput)
    dims = Flexpart.remdim(fpoutput, "spec001_mr")
    globattr = Flexpart.attrib(fpoutput)
    specattrib = Flexpart.attrib(fpoutput, "spec001_mr")
    adims = Flexpart.alldims(fpoutput, "spec001_mr")

    spec001 = Flexpart.select(fpoutput, "spec001_mr");
    oro = Flexpart.select(fpoutput, "ORO")

    spec001_alltimes = Flexpart.select(fpoutput, "spec001_mr", (time=:, height=1, pointspec=1, nageclass=1))
    spec001_2d_1 = Flexpart.select(fpoutput, "spec001_mr",
        Dict(:time=>adims[:time][4], :height=>adims[:height][2], :pointspec=>1, :nageclass=>1))

    spec001_2d_2 = Flexpart.select(fpoutput, "spec001_mr",
        Dict(:time=>adims[:time][7], :height=>adims[:height][2], :pointspec=>1, :nageclass=>1))

    spec001_2d_3 = Flexpart.select(fpoutput, "spec001_mr",
        Dict(:time=>adims[:time][7], :height=>:, :pointspec=>1, :nageclass=>1))

    seltime = Flexpart.selected(spec001_alltimes)
    sel = Flexpart.selected(spec001_2d)

    dayav_error = Flexpart.daily_average(spec001)
    dayav = Flexpart.daily_average(spec001_alltimes)

    compdim = Flexpart.completedim(spec001_alltimes)

    Flexpart.addable(spec001_2d_1, spec001_2d_2)
    Flexpart.addable(spec001_2d_1, spec001_alltimes)

    added = spec001_2d_1 + spec001_2d_2
    diff = spec001_2d_1 - spec001_2d_2

    Flexpart.write_daily_average(spec001, copy=true)

    ###################################
    ###### TEST FLEXEXTRACT ###########
    ###################################
    installpath = "/home/tcarion/flexpart/flex_extract_v7.1.2"
    defaultcontrol = "./test/fe_template/CONTROL_OD.OPER.FC.eta.highres.app"
    fepath = "test/fe_template/fedir"
    pythonpath = "/opt/anaconda3/bin/python3"

    fcontrol = FlexControl(defaultcontrol)
    area = [52.2, 4, 49, 6]
    fcontrol[:REQUEST] = 1
    # set!(fcontrol, Dict(:CLASS => "foo", :ETA => 2))
    set_area!(fcontrol, area)
    set_steps!(fcontrol, DateTime("2021-09-05T00:00:00"), DateTime("2021-09-07T00:00:00"), 1)

    fedir = FlexextractDir(fepath, fcontrol)
    fesource = FeSource(installpath, pythonpath)

    cmd = Flexpart.getcmd(fedir, fesource)

    # pip = pipeline(cmd, `sleep 3`, `echo COUCOU`)
    # logf = open("log.log", "w")
    # open(pip) do io
    #     lines = readlines(io, keep=true)
    #     for line in lines
    #         Base.write(logf, line)
    #         flush(logf)
    #     end
    # end
    # close(logf)

    # write(fcontrol, "test/fe_template/fe_output")

    ###################################
    ######## TEST MARSREQUESTS ########
    ###################################
    destmars = "test/fe_template/"
    csvpath = "test/fe_template/²fedir/input/mars_requests.csv"
    tmpdir = "test/tmp"
    requests = MarsRequest(csvpath)
    req = requests[1]
    push!(req.dict, :dqsdqsdq => "dqsdqsdq")
    Flexpart.retrieve(fesource, [req])
    # cmd = Flexpart.retrievecmd(fesource, req, tmpdir)

    open("log.log", "w") do logf
        Flexpart.retrieve(fesource, [req]) do stream
            data = readline(stream, keep=true)
            Base.write(logf, data)
            flush(logf)
        end
    end
    write(destmars, [requests[1]])

    ###################################
    ###### TEST PREPARE FLEXPART ######
    ###################################
    fedir = FlexextractDir(fepath)
    cmd = Flexpart.preparecmd(fedir, fesource)
    run(cmd)

    cmd = Flexpart.preparecmd(FlexextractDir(), fesource)
    Flexpart.prepare(fedir, fesource) do stream
        data = readline(stream, keep=true)
        Base.write(logf, data)
        flush(logf)
    end
    ###### TEST FLEXPART OLD ######
    close(logf)
    Flexpart.relloc(filename)
    Flexpart.start_dt(filename)
    Flexpart.times_dt(filename)
    output1 = Flexpart.FlexpartOutput(filename);
    Flexpart.select!(output1, "spec001_mr");
    output2 = Flexpart.FlexpartOutput(filename);
    @test Flexpart.hasselection(output1) == true
    @test Flexpart.hasselection(output2) == false
    @test Flexpart.hasfield2d(output1) == true

    Flexpart.isspatial(output1)
    Flexpart.isspatial(output2)
    Flexpart.remdim(output1)
    Flexpart.remdim(output2)

    dims = (time=:, height=1, pointspec=1, nageclass=1)
    dims = Dict(:height=>15., :time=>"20200203T080000", :pointspec=>1, :nageclass=>1)
    Flexpart.select!(output1, dims);
    time_av = Flexpart.faverage(output1);
    daily_av, days = Flexpart.daily_average(output1);
    Flexpart.attrib(output1)

    out_daily = Flexpart.write_daily_average(output1)
    Flexpart.select!(out_daily, "daily_av")
    Flexpart.select!(out_daily, (day=:,))
    out_worst = Flexpart.write_worst_day(out_daily)
    Flexpart.alldims(output2)
    Flexpart.close(output1)
    Flexpart.close(output2)

# end


text:
using Distributions
using Base.Test

d = Kolmogorov()
# compare to Smirnov's tabulated values
@test round(cdf(d,0.28),6) == .000_001
@test round(cdf(d,0.50),6) == .036_055
@test round(cdf(d,0.99),6) == .719_126
@test round(cdf(d,1.00),6) == .730_000
@test round(cdf(d,1.01),6) == .740_566
# @test round(cdf(d,1.04),6) == .770_434 # table value appears to be wrong
@test round(cdf(d,1.50),6) == .977_782
@test round(cdf(d,2.00),6) == .999_329
@test round(cdf(d,2.50),7) == .999_992_5
@test round(cdf(d,3.00),8) == .999_999_97


text:
module TestMagick
using Test

const msg_gz = [
    0x1F, 0x8B, 0x08, 0x08, 0x36, 0xBF, 0x68, 0x5E,
    0x00, 0x03, 0x6D, 0x73, 0x67, 0x2E, 0x74, 0x78,
    0x74, 0x00, 0xCB, 0x48, 0xCD, 0xC9, 0xC9, 0x57,
    0x48, 0x2F, 0x4A, 0x2C, 0xC8, 0xC8, 0x4C, 0xA9,
    0xE0, 0x02, 0x00, 0xCC, 0x85, 0x6F, 0x30, 0x0F,
    0x00, 0x00, 0x00
]

const msg_bz2 = [
    0x42, 0x5A, 0x68, 0x39, 0x31, 0x41, 0x59, 0x26,
    0x53, 0x59, 0x0D, 0x64, 0xB8, 0xC6, 0x00, 0x00,
    0x03, 0x51, 0x80, 0x00, 0x10, 0x40, 0x00, 0x26,
    0xE4, 0xD0, 0x40, 0x20, 0x00, 0x31, 0x00, 0xD0,
    0x01, 0x43, 0x46, 0x87, 0xEA, 0x9F, 0xA8, 0x88,
    0x4C, 0xB0, 0x36, 0x16, 0x85, 0xDC, 0x91, 0x4E,
    0x14, 0x24, 0x03, 0x59, 0x2E, 0x31, 0x80
]

const msg_txt = "hello graphidx\n"

is_gzip(io::IO)::Bool =
    read(io, 2) == [0x1F, 0x8B]

is_gzip(path::String)::Bool =
    is_gzip(open(path))

is_bzip2(io::IO)::Bool =
    read(io, 3) == UInt8['B', 'Z', 'h']

is_bzip2(path::String)::Bool =
    is_bzip2(open(path))

@testset "Magick Bytes                   " begin
    @test is_gzip(IOBuffer(msg_gz))
    @test !is_gzip(IOBuffer(msg_bz2))
    @test !is_gzip(IOBuffer(msg_txt))
    @test !is_bzip2(IOBuffer(msg_gz))
    @test is_bzip2(IOBuffer(msg_bz2))
    @test !is_bzip2(IOBuffer(msg_txt))
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
include(joinpath(dirname(@__FILE__), "Tree0" * (VERSION < v"0.6-" ? "5" : "6") * ".jl"))

Tree{T}(::Type{T}) = Tree(Node(Void, T, nothing))
AnyTree() = Tree(Any)
Tree() = AnyTree()

StructItems(items::Vector; idfield::Symbol = :id, parentfield::Symbol = :parent) = StructItems(items, idfield, parentfield)

KeyedTable{K, V}(data::Associative{K, V}, key) = Table{K, V}(data, e -> e[key])
StructTable{K, V}(data::Associative{K, V}, field::Symbol) = Table{K, V}(data, e -> getfield(e, field))

newnode{T}(::Tree{T}, data::T) = Node(T, T, data)
nodetype{T}(t::Tree{T}) = Node{T, T}

"""
  convert(Tree, items::Table{K, T})::Tree{T}

Converts a table to Tree.

# Examples

```jldoctest
julia> convert(Tree, KeyedTable(Dict(
                "ELECTRONICS" => ("ELECTRONICS", nothing),
                "TELEVISIONS" => ("TELEVISIONS", "ELECTRONICS"),
                "A" => ("A", nothing),
                "B" => ("B", "ELECTRONICS"),
              ), 2))
ImmutableNestedSetModel.Tree{Any}:
 1  2  ("A",nothing)
 3  8  ("ELECTRONICS",nothing)
 4  5  - ("B","ELECTRONICS")
 6  7  - ("TELEVISIONS","ELECTRONICS")

julia> type E
          id
          parent
       end

julia> convert(Tree, StructItems([
               E("ELECTRONICS", nothing),
               E("TELEVISIONS", "ELECTRONICS"),
               E("A", nothing),
               E("B", "ELECTRONICS"),
             ], :id, :parent))
ImmutableNestedSetModel.Tree{E}:
 1  2  E("A",nothing)
 3  8  E("ELECTRONICS",nothing)
 4  5  - E("B","ELECTRONICS")
 6  7  - E("TELEVISIONS","ELECTRONICS")

julia> convert(Tree, StructTable(Dict(
                  "ELECTRONICS" => E("ELECTRONICS", nothing),
                  "TELEVISIONS" => E("TELEVISIONS", "ELECTRONICS"),
                  "A" => E("A", nothing),
                  "B" => E("B", "ELECTRONICS"),
                ), :parent))
ImmutableNestedSetModel.Tree{E}:
 1  2  E("A",nothing)
 3  8  E("ELECTRONICS",nothing)
 4  5  - E("B","ELECTRONICS")
 6  7  - E("TELEVISIONS","ELECTRONICS")

```
"""
Base.convert{K, T}(::Type{Tree}, items::Table{K, T}) = prepare(append!(Tree(T), items))

"""
  convert(Tree, items::StructItems{T})::Tree{T}

Converts a StructItems to Tree.

# Examples

```jldoctest
julia> type E
          id
          parent
       end

julia> convert(Tree, StructItems([
               E("ELECTRONICS", nothing),
               E("TELEVISIONS", "ELECTRONICS"),
               E("A", nothing),
               E("B", "ELECTRONICS"),
             ], :id, :parent))
ImmutableNestedSetModel.Tree{E}:
 1  2  E("A",nothing)
 3  8  E("ELECTRONICS",nothing)
 4  5  - E("B","ELECTRONICS")
 6  7  - E("TELEVISIONS","ELECTRONICS")

```
"""
Base.convert{T}(::Type{Tree}, items::StructItems{T}) = prepare(append!(Tree(T), items))
Base.convert{T}(::Type{KeyedItems}, v::Vector{T}) = KeyedItems(v, 1, 2)

"""
  convert(Tree, items::Table{K, T})::Tree{T}

Convert default elements to Tree.

# Examples

```jldoctest
julia> convert(Tree, [
         ("ELECTRONICS", nothing),
         ("TELEVISIONS", "ELECTRONICS"),
         ("A", nothing),
         ("B", "ELECTRONICS"),
       ])
ImmutableNestedSetModel.Tree{Tuple{String,Any}}:
 1  2  ("A",nothing)
 3  8  ("ELECTRONICS",nothing)
 4  5  - ("B","ELECTRONICS")
 6  7  - ("TELEVISIONS","ELECTRONICS")

```
"""
Base.convert{T<:Union{Tuple, Vector, Associative}}(::Type{Tree}, v::Vector{T}) = prepare(append!(Tree(T), convert(KeyedItems, v)))

iter(t::Tree) = NodeIterator(t.root, true)

Base.push!{T}(t::Tree{T}, n::Node{T}) = push!(t.root, n)

function Base.append!{K, T}(t::Tree{T}, table::Table{K, T})
  nodes = Dict{Any, Node{T}}()

  for (id, item) = table.data
    nodes[id] = newnode(t, item)
  end

  for n = values(nodes)
    parentid = table.getparent(n.data)

    if parentid == nothing
      push!(t, n)
    else
      push!(nodes[parentid], n)
    end
  end
  t
end

function Base.append!{T}(t::Tree{T}, lines::StructItems{T})
  nodes = Dict{Any, nodetype(t)}()

  for item in lines.items
    nodes[getfield(item, lines.idfield)] = newnode(t, item)
  end

  for n in values(nodes)
    parentid = getfield(n.data, lines.parentfield)

    if parentid == nothing
      push!(t, n)
    else
      push!(nodes[parentid], n)
    end
  end
  t
end

function Base.append!{T}(t::Tree{T}, lines::KeyedItems)
  nodes = Dict{Any, nodetype(t)}()

  for item in lines.items
    nodes[item[lines.idkey]] = newnode(t, item)
  end

  for n in values(nodes)
    parentid = n.data[lines.parentkey]

    if parentid == nothing
      push!(t, n)
    else
      push!(nodes[parentid], n)
    end
  end
  t
end

function prepare(t::Tree)
  prepare(t.root, 0)
  t
end

tos(tree::Tree) = show(tree)

function Base.show{T}(io::IO, tree::Tree{T})
  write(io, repr(typeof(tree)), ":\n")
  for node = iter(tree)
    c = parentscount(node)
    write(io, @sprintf "%2d %2d %s %s\n" node.lft node.rgt (c == 0 ? "" : " " * repeat("-", c)) node.data)
  end
end

tov(tree::Tree) = String[tos(node) for node = iter(tree)]


text:
using Images

using PCG.Types
using PCG.Geometry
using PCG.Universes
using PCG.Topologies.HexGreedTopologies
using PCG.Neighborhoods
using PCG.Operations
using PCG.Spaces
using PCG.Recorders.GreedImages
using PCG.Recorders.HexGreedImages
using PCG.Recorders.TurnsLogger

const RADIUS = 20
const CELL_SIZE = Size(7, 7)
const DEBUG = false

struct Terrain <: Checkable
    value::Int64
end

function Operations.check(element::Element, parameters::Terrain)
    return element.properties.terrain == parameters
end

const GRASS = Terrain(1)
const WATER = Terrain(2)
const SAND = Terrain(3)
const FOREST = Terrain(4)

struct Properties <: AbstractProperties
    terrain::Terrain
end

const SPRITE_GRASS = HexSprite(RGBA{N0f8}(0.0, 1.0, 0.0, 1.0), CELL_SIZE)
const SPRITE_WATER = HexSprite(RGBA{N0f8}(0.0, 0.0, 1.0, 1.0), CELL_SIZE)
const SPRITE_SAND = HexSprite(RGBA{N0f8}(1.0, 1.0, 0.0, 1.0), CELL_SIZE)
const SPRITE_FOREST = HexSprite(RGBA{N0f8}(0.0, 0.5, 0.0, 1.0), CELL_SIZE)
const SPRITE_ERROR = HexSprite(RGBA{N0f8}(1.0, 0.0, 0.0, 1.0), CELL_SIZE)

function GreedImages.choose_sprite(recorder::GreedImageRecorder, element)
    if element |> GRASS |> exists
        return SPRITE_GRASS
    end

    if element |> WATER |> exists
        return SPRITE_WATER
    end

    if element |> SAND |> exists
        return SPRITE_SAND
    end

    if element |> FOREST |> exists
        return SPRITE_FOREST
    end

    return SPRITE_ERROR
end

function process(universe::Universe)

    ring = Neighborhood(universe.topology, ring_distance)
    euclidean = Neighborhood(universe.topology, euclidean_distance)

    complete_turn!(universe)

    universe() do element
        if element |> Fraction(0.01) |> exists
            element << (terrain=WATER,)
        end
    end

    universe() do element
        if element |> Fraction(0.8) |> GRASS |> euclidean(1, 5) |> WATER |> exists
            element << (terrain=WATER,)
        end
    end

    universe() do element
        if element |> GRASS |> ring() |> WATER |> exists
            element << (terrain=SAND,)
        end
    end

    universe(turns=3) do element
        if element |> Fraction(0.1) |> GRASS |> ring() |> SAND |> exists
            element << (terrain=SAND,)
        end
    end

    universe(turns=4) do element
        if element |> SAND |> ring() |> WATER |> count >= 5
            element << (terrain=WATER,)
        end
    end

    universe() do element
        if element |> GRASS |> Fraction(0.03) |> exists
            element << (terrain=FOREST,)
        end
    end

    universe() do element
        if (element |> GRASS |> Fraction(0.1) |> exists &&
            element |> ring(2, 2) |> FOREST |> exists &&
            element |> ring() |> new |> FOREST |> not_exists)

            element << (terrain=FOREST,)
        end
    end

end

topology = HexGreedTopology(RADIUS)

universe = initialize(topology,
                      Properties(GRASS),
                      DEBUG ? Recorder[] : [TurnsLoggerRecorder(),
                                            GreedImageRecorder(CELL_SIZE, 1000, "output.gif")])

@time process(universe)

finish_recording!(universe)

max_stars_repo_name: Tiendil/pcg-julia

max_issues_repo_name: Tiendil/pcg-julia

max_forks_repo_name: Tiendil/pcg-julia


text:
#function to calculate volume
function sphere_vol(r)

    return 4/3*pi*r^3
end

quadratic(a, sqr_term, b) = (-b + sqr_term) / 2a

#quadratic funtion calculation

function quadratic2(a::Float64, b::Float64, c::Float64)

    sqr_term = sqrt(b^2-4a*c)
    r1 = quadratic(a, sqr_term, b)
    r2 = quadratic(a, -sqr_term, b)

    r1, r2
end

vol = sphere_vol(3)

@printf "volume = %0.3f\n" vol
#prints volume = 113.097

quad1, quad2 = quadratic2(2.0, -2.0, -12.0)
println("result 1: ", quad1)
#> result 1: 3.0
println("result 2: ", quad2)
#> result 2: -2.0


text:
# convenience macro
macro wrap_h(signature)
  Meta.isexpr(signature, :call) ||
    throw(ArgumentError("signature has to be a function call expression"))

  name = signature.args[1]
  args = signature.args[2:end]
  args_wo_h = [arg for arg in args if arg !== :h]

  quote
    if f.$name === nothing
      nothing
    else
      if isinplace(f)
        let _f = f.$name, h = h
          ($(args_wo_h...),) -> _f($(args...))
        end
      else
        let _f = f.$name, h = h
          ($(args_wo_h[2:end]...),) -> _f($(args[2:end]...))
        end
      end
    end
  end |> esc
end

struct ODEFunctionWrapper{iip,F,H,TMM,Ta,Tt,TJ,JP,SP,TW,TWt,TPJ,S,TCV} <: DiffEqBase.AbstractODEFunction{iip}
  f::F
  h::H
  mass_matrix::TMM
  analytic::Ta
  tgrad::Tt
  jac::TJ
  jac_prototype::JP
  sparsity::SP
  Wfact::TW
  Wfact_t::TWt
  paramjac::TPJ
  syms::S
  colorvec::TCV
end

function ODEFunctionWrapper(f::DiffEqBase.AbstractDDEFunction, h)
  # wrap functions
  jac = @wrap_h jac(J, u, h, p, t)
  Wfact = @wrap_h Wfact(W, u, h, p, dtgamma, t)
  Wfact_t = @wrap_h Wfact_t(W, u, h, p, dtgamma, t)

  ODEFunctionWrapper{isinplace(f),typeof(f.f),typeof(h),typeof(f.mass_matrix),
                     typeof(f.analytic),typeof(f.tgrad),typeof(jac),
                     typeof(f.jac_prototype),typeof(f.sparsity),
                     typeof(Wfact),typeof(Wfact_t),
                     typeof(f.paramjac),typeof(f.syms),typeof(f.colorvec)}(
                       f.f, h, f.mass_matrix, f.analytic, f.tgrad, jac,
                       f.jac_prototype, f.sparsity, Wfact, Wfact_t, f.paramjac,
                       f.syms, f.colorvec)
end

(f::ODEFunctionWrapper{true})(du, u, p, t) = f.f(du, u, f.h, p, t)
(f::ODEFunctionWrapper{false})(u, p, t) = f.f(u, f.h, p, t)

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
function IncrementedValue(value, step)
    increment = value
    if increment == 0
        increment = 1.
    end
    myvalue = value + abs(increment) * (step)
    return myvalue
end

macro IterateParameters(ReferenceDict, VariedDict, ReadMethod, Steps, Model)
    quote
        for (key, value) in $(esc(ReferenceDict))
            if value[2] == "present"
                for (index, mystep) in enumerate($(esc(Steps)))
                    CopyDictCosmo = deepcopy($(esc(ReferenceDict)))
                    myvalue = IncrementedValue(value[1], mystep)
                    CopyDictCosmo[key] = [myvalue]
                    iterated = $(ReadMethod)(CopyDictCosmo, $(esc(Model)))
                    $(esc(VariedDict))["dvar_"*key*"_step_p_"*string(index)] = [iterated]
                    myvalue = IncrementedValue(value[1], -mystep)
                    CopyDictCosmo[key] = [myvalue]
                    iterated = $(ReadMethod)(CopyDictCosmo, $(esc(Model)))
                    $(esc(VariedDict))["dvar_"*key*"_step_m_"*string(index)] = [iterated]
                end
            end
        end
        iterated = $(ReadMethod)($(esc(ReferenceDict)), $(esc(Model)))
        $(esc(VariedDict))["dvar_central_step_0"] = [iterated]
    end
end

"""
    CreateCosmologies(DictCosmo::Dict, steps::Array)

This function creates a dictionary containing all combinations required for forecasts.
"""
function CreateCosmologies(DictCosmo::Dict, CosmoModel::String, IADict::Dict,
    IAModel::String, BiasDict::Dict, BiasModel::String, steps::Array)
    MyDictCosmo = Dict{String,Array{Any,1}}()
    @IterateParameters DictCosmo MyDictCosmo ReadCosmologyForecast steps CosmoModel

    MyDictIA = Dict{String,Array{Any,1}}()
    @IterateParameters IADict MyDictIA ReadIntrinsicAlignmentForecast steps IAModel

    MyDictBias = Dict{String,Array{Any,1}}()
    @IterateParameters BiasDict MyDictBias ReadBiasForecast steps BiasModel
    return MyDictCosmo, MyDictIA, MyDictBias
end

function CreateDirectoriesForecast!(Cosmologies::Dict, DictCosmo::Dict, IntrinsicAlignment::Dict,
    DictIA::Dict, Bias::Dict, DictBias::Dict, path::String)
    mkdir(path)
    mkdir(path*"PowerSpectrum")
    CreateDirectoriesPmmCℓ(Cosmologies, path, "PowerSpectrum", true)
    mkdir(path*"Angular")
    CreateDirectoriesPmmCℓ(Cosmologies, path, "Angular", true)
    CreateDirectoriesPmmCℓ(IntrinsicAlignment, path, "Angular", false)
    CreateDirectoriesPmmCℓ(Bias, path, "Angular", false)
    mkdir(path*"Derivative/")
    CreateDirectoriesDerivative(DictCosmo, path)
    CreateDirectoriesDerivative(DictIA, path)
    CreateDirectoriesDerivative(DictBias, path)
end

function CreateDirectoriesPmmCℓ(DictionaryParams::Dict, PathFolder::String, NameFolder::String,
    central::Bool)
    if central
        mkdir(PathFolder*"/"*NameFolder*"/"*"dvar_central_step_0")
    end
    for (key, value) in DictionaryParams
        if key != "dvar_central_step_0"
            mkdir(PathFolder*"/"*NameFolder*"/"*key)
        end
    end
end

function CreateDirectoriesDerivative(DictionaryParams::Dict, PathFolder::String)
    for (key, value) in DictionaryParams
        if value[2] == "present"
            mkdir(PathFolder*"Derivative/"*key)
        end
    end
end

function ForecastPowerSpectra!(Cosmologies::Dict, Path::String,
    CosmologicalGrid::CosmologicalGrid)
    for (key, value) in Cosmologies
        backgroundquantities = BackgroundQuantities(
        HZArray = zeros(length(CosmologicalGrid.ZArray)),
        rZArray=zeros(length(CosmologicalGrid.ZArray)))
        ComputeBackgroundQuantitiesGrid!(CosmologicalGrid,
        backgroundquantities, value[1])
        ClassyParams = Initializeclassy(value[1])
        powerspectrum = PowerSpectrum(PowerSpectrumLinArray =
        zeros(length(CosmologicalGrid.KArray), length(CosmologicalGrid.ZArray)),
        PowerSpectrumNonlinArray = zeros(length(CosmologicalGrid.KArray),
        length(CosmologicalGrid.ZArray)),
        InterpolatedPowerSpectrum = zeros(length(
        CosmologicalGrid.ℓBinCenters), length(CosmologicalGrid.ZArray)))
        EvaluatePowerSpectrum!(ClassyParams, CosmologicalGrid, powerspectrum)
        WritePowerSpectrumBackground(powerspectrum, backgroundquantities,
        CosmologicalGrid, Path*key*"/p_mm")
        WriteCosmology!(value[1], Path*key)
    end
end

function InstantiateComputeWeightFunctionGrid(
    ConvolvedDensity::AbstractConvolvedDensity,
    Cosmology::AbstractCosmology, CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities,
    GCWeightFunction::GCWeightFunction)
    ComputeBiasGrid!(CosmologicalGrid, GCWeightFunction, ConvolvedDensity)
    ComputeWeightFunctionGrid!(GCWeightFunction, ConvolvedDensity,
    CosmologicalGrid, BackgroundQuantities, Cosmology)
    return GCWeightFunction
end

function InstantiateComputeWeightFunctionGrid(
    ConvolvedDensity::AbstractConvolvedDensity,
    Cosmology::AbstractCosmology,
    CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities,
    LensingFunction::WLWeightFunction)
    ComputeLensingEfficiencyGrid!(LensingFunction, ConvolvedDensity,
    CosmologicalGrid, BackgroundQuantities, Cosmology, CustomLensingEfficiency())
    ComputeIntrinsicAlignmentGrid!(CosmologicalGrid, LensingFunction, ConvolvedDensity,
    BackgroundQuantities, Cosmology)
    ComputeWeightFunctionGrid!(LensingFunction, ConvolvedDensity,
    CosmologicalGrid, BackgroundQuantities, Cosmology)
    return LensingFunction
end

function Forecast∂Cℓ!(DictCosmo::Dict, IADict::Dict, BiasDict::Dict, PathInput::String,
    PathConfig::String, Steps::Array)
    ProbesDict = JSON.parsefile(PathConfig)
    CoefficientsArray = GetProbesArray(ProbesDict)
    for Coefficient in CoefficientsArray
        CentralCosmologyCL = ReadCℓ(
        PathInput*"/Angular/dvar_central_step_0/cl", Coefficient)
        CℓArray = zeros(size(CentralCosmologyCL.CℓArray, 1),
        size(CentralCosmologyCL.CℓArray, 2),
        size(CentralCosmologyCL.CℓArray, 3),
        2*length(Steps)+1)
        DerivativeArray = similar(CentralCosmologyCL.CℓArray)
        CℓArray[:, :, :, length(Steps) + 1] .=
        CentralCosmologyCL.CℓArray[:,:,:]
        stepvalues = zeros(2*length(Steps)+1)
        ∂cℓ = zeros(size(CℓArray, 1),
        size(CℓArray, 2), size(CℓArray, 3))
        #TODO: those three loops are quite similar. Macro?
        for (key, value) in DictCosmo
            if value[2] == "present"
                stepvalues[length(Steps) + 1] = value[1]
                for (index, mystep) in enumerate(Steps)
                    cℓminus = ReadCℓ(
                    PathInput*"/Angular/dvar_"*key*"_step_m_"*string(index)*"/cl",
                    Coefficient)
                    cℓplus = ReadCℓ(
                    PathInput*"/Angular/dvar_"*key*"_step_p_"*string(index)*"/cl",
                    Coefficient)
                    CℓArray[:, :, :, length(Steps) + 1 - index] .=
                    cℓminus.CℓArray
                    CℓArray[:, :, :, length(Steps) + 1 + index] .=
                    cℓplus.CℓArray
                    stepvalues[length(Steps) + 1 - index] =
                    IncrementedValue(value[1], -mystep)
                    stepvalues[length(Steps) + 1 + index] =
                    IncrementedValue(value[1],  mystep)
                end
                for idx_a in 1:size(CentralCosmologyCL.CℓArray, 2)
                    for idx_b in 1:size(CentralCosmologyCL.CℓArray, 3)
                        for idx_l in 1:size(CentralCosmologyCL.CℓArray, 1)
                            y = CℓArray[idx_l, idx_a, idx_b, :]
                            der = SteMDerivative(stepvalues, y)
                            ∂cℓ[idx_l, idx_a, idx_b] = der
                        end
                    end
                end
                Write∂Cℓ!(∂cℓ,
                PathInput*"/Derivative/"*key*"/"*key, Coefficient)
            end
        end
        for (key, value) in IADict
            if value[2] == "present"
                stepvalues[length(Steps) + 1] = value[1]
                for (index, mystep) in enumerate(Steps)
                    cℓminus = ReadCℓ(
                    PathInput*"/Angular/dvar_"*key*"_step_m_"*string(index)*"/cl",
                    Coefficient)
                    cℓplus = ReadCℓ(
                    PathInput*"/Angular/dvar_"*key*"_step_p_"*string(index)*"/cl",
                    Coefficient)
                    CℓArray[:, :, :, length(Steps) + 1 - index] .=
                    cℓminus.CℓArray
                    CℓArray[:, :, :, length(Steps) + 1 + index] .=
                    cℓplus.CℓArray
                    stepvalues[length(Steps) + 1 - index] =
                    IncrementedValue(value[1], -mystep)
                    stepvalues[length(Steps) + 1 + index] =
                    IncrementedValue(value[1],  mystep)
                end
                for idx_a in 1:size(CentralCosmologyCL.CℓArray, 2)
                    for idx_b in 1:size(CentralCosmologyCL.CℓArray, 3)
                        for idx_l in 1:size(CentralCosmologyCL.CℓArray, 1)
                            y = CℓArray[idx_l, idx_a, idx_b, :]
                            der = SteMDerivative(stepvalues, y)
                            ∂cℓ[idx_l, idx_a, idx_b] = der
                        end
                    end
                end
                Write∂Cℓ!(∂cℓ,
                PathInput*"/Derivative/"*key*"/"*key, Coefficient)
            end
        end
        for (key, value) in BiasDict
            if value[2] == "present"
                stepvalues[length(Steps) + 1] = value[1]
                for (index, mystep) in enumerate(Steps)
                    cℓminus = ReadCℓ(
                    PathInput*"/Angular/dvar_"*key*"_step_m_"*string(index)*"/cl",
                    Coefficient)
                    cℓplus = ReadCℓ(
                    PathInput*"/Angular/dvar_"*key*"_step_p_"*string(index)*"/cl",
                    Coefficient)
                    CℓArray[:, :, :, length(Steps) + 1 - index] .=
                    cℓminus.CℓArray
                    CℓArray[:, :, :, length(Steps) + 1 + index] .=
                    cℓplus.CℓArray
                    stepvalues[length(Steps) + 1 - index] =
                    IncrementedValue(value[1], -mystep)
                    stepvalues[length(Steps) + 1 + index] =
                    IncrementedValue(value[1],  mystep)
                end
                for idx_a in 1:size(CentralCosmologyCL.CℓArray, 2)
                    for idx_b in 1:size(CentralCosmologyCL.CℓArray, 3)
                        for idx_l in 1:size(CentralCosmologyCL.CℓArray, 1)
                            y = CℓArray[idx_l, idx_a, idx_b, :]
                            der = SteMDerivative(stepvalues, y)
                            ∂cℓ[idx_l, idx_a, idx_b] = der
                        end
                    end
                end
                Write∂Cℓ!(∂cℓ,
                PathInput*"/Derivative/"*key*"/"*key, Coefficient)
            end
        end
    end
end

#TODO: maybe we can create a macro?

function InstantiateWL(DictInput::Dict)
    LensingFunction = WLWeightFunction()
    return LensingFunction
end

function InstantiateWL(DictInput::Dict, IntrinsicAlignment::AbstractIntrinsicAlignment)
    LensingFunction = WLWeightFunction()
    LensingFunction.IntrinsicAlignmentModel = IntrinsicAlignment
    return LensingFunction
end

function InstantiateGC(DictInput::Dict)
    gcWeightFunction = GCWeightFunction()
    InstantiateBias(DictInput, gcWeightFunction)
    return gcWeightFunction
end

function InstantiateGC(DictInput::Dict, Bias::AbstractBias)
    gcWeightFunction = GCWeightFunction()
    gcWeightFunction.BiasKind = Bias
    return gcWeightFunction
end

function InstantiateBias(DictInput::Dict,
    GCWeightFunction::GCWeightFunction)
    if DictInput["PhotometricGalaxy"]["Bias"] == "PiecewiseBias"
        GCWeightFunction.BiasKind = PiecewiseBias()
    elseif DictInput["PhotometricGalaxy"]["Bias"] == "EuclidBias"
        GCWeightFunction.BiasKind = EuclidBias()
    else
        println("Bias must be correctly specified!")
    end
end

function InitializeProbes(DictInput::Dict,
    ConvolvedDensity::AbstractConvolvedDensity,
    Cosmology::AbstractCosmology,
    CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities)
    DictProbes = Dict()
    if DictInput["Lensing"]["present"]
        LensingFunction = InstantiateWL(DictInput::Dict)
        LensingFunction = InstantiateComputeWeightFunctionGrid(ConvolvedDensity,
        Cosmology, CosmologicalGrid, BackgroundQuantities,
        LensingFunction)
        push!(DictProbes, "Lensing" => LensingFunction)
    end
    if DictInput["PhotometricGalaxy"]["present"]
        GCWeightFunction = InstantiateGC(DictInput::Dict)
        GCWeightFunction =
        InstantiateComputeWeightFunctionGrid(ConvolvedDensity,
        Cosmology, CosmologicalGrid, BackgroundQuantities,
        GCWeightFunction)
        push!(DictProbes, "PhotometricGalaxy" => GCWeightFunction)
    end
    return DictProbes
end

function InitializeProbes(DictInput::Dict, ConvolvedDensity::AbstractConvolvedDensity,
    Cosmology::AbstractCosmology, IntrinsicAlignment::AbstractIntrinsicAlignment,
    Bias::AbstractBias, CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities, PathInput::String)
    DictProbes = Dict()
    if DictInput["Lensing"]["present"]
        LensingFunction = InstantiateWL(DictInput::Dict)
        LensingFunction.IntrinsicAlignmentModel = IntrinsicAlignment
        LensingFunction = InstantiateComputeWeightFunctionGrid(ConvolvedDensity,
        Cosmology, CosmologicalGrid, BackgroundQuantities,
        LensingFunction)
        push!(DictProbes, "Lensing" => LensingFunction)
    end
    if DictInput["PhotometricGalaxy"]["present"]
        GCWeightFunction = InstantiateGC(DictInput::Dict)
        GCWeightFunction.BiasKind = Bias
        GCWeightFunction = InstantiateComputeWeightFunctionGrid(ConvolvedDensity,
        Cosmology, CosmologicalGrid, BackgroundQuantities, GCWeightFunction)
        push!(DictProbes, "PhotometricGalaxy" => GCWeightFunction)
    end
    return DictProbes
end

function GetProbesArray(DictInput::Dict)
    ProbesArray = []
    CoefficientsArray = []
    if DictInput["Lensing"]["present"]
        push!(ProbesArray, "Lensing")
    end
    if DictInput["PhotometricGalaxy"]["present"]
        push!(ProbesArray, "PhotometricGalaxy")
    end
    sort!(ProbesArray)
    for key_A in ProbesArray
        for key_B in ProbesArray
            if key_B*"_"*key_A in CoefficientsArray
            else
                push!(CoefficientsArray, key_A*"_"*key_B)
            end
        end
    end
    return CoefficientsArray
end

function InitializeForecastCℓ(ProbesDict::Dict,
    BackgroundQuantities::BackgroundQuantities,
    Cosmology::AbstractCosmology,
    CosmologicalGrid::CosmologicalGrid, PowerSpectrum::PowerSpectrum,
    PathOutput::String, key::String)
    ProbesArray = []
    CoefficientsArray = []
    for (key, value) in ProbesDict
        push!(ProbesArray, key)
    end
    sort!(ProbesArray)
    for key_A in ProbesArray
        for key_B in ProbesArray
            if key_B*"_"*key_A in CoefficientsArray
            else
                push!(CoefficientsArray, key_A*"_"*key_B)
                cℓ = Cℓ(CℓArray =
                zeros(length(CosmologicalGrid.ℓBinCenters),
                length(ProbesDict[key_A].WeightFunctionArray[:, 1]),
                length(ProbesDict[key_B].WeightFunctionArray[:, 1])))
                ComputeCℓ!(cℓ,
                ProbesDict[key_A], ProbesDict[key_B], BackgroundQuantities,
                Cosmology, CosmologicalGrid, PowerSpectrum,
                CosmoCentral.CustomSimpson())
                WriteCℓ!(key_A*"_"*key_B,
                cℓ, PathOutput*key*"/cl")
                WriteCosmology!(Cosmology, PathOutput*key)
            end
        end
        WriteWeightFunctions!(key_A, ProbesDict[key_A], PathOutput*key*"/cl")
    end
    WriteCosmologicalGrid!(PathOutput*key*"/cl", CosmologicalGrid)
end

function InitializeForecastCℓ(ProbesDict::Dict,
    BackgroundQuantities::BackgroundQuantities,
    Cosmology::AbstractCosmology,
    CosmologicalGrid::CosmologicalGrid, PowerSpectrum::PowerSpectrum,
    PathOutput::String, CosmoDict::Dict, key::String)
    ProbesArray = []
    CoefficientsArray = []
    for (key, value) in ProbesDict
        push!(ProbesArray, key)
    end
    sort!(ProbesArray)
    for key_A in ProbesArray
        for key_B in ProbesArray
            #TODO is this if-else necessary?
            if key_B*"_"*key_A in CoefficientsArray
            else
                push!(CoefficientsArray, key_A*"_"*key_B)
                cℓ = Cℓ(CℓArray = zeros(length(CosmologicalGrid.ℓBinCenters),
                length(ProbesDict[key_A].WeightFunctionArray[:, 1]),
                length(ProbesDict[key_B].WeightFunctionArray[:, 1])))
                ComputeCℓ!(cℓ, ProbesDict[key_A], ProbesDict[key_B], BackgroundQuantities,
                Cosmology, CosmologicalGrid, PowerSpectrum,
                CosmoCentral.CustomSimpson())
                WriteCℓ!(key_A*"_"*key_B, cℓ, PathOutput*key*"/cl")
                WriteParameters!(CosmoDict, PathOutput*key)
            end
        end
    end
end

function ForecastCℓ!(Cosmologies::Dict, IntrinsicAlignment::Dict, Bias::Dict,
    PathInputPmm::String, PathOutputCℓ::String, CosmologicalGrid::CosmologicalGrid,
    PathConfigCℓ::String, PathInputCℓ::String)
    ProbesDict = JSON.parsefile(PathConfigCℓ)
    analyticaldensity = AnalitycalDensity()
    NormalizeAnalitycalDensity!(analyticaldensity)
    instrumentresponse = InstrumentResponse()
    convolveddensity = ConvolvedDensity(DensityGridArray =
    ones(10, length(CosmologicalGrid.ZArray)))
    NormalizeConvolvedDensity!(convolveddensity, analyticaldensity,
    instrumentresponse, CosmologicalGrid)
    ComputeConvolvedDensityGrid!(CosmologicalGrid, convolveddensity,
    analyticaldensity, instrumentresponse)
    #TODO: these three for loops are quite similar. The only difference is the cycled
    #dictionary. Maybe we can create a macro? The only difficulty is that we are iterating
    #over different dictionaries. maybe we can write some if in the dficitonaries to decide
    #which expr return?
    for (key, value) in Cosmologies
        Cosmology = value[1]
        IA = IntrinsicAlignment["dvar_central_step_0"][1]
        bias = Bias["dvar_central_step_0"][1]
        PowerSpectrum, BackgroundQuantities, CosmologicalGrid =
        ReadPowerSpectrumBackground(PathInputPmm*key*"/p_mm",
        CosmologicalGrid.ℓBinCenters, CosmologicalGrid.ℓBinWidths)
        ExtractGrowthFactor!(BackgroundQuantities, PowerSpectrum)
        CosmologicalGrid.ℓBinCenters[1,1] #TODO wtf represents this???
        DictProbes = InitializeProbes(ProbesDict, convolveddensity, Cosmology, IA,
        bias, CosmologicalGrid, BackgroundQuantities, PathInputCℓ)
        ComputeLimberArray!(CosmologicalGrid, BackgroundQuantities)
        InterpolatePowerSpectrumLimberGrid!(CosmologicalGrid, BackgroundQuantities,
        PowerSpectrum, CosmoCentral.BSplineCubic())
        InitializeForecastCℓ(DictProbes, BackgroundQuantities, Cosmology,
        CosmologicalGrid, PowerSpectrum, PathOutputCℓ, key)
    end

    for (key, value) in IntrinsicAlignment
        if key != "dvar_central_step_0"
            Cosmology = Cosmologies["dvar_central_step_0"][1]
            IA = value[1]
            bias = Bias["dvar_central_step_0"][1]
            PowerSpectrum, BackgroundQuantities, CosmologicalGrid =
            ReadPowerSpectrumBackground(PathInputPmm*"dvar_central_step_0"*"/p_mm",
            CosmologicalGrid.ℓBinCenters, CosmologicalGrid.ℓBinWidths)
            ExtractGrowthFactor!(BackgroundQuantities, PowerSpectrum)
            CosmologicalGrid.ℓBinCenters[1,1]
            DictProbes = InitializeProbes(ProbesDict, convolveddensity,
            Cosmology, IA, bias, CosmologicalGrid, BackgroundQuantities,
            PathInputCℓ)
            ComputeLimberArray!(CosmologicalGrid, BackgroundQuantities)
            InterpolatePowerSpectrumLimberGrid!(CosmologicalGrid,
            BackgroundQuantities, PowerSpectrum, CosmoCentral.BSplineCubic())
            InitializeForecastCℓ(DictProbes, BackgroundQuantities,
            Cosmology, CosmologicalGrid, PowerSpectrum, PathOutputCℓ, key)
        end
    end

    for (key, value) in Bias
        if key != "dvar_central_step_0"
            Cosmology = Cosmologies["dvar_central_step_0"][1]
            IA = IntrinsicAlignment["dvar_central_step_0"][1]
            bias = value[1]
            PowerSpectrum, BackgroundQuantities, CosmologicalGrid =
            ReadPowerSpectrumBackground(PathInputPmm*"dvar_central_step_0"*"/p_mm",
            CosmologicalGrid.ℓBinCenters, CosmologicalGrid.ℓBinWidths)
            ExtractGrowthFactor!(BackgroundQuantities, PowerSpectrum)
            CosmologicalGrid.ℓBinCenters[1,1]
            DictProbes = InitializeProbes(ProbesDict, convolveddensity,
            Cosmology, IA, bias, CosmologicalGrid, BackgroundQuantities,
            PathInputCℓ)
            ComputeLimberArray!(CosmologicalGrid, BackgroundQuantities)
            InterpolatePowerSpectrumLimberGrid!(CosmologicalGrid,
            BackgroundQuantities, PowerSpectrum, CosmoCentral.BSplineCubic())
            InitializeForecastCℓ(DictProbes, BackgroundQuantities,
            Cosmology, CosmologicalGrid, PowerSpectrum, PathOutputCℓ, key)
        end
    end
end

function ExtractVariedParameters(InputList::Vector{Dict{String, Vector{Any}}})
    OutputList = []
    for myDict in InputList
        for (key, value) in myDict
            if value[2] == "present"
                push!(OutputList, key)
            end
        end
    end
    return OutputList
end

"""
    ForecastFisherαβ(PathCentralCℓ::String, Path∂Cℓ::String,
    InputList::Vector{Dict{String, Vector{Any}}}, CosmologicalGrid::CosmologicalGrid)

This function evaluate the Fisher Matrix according to the following formula:
```math
F_{\\alpha \\beta}=\\sum_{\\ell=\\ell_{\\min }}^{\\ell_{\\max }} \\operatorname{Tr}\\left
\\{[\\boldsymbol{\\Sigma}(\\ell)]^{-1} \\frac{\\partial \\mathbf{C}(\\ell)}{\\partial
\\alpha}[\\boldsymbol{\\Sigma}(\\ell)]^{-1} \\frac{\\partial \\mathbf{C}(\\ell)}{\\partial
\\beta}\\right\\},
```
where ``\\alpha`` and ``\\beta`` are parameters of the Fisher Matrix, ``\\ell`` are the
multipoles, ``\\Sigma`` is the Covariance Matrix of the  Field Approach and
``\\frac{\\partial \\mathbf{C}(\\ell)}{\\partial \\alpha}`` is the Matrix of the derivatives
of the ``C_\\ell`` wrt parameter ``\\alpha``.
"""
function ForecastFisherαβ(PathCentralCℓ::String, Path∂Cℓ::String,
    InputList::Vector{Dict{String, Vector{Any}}}, CosmologicalGrid::CosmologicalGrid)
    Fisher = Fisherαβ()
    VariedParameters = ExtractVariedParameters(InputList)
    Fisher.FisherMatrix = zeros(length(VariedParameters), length(VariedParameters))
    Fisher.ParametersList = VariedParameters
    Fisher.SelectedParametersList = VariedParameters
    AnalitycalDensity = CosmoCentral.AnalitycalDensity()
    NormalizeAnalitycalDensity!(AnalitycalDensity)
    InstrumentResponse = CosmoCentral.InstrumentResponse()
    ConvolvedDensity = CosmoCentral.ConvolvedDensity(DensityGridArray =
    ones(10, length(CosmologicalGrid.ZArray)))
    NormalizeConvolvedDensity!(ConvolvedDensity, AnalitycalDensity, InstrumentResponse,
    CosmologicalGrid)
    ComputeConvolvedDensityGrid!(CosmologicalGrid, ConvolvedDensity, AnalitycalDensity,
    InstrumentResponse)
    ComputeSurfaceDensityBins!(ConvolvedDensity, AnalitycalDensity)

    Cℓ = ReadCℓ(PathCentralCℓ, "Lensing_Lensing")
    #TODO now only LL, but this need to be more flexible...maybe list with probes?
    Cov = InstantiateEvaluateCovariance(Cℓ, ConvolvedDensity, CosmologicalGrid, "Lensing",
    "Lensing")
    EvaluateFisherMatrix!(VariedParameters, Fisher, Path∂Cℓ, Cov)
    SelectMatrixAndMarginalize!(VariedParameters, Fisher)
    return Fisher
end

"""
    ForecastFisherαβ(PathCentralCℓ::String, Path∂Cℓ::String,
    InputList::Vector{Dict{String, Vector{Any}}}, CosmologicalGrid::CosmologicalGrid,
    ciccio::String)

This function evaluate the Fisher Matrix according to the following formula:
```math
F_{\\alpha \\beta}=\\sum_{\\ell=\\ell_{\\min }}^{\\ell_{\\max }} \\operatorname{vecp}
\\left(\\frac{\\partial \\mathbf{C}(\\ell)}{\\partial \\alpha}\\right)^{T}
\\left(\\boldsymbol{\\Xi}(\\ell)\\right)^{-1} \\operatorname{vecp}\\left(\\frac{\\partial
\\mathbf{C}(\\ell)} {\\partial \\beta}\\right),
```
where ``\\alpha`` and ``\\beta`` are parameters of the Fisher Matrix, ``\\ell`` are the
multipoles, ``\\Xi`` is the Covariance Matrix of the  Field Approach
[`CℓCovariance`](@ref) and ``\\frac{\\partial \\mathbf{C}(\\ell)}{\\partial \\alpha}``
is the Matrix of the derivatives of the ``C_\\ell`` wrt parameter ``\\alpha``.
"""
function ForecastFisherαβ(PathCentralCℓ::String, Path∂Cℓ::String,
    InputList::Vector{Dict{String, Vector{Any}}}, CosmologicalGrid::CosmologicalGrid,
    ciccio::String)
    Fisher = Fisherαβ()
    VariedParameters = ExtractVariedParameters(InputList)
    Fisher.FisherMatrix = zeros(length(VariedParameters), length(VariedParameters))
    Fisher.ParametersList = VariedParameters
    Fisher.SelectedParametersList = VariedParameters
    #here we instantiate the density again to evaluate the noise. Maybe it could be better
    #if we evaluated the density once for all, we passed it to Cℓ evaluator and then to
    #Fisher, in order to be more safe.
    AnalitycalDensity = CosmoCentral.AnalitycalDensity()
    NormalizeAnalitycalDensity!(AnalitycalDensity)
    InstrumentResponse = CosmoCentral.InstrumentResponse()
    ConvolvedDensity = CosmoCentral.ConvolvedDensity(DensityGridArray =
    ones(10, length(CosmologicalGrid.ZArray)))
    NormalizeConvolvedDensity!(ConvolvedDensity, AnalitycalDensity, InstrumentResponse,
    CosmologicalGrid)
    ComputeConvolvedDensityGrid!(CosmologicalGrid, ConvolvedDensity, AnalitycalDensity,
    InstrumentResponse)
    ComputeSurfaceDensityBins!(ConvolvedDensity, AnalitycalDensity)

    Cℓ = ReadCℓ(PathCentralCℓ, "Lensing_Lensing")
    #TODO now only LL, but this need to be more flexible...maybe list with probes?
    Covaₗₘ = InstantiateEvaluateCovariance(Cℓ, ConvolvedDensity, CosmologicalGrid, "Lensing",
    "Lensing")
    CovCℓ = InstantiateEvaluateCovariance(Covaₗₘ)
    EvaluateFisherMatrix!(VariedParameters, Fisher, Path∂Cℓ, CovCℓ)
    SelectMatrixAndMarginalize!(VariedParameters, Fisher)
    return Fisher
end

function EvaluateFisherMatrix!(VariedParameters::Vector{}, Fisher::AbstractFisher,
    Path∂Cℓ::String, Cov::AbstractCovariance)
    for (indexα, Parα) in enumerate(VariedParameters)
        ∂Cℓα = Read∂Cℓ(Path∂Cℓ*"/"*Parα*"/"*Parα, "Lensing_Lensing")
        for (indexβ, Parβ) in enumerate(VariedParameters)
            ∂Cℓβ = Read∂Cℓ(Path∂Cℓ*"/"*Parβ*"/"*Parβ, "Lensing_Lensing")
            EvaluateFisherMatrixElement!(Fisher, Cov, ∂Cℓα, ∂Cℓβ, Parα, Parβ)
            Fisher.FisherMatrix[indexα, indexβ] = Fisher.FisherDict[Parα*"_"*Parβ]
        end
    end
end

function SelectMatrixAndMarginalize!(VariedParameters::Vector{}, Fisher::AbstractFisher)
    for (idx, Par) in enumerate(reverse(VariedParameters))
        reverse_index = length(VariedParameters)-idx+1
        if Fisher.FisherMatrix[reverse_index, reverse_index] == 0
            Fisher.FisherMatrix = Fisher.FisherMatrix[1:end .!= reverse_index,
            1:end .!= reverse_index]
            Fisher.SelectedParametersList =
            Fisher.SelectedParametersList[1:end .!= reverse_index]
        end
    end
    Fisher.CorrelationMatrix = inv(Fisher.FisherMatrix)
    for (idxα, Parα) in enumerate(Fisher.SelectedParametersList)
        Fisher.MarginalizedErrors[Parα] = sqrt(Fisher.CorrelationMatrix[idxα, idxα])
    end
end


text:
#############################################################################
# geom_mean_epicone.jl
# Constrains T to
#   A #_t B ⪯ T
# where:
#   * A #_t B is the t-weighted geometric mean of A and B:
#     A^{1/2} (A^{-1/2} B A^{-1/2})^t A^{1/2}
#     Parameter t should be in [-1, 0] or [1, 2].
#   * Constraints A ⪰ 0, B ⪰ 0 are added.
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#
# Note on fullhyp: GeomMeanEpiCone will always return a full epigraph cone
# (unlike GeomMeanHypoCone) and so this parameter is not really used.  It is
# here just for consistency with the GeomMeanHypoCone function.
#
#REFERENCE
#   Ported from CVXQUAD which is based on the paper: "Lieb's concavity
#   theorem, matrix geometric means and semidefinite optimization" by Hamza
#   Fawzi and James Saunderson (arXiv:1512.03401)
#############################################################################

struct GeomMeanEpiCone
    A::AbstractExpr
    B::AbstractExpr
    t::Rational
    size::Tuple{Int, Int}

    function GeomMeanEpiCone(A::AbstractExpr, B::AbstractExpr, t::Rational, fullhyp::Bool=true)
        if size(A) != size(B)
            throw(DimensionMismatch("A and B must be the same size"))
        end
        n = size(A)[1]
        if size(A) != (n, n)
            throw(DimensionMismatch("A and B must be square"))
        end
        if t < -1 || (t > 0 && t < 1) || t > 2
            throw(DomainError(t, "t must be in the range [-1, 0] or [1, 2]"))
        end
        return new(A, B, t, (n, n))
    end

    GeomMeanEpiCone(A::Value,        B::AbstractExpr, t::Rational, fullhyp::Bool=true) = GeomMeanEpiCone(Constant(A), B, t, fullhyp)
    GeomMeanEpiCone(A::AbstractExpr, B::Value,        t::Rational, fullhyp::Bool=true) = GeomMeanEpiCone(A, Constant(B), t, fullhyp)
    GeomMeanEpiCone(A::Value,        B::Value,        t::Rational, fullhyp::Bool=true) = GeomMeanEpiCone(Constant(A), Constant(B), t, fullhyp)

    GeomMeanEpiCone(A::AbstractExprOrValue, B::AbstractExprOrValue, t::Integer, fullhyp::Bool=true) = GeomMeanEpiCone(A, B, t//1, fullhyp)
end

struct GeomMeanEpiConeConstraint <: Constraint
    head::Symbol
    id_hash::UInt64
    T::AbstractExpr
    cone::GeomMeanEpiCone

    function GeomMeanEpiConeConstraint(T::AbstractExpr, cone::GeomMeanEpiCone)
        if size(T) != cone.size
            throw(DimensionMismatch("T must be size $(cone.size)"))
        end
        id_hash = hash((cone.A, cone.B, cone.t, :GeomMeanEpiCone))
        return new(:GeomMeanEpiCone, id_hash, T, cone)
    end

    GeomMeanEpiConeConstraint(T::Value, cone::GeomMeanEpiCone) = GeomMeanEpiConeConstraint(Constant(T), cone)
end

in(T, cone::GeomMeanEpiCone) = GeomMeanEpiConeConstraint(T, cone)

function AbstractTrees.children(constraint::GeomMeanEpiConeConstraint)
    return (constraint.T, constraint.cone.A, constraint.cone.B, "t=$(constraint.cone.t)")
end

# For t ∈ [-1, 0] ∪ [1, 2] the t-weighted matrix geometric mean is matrix convex (arxiv:1512.03401).
# So if A and B are convex sets, then T ⪰ A #_t B will be a convex set.
function vexity(constraint::GeomMeanEpiConeConstraint)
    A = vexity(constraint.cone.A)
    B = vexity(constraint.cone.B)
    T = vexity(constraint.T)

    # NOTE: can't say A == NotDcp() because the NotDcp constructor prints a warning message.
    if typeof(A) == ConcaveVexity || typeof(A) == NotDcp
        return NotDcp()
    end
    if typeof(B) == ConcaveVexity || typeof(B) == NotDcp
        return NotDcp()
    end
    # Copied from vexity(c::GtConstraint)
    vex = ConvexVexity() + (-T)
    if vex == ConcaveVexity()
        vex = NotDcp()
    end
    return vex
end

function conic_form!(constraint::GeomMeanEpiConeConstraint, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, constraint)
        A = constraint.cone.A
        B = constraint.cone.B
        t = constraint.cone.t
        T = constraint.T
        is_complex = sign(A) == ComplexSign() || sign(B) == ComplexSign() || sign(T) == ComplexSign()
        if is_complex
            make_temporary = () -> HermitianSemidefinite(size(A)[1])
        else
            make_temporary = () -> Semidefinite(size(A)[1])
        end

        Z = make_temporary()

        if t <= 0
            conic_form!([T A; A Z] ⪰ 0, unique_conic_forms)
            conic_form!(Z in GeomMeanHypoCone(A, B, -t, false), unique_conic_forms)
        else
            @assert t >= 1 # range of t checked in GeomMeanEpiCone constructor
            conic_form!([T B; B Z] ⪰ 0, unique_conic_forms)
            conic_form!(Z in GeomMeanHypoCone(A, B, 2-t, false), unique_conic_forms)
        end

        cache_conic_form!(unique_conic_forms, constraint, Array{Convex.ConicConstr,1}())
    end
    return get_conic_form(unique_conic_forms, constraint)
end