text:
function RandomMDP(ns::Int64, na::Int64, discount::Float64; rng::AbstractRNG=Random.GLOBAL_RNG)
    # random dynamics
    T = rand(rng, ns, na, ns)
    # normalize
    for i = 1:ns, j = 1:na
        T[:,j,i] /= sum(T[:,j,i])
    end
    # random rewards [-0.5, 0.5]
    R = rand(rng, ns, na) .- 0.5
    return TabularMDP(T, R, discount)
end
RandomMDP() = RandomMDP(100, 5, 0.9)

function RandomPOMDP(ns::Int64, na::Int64, no::Int64, discount::Float64; rng::AbstractRNG=Random.GLOBAL_RNG)
    # random dynamics
    T = rand(rng, ns, na, ns)
    # random observation model
    O = rand(rng, no, na, ns)
    # normalize
    for i = 1:ns, j = 1:na
        T[:,j,i] /= sum(T[:,j,i])
        O[:,j,i] /= sum(O[:,j,i])
    end
    # random rewards [-0.5, 0.5]
    R = rand(rng, ns, na) .- 0.5
    return TabularPOMDP(T, R, O, discount)
end
RandomPOMDP() = RandomPOMDP(100, 5, 20, 0.9)


text:
using DashHtmlComponents

function Syntax(children::Component; style = code_container, summary = "")
    if (summary != "")
        return(
            html_details(
                (
                    html_summary(children = summary),
                    html_div(children)
                ),
                open = true
            )
        )
    else
        return(
            dcc_markdown(
                children,
                style = style
            )
        )
    end
end

ext: jl

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
type FakeDriverModel <: DriverModel{FakeDriveAction} end

let
    roadway = get_test_roadway()
    trajdata = get_test_trajdata(roadway)
    veh = get(trajdata, 1, 1)

    model = FakeDriverModel()
    reset_hidden_state!(model)
    observe!(model, Scene(), roadway, 1)
    prime_with_history!(model, trajdata, roadway, 1, 2, 1)

    @test get_name(model) == "???"
    @test action_type(model) <: FakeDriveAction
    @test_throws ErrorException rand(model)
    @test_throws ErrorException pdf(model, FakeDriveAction())
    @test_throws ErrorException logpdf(model, FakeDriveAction())
end


text:
"""
# Method p45

Analysis of elasto-plastic beams or rigid-jointed frames using a 2-node
Frame structural element in 1, 2 or 3 dimensions.

### Constructors
```julia
p45(data)
```
### Arguments
```julia
* `data::Dict{Symbol, Any}`  : Dictionary containing all input data
```

### Required data dictionary keys
```julia
* struc_el::StructuralElement                          : Type of  structural element
* support::Array{Tuple{Int,Array{Int,1}},1}        : Fixed-displacements vector
* loaded_nodes::Array{Tuple{Int,Array{Float64,1}},1} : Node load vector
* properties::Vector{Float64}                          : Material properties
* x_coords::FloatRange{Float64}                        : x-coordinate vector
* dload::FloatRange{Float64}                           : load steps
```

### Optional additional data dictionary keys
```julia
* penalty = 1e20                 : Penalty used for fixed degrees of freedoms
* etype::Vector{Int}           : Element material vector if np_types > 1
* y_coords::FloatRange{Float64}  : y-coordinate vector (2D)
* z_coords::FloatRange{Float64}  : x-coordinate vector (3D)
* limit = 250                    : Iteration limit
* tol = 0.0001                   : Tolerance for iteration convergence
```

### Related help
```julia
?StructuralElement             : List of available structural element types
?Frame                         : Help on a Frame structural element
?FiniteElement                 : List finite element types
?Line                          : Help on Line finite element
```
"""
function p45(data::Dict{Symbol, Any})

  # Parse & check FEdict data

  if :struc_el in keys(data)
    struc_el = data[:struc_el]
  else
    println("No fin_el type specified.")
    return
  end

  ndim = struc_el.ndim
  nst = struc_el.nst

  # Add radial stress
  if ndim == 3 && typeof(struc_el) !== Frame && struc_el.axisymmetric
    nst = 4
  end

  fin_el = struc_el.fin_el
  @assert typeof(fin_el) <: FiniteElement

  if typeof(fin_el) == Line
    if typeof(struc_el) == Frame
      nels = struc_el.nels
      nn = struc_el.nn
    else
      (nels, nn) = mesh_size(fin_el, struc_el.nxe)
    end
  elseif typeof(fin_el) == Triangle || typeof(fin_el) == Quadrilateral
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye)
  elseif typeof(fin_el) == Hexahedron
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye, struc_el.nze)
  else
    println("$(typeof(fin_el)) is not a known finite element.")
    return
  end

  if ndim == 1
    nodof = 2
  else
    nodof = ndim == 2 ? 3 : 6     # Degrees of freedom per node
  end
  ndof = fin_el.nod * nodof      # Degrees of freedom per fin_el

  # Update penalty if specified in FEdict

  penalty = 1e20
  if :penalty in keys(data)
    penalty = data[:penalty]
  end

  # Allocate all arrays

  # Start with arrays to be initialized from FEdict

  if :properties in keys(data)
    prop = zeros(size(data[:properties], 1), size(data[:properties], 2))
    for i in 1:size(data[:properties], 1)
      prop[i, :] = data[:properties][i, :]
    end
  else
    println("No :properties key found in FEdict")
  end

  nf = ones(Int, nodof, nn)

  if :support in keys(data)
    for i in 1:size(data[:support], 1)
      nf[:, data[:support][i][1]] = data[:support][i][2]
    end
  end

  x_coords = zeros(nn)
  if :x_coords in keys(data)
    x_coords = data[:x_coords]
  end

  y_coords = zeros(nn)
  if :y_coords in keys(data)
    y_coords = data[:y_coords]
  end

  z_coords = zeros(nn)
  if :z_coords in keys(data)
    z_coords = data[:z_coords]
  end

  etype = ones(Int, nels)
  if :etype in keys(data)
    etype = data[:etype]
  end

  g_num = zeros(Int, fin_el.nod, nels)
  if :g_num in keys(data)
    g_num = data[:g_num]
  end

  gamma = zeros(nels)
  if :gamma in keys(data)
    gamma = data[:gamma]
  end

  # All other arrays

  points = zeros(struc_el.nip, ndim)
  g = zeros(Int, ndof)
  g_coord = zeros(ndim,nn)
  fun = zeros(fin_el.nod)
  coord = zeros(fin_el.nod, ndim)
  jac = zeros(ndim, ndim)
  der = zeros(ndim, fin_el.nod)
  deriv = zeros(ndim, fin_el.nod)
  bee = zeros(nst,ndof)
  km = zeros(ndof, ndof)
  mm = zeros(ndof, ndof)
  gm = zeros(ndof, ndof)
  kg = zeros(ndof, ndof)
  eld = zeros(ndof)
  weights = zeros(struc_el.nip)
  g_g = zeros(Int, ndof, nels)
  num = zeros(Int, fin_el.nod)
  actions = zeros(ndof, nels)
  displacements = zeros(size(nf, 1), ndim)
  gc = ones(ndim, ndim)
  dee = zeros(nst,nst)
  sigma = zeros(nst)
  axial = zeros(nels)

  # Set global coordinates

  g_coord[1,:] = data[:x_coords]
  if ndim > 1
    g_coord[2,:] = data[:y_coords]
  end
  if ndim > 2
    g_coord[3,:] = data[:z_coords]
  end

  PtFEM.formnf!(nodof, nn, nf)
  neq = maximum(nf)

  ell = zeros(nels)
  if :x_coords in keys(data)
    for i in 1:length(data[:x_coords])-1
      ell[i] = data[:x_coords][i+1] - data[:x_coords][i]
    end
  end

  for i in 1:nels
    num = g_num[:, i]
    PtFEM.num_to_g!(num, nf, g)
    g_g[:, i] = g
  end

  println("There are $(neq) equations.\n")

  gsm = spzeros(neq, neq)
  for i in 1:nels
    num = g_num[:, i]
    coord = g_coord[:, num]'
    km = PtFEM.rigid_jointed!(km, prop, gamma, etype, i, coord)
    g = g_g[:, i]
    PtFEM.fsparm!(gsm, g, km)
  end

  limit = 10
  if :limit in keys(data)
    limit = data[:limit]
  end

  tol = 1.0e-5
  if :tol in keys(data)
    tol = data[:tol]
  end

  incs = 0
  if :incs in keys(data)
    incs = data[:incs]
  else
    println("No increments for loads specified.")
    return
  end

  dload = zeros(incs)
  if :dload in keys(data)
    dload = data[:dload]
  end

  loads = zeros(neq+1)
  if :loaded_nodes in keys(data)
    for i in 1:size(data[:loaded_nodes], 1)
      loads[nf[:, data[:loaded_nodes][i][1]] .+ 1] = data[:loaded_nodes][i][2]
    end
  end

  inode = zeros(Int, size(data[:loaded_nodes], 1))
  ival = zeros(size(data[:loaded_nodes], 1), size(data[:loaded_nodes][1][2], 2))
  if :loaded_nodes in keys(data)
    for i in 1:size(data[:loaded_nodes], 1)
      inode[i] = data[:loaded_nodes][i][1]
      ival[i, :] = data[:loaded_nodes][i][2]
    end
  end

  fixed_freedoms = 0
  if :fixed_freedoms in keys(data)
    fixed_freedoms = size(data[:fixed_freedoms], 1)
  end
  no = zeros(Int, fixed_freedoms)
  node = zeros(Int, fixed_freedoms)
  sense = zeros(Int, fixed_freedoms)
  value = zeros(Float64, fixed_freedoms)
  if fixed_freedoms > 0
    for i in 1:fixed_freedoms
      node[i] = data[:fixed_freedoms][i][1]
      sense[i] = data[:fixed_freedoms][i][2]
      no[i] = nf[sense[i], node[i]]
      value[i] = data[:fixed_freedoms][i][3]
      gsm[no[i], no[i]] += penalty
      loads[no[i] + 1] = gsm[no[i], no[i]] .* value[i]
    end
  end

  cfgsm = cholesky(gsm)

  # Special arrays
  eldtot = zeros(neq+1)
  holdr = zeros(ndof, nels)
  react = zeros(ndof)
  action = zeros(ndof)

  local converged = true

  total_load = 0.0
  for iy in 1:incs
    total_load += dload[iy]
    println("\nLoad step: $(iy), load factor: $total_load")
    bdylds = zeros(neq+1)
    oldlds = zeros(neq+1)
    iters = 0
    while true
      iters += 1
      print(".")
      loads = zeros(neq+1)
      for i in 1:size(inode, 1)
        loads[nf[:, data[:loaded_nodes][i][1]] .+ 1] = dload[iy] * ival[i, :]
      end
      loads += bdylds
      bdylds = zeros(neq+1)
      loads[2:end] = cfgsm \ loads[2:end]
      converged = checon(loads, oldlds, tol)
      for iel in 1:nels
        num = g_num[:, iel]
        coord = g_coord[:, num]'
        g = g_g[:, iel]
        eld = loads[g .+ 1]
        km = rigid_jointed!(km, prop, gamma, etype, iel, coord)
        action = km * eld
        react = zeros(ndof)
        if limit !== 1
          hinge!(coord, holdr, action, react, prop, iel, etype, gamma)
          #@show react
          bdylds[g .+ 1] -= react
          bdylds[1] = 0.0
        end
        if iters == limit || converged
          holdr[:, iel] += react[:] + action[:]
        end
      end
      if iters == limit || converged
        break
      end
    end
    eldtot += loads
    println("\n   Node     Displacement(s) and Rotation(s)  (Iterations: $(iters))")
    for i in 1:size(inode, 1)
      println("    $(inode[i])  $(eldtot[nf[:, inode[i]] .+ 1])")
    end
    if iters == limit && limit !== 1
      break
    end
  end
  (inode[size(inode, 1)], eldtot[nf[:, inode[size(inode, 1)]+1]])
end

max_stars_repo_path: src/4 Static Equilibrium/p45.jl

max_issues_repo_path: src/4 Static Equilibrium/p45.jl

max_forks_repo_path: src/4 Static Equilibrium/p45.jl


text:
using Gadfly
using Colors

# inspirado no estilo de plot da BBC, documentado em: https://bbc.github.io/rcookbook/#how_to_create_bbc_style_graphics
# (repo: https://github.com/bbc/bbplot)
function estilo_plot()
    Theme(
        plot_padding=[4mm, 4mm, 9mm, 9mm],
        bar_spacing=1mm,
        boxplot_spacing=1mm,
        background_color=colorant"white",
        panel_fill=colorant"white",

        guide_title_position=:left,
        key_position=:top,
        key_title_font_size=12pt,
        key_title_color=colorant"#222222",
        key_label_font_size=12pt,
        key_label_color=colorant"#222222",
        key_max_columns=6,

        major_label_font_size=13pt,
        major_label_color=colorant"#222222",

        minor_label_font_size=11pt,
        minor_label_color=colorant"#222222",

        point_label_font_size=10pt,
        point_label_color=colorant"#222222",

        grid_color=colorant"#cbcbcb"
    )
end

max_stars_repo_name: llrt/eda-qsa-rfb-jl

max_issues_repo_name: llrt/eda-qsa-rfb-jl

max_forks_repo_name: llrt/eda-qsa-rfb-jl


text:
# Demo the usage of IPPCore functions

using IPPCore

# library version

println(ippversion())
println()

# status -> message

for i = 0:9
	msg = ippstatus_string(i)
	println("Status code $i ==> $msg")
end
println()

# cpu information

println("CPU: $(ippcpudescr())")
# println("CPU Freq: $(ippcpufreq()) MHz")
println()


text:
"""
Creates a default visualization for any value.
The defaults can be customized via the key word arguments and the style parameter.
The style can change the the look completely (e.g points displayed as lines, or particles),
while the key word arguments just alter the parameters of one visualization.
Always returns a context, which can be displayed on a window via view(::Context, [display]).
"""
visualize(@nospecialize(main), s::Symbol=:default; kw_args...) = visualize(main, Style{s}(), Dict{Symbol, Any}(kw_args))

function visualize(@nospecialize(main), @nospecialize(s), @nospecialize(data))
    data = _default(main, s, copy(data))
    @gen_defaults! data begin # make sure every object has these!
        model = Mat4f0(I)
    end
    return assemble_shader(data)
end


text:
function okb_DA_algo(
    prop_pref::Matrix{Int}, # preference of proposers
    resp_pref::Matrix{Int}, # preference of respondents
    capacities::Vector{Int}) # capacities of respondents

    prop_n = size(prop_pref, 2)
    resp_n = size(resp_pref, 2)

    prop_match = zeros(Int, prop_n)
    resp_match = zeros(Int, sum(capacities))

    indicator = Array(Int, resp_n + 1) # for making resp_match
    indicator[1] = 1
    for i in 1:resp_n
        indicator[i+1] = indicator[i] + capacities[i]
    end

    resp_rank = _prefs2ranks(resp_pref) # [1 0 2] -> [1 3 2] in order not to use find!()

    occupied = zeros(Int, resp_n) # Numbers of respondendts' occupied seats

    propose_count = ones(Int, prop_n)  # Array which means the number of times of propose

    single_prop = trues(prop_n) # Bool values which mean unmatched(0) or matched(1)

    while any(single_prop)
        for p in 1:prop_n
            if single_prop[p]
                r = prop_pref[propose_count[p], p] # p proposes r

                if r == 0  # Who is proposed is none
                    single_prop[p] = false

                elseif resp_rank[p, r] > resp_rank[prop_n + 1, r]
                    # Proposal fails

                elseif occupied[r] < capacities[r]
                    # some seats are vacant
                    prop_match[p] = r
                    resp_match[indicator[r] + occupied[r]] = p
                    single_prop[p] = false
                    occupied[r] += 1

                else
                    #all seats are occupied
                    worst_rank = resp_match[indicator[r]] # Who is the least preferred among the currently accepted
                    worst_index = 0
                    for i in 1 : capacities[r] - 1
                        compared = resp_match[indicator[r] + i]
                        if resp_rank[worst_rank, r] < resp_rank[compared, r]
                            worst_rank = compared
                            worst_index = i
                        end
                    end

                    if resp_rank[p, r] < resp_rank[worst_rank, r]
                        # preferred more than worst_rank
                        prop_match[p] = r
                        prop_match[worst_rank] = 0
                        resp_match[indicator[r] + worst_index] = p
                        single_prop[p] = false
                        single_prop[worst_rank] = true
                    end

                end
                propose_count[p] += 1
            end
        end
    end
    return prop_match, resp_match, indicator
end

# case of one to one
function okb_DA_algo(
    prop_pref::Matrix{Int}, # preference of proposers
    resp_pref::Matrix{Int}) # preference of respondents

    capacities = ones(Int,size(resp_pref, 2))
    prop_match, resp_match, indicator =
    okb_DA_algo(prop_pref, resp_pref, capacities)
    return prop_match, resp_match
end

#from oyamad's sourse code(https://github.com/oyamad/Matching.jl/blob/master/src/deferred_acceptance.jl)
function _prefs2ranks{T<:Integer}(prefs::Matrix{T})
    unmatched = 0
    ranks = similar(prefs)
    m, n = size(prefs)
    for j in 1:n
        for i in 1:m
            k = prefs[i, j]
            if k == unmatched
                ranks[end, j] = i
            else
                ranks[k, j] = i
            end
        end
    end
    return ranks
end

max_stars_repo_name: M-okb/DA-Matching

max_issues_repo_name: M-okb/DA-Matching

max_forks_repo_name: M-okb/DA-Matching


text:
################################################################################
# add_port
################################################################################

function add_port(component::Component, name, class; metadata = emptymeta())
    return add_port(component, Port(name, class; metadata = metadata))
end

function add_port(component::Component, name, class, number; metadata = emptymeta())
    for i in 0:(number - 1)
        portname = "$name[$(i)]"
        # Choose whether to iterate through metadata or not.
        if typeof(metadata) <: Dict
            port = Port(portname, class; metadata = metadata)
        else
            port = Port(portname, class; metadata = metadata[i + 1])
        end
        add_port(component, port)
    end
    return nothing
end

function add_port(component::Component, port::Port)
    if haskey(component.ports, port.name)
        error("Port: $(port.name) already exists in component $(component.name)")
    end
    return component.ports[port.name] = port
end

################################################################################
# add_child
################################################################################

function add_child(component::Component, child::Component, name::String, number)
    for i in 0:(number - 1)
        add_child(component, child, "$name[$i]")
    end
    return nothing
end

function add_child(component::Component, child::Component, name::String)
    if haskey(component.children, name)
        error("Component $(c.name) already has a child $name.")
    end
    component.children[name] = deepcopy(child)
    return nothing
end

function add_child(
    toplevel::TopLevel, child::Component, address::Address, name::String = string(address)
)
    if haskey(toplevel.children, name)
        error("TopLevel $(toplevel.name) already has a child named $name.")
    end
    if haskey(toplevel.address_to_child, address)
        error("""
              TopLevel $(toplevel.name) already has a child assigned to
              address $address.
              """)
    end

    toplevel.children[name] = deepcopy(child)
    # Maintain cross references between child names and addresses
    toplevel.child_to_address[name] = address
    toplevel.address_to_child[address] = name
    return nothing
end

################################################################################
# add_link
################################################################################

# Convenience functions allowing add_link to be called with
#   - String, vector of string, PortPath, vector of PortPath and have
#   everything just work correctly. Maps all of these to a vector of portpaths.
portpath_promote(s::Vector{Path{Port}}) = identity(s)
portpath_promote(s::Path{Port}) = [s]
portpath_promote(s::String) = [Path{Port}(s)]
portpath_promote(s::Vector{String}) = [Path{Port}(i) for i in s]

function check_directions(component, sources, sinks)
    for source in sources
        port = relative_port(component, source)
        if !checkclass(port, Source)
            error("$source is not a valid source port for component $(component.name).")
        end
    end
    for sink in sinks
        port = relative_port(component, sink)
        if !checkclass(port, Sink)
            error("$sink is not a valid sink port for component $(component.name).")
        end
    end
end

function add_link(
    component::AbstractComponent,
    src_any,
    dest_any,
    safe = false;
    metadata = emptymeta(),
    linkname = nothing,
)

    # Promote the passed sources and destinations to Vector{Path{Port}}
    sources = portpath_promote(src_any)
    sinks = portpath_promote(dest_any)

    # Check port directions
    check_directions(component, sources, sinks)

    # check if port is available for a new link
    for port in Iterators.flatten((sources, sinks))
        if !isfree(component, port)
            safe ? (return false) : error("$port already has a connection.")
        end
    end

    # Check if we're trying to connect ports beyond just immedate hierarchy.
    for port in Iterators.flatten((sources, sinks))
        if length(port) == 0 || length(port) > 2
            error("""
                Links must connect ports defined within a component or to IO
                ports of an immediate child.

                $port violates this rule!.
                """)
        end
    end

    # create link

    # provide default name
    if linkname == nothing
        linkname = "link[$(length(component.links)+1)]"
    end

    if haskey(component.links, linkname)
        if safe
            return false
        else
            error("Link name: $linkname already exists in component $(component.name)")
        end
    end

    newlink = Link(linkname, sources, sinks, metadata)
    # Create a key for this link and add it to the component.
    component.links[linkname] = newlink
    # Assign the link to all top level ports.
    for port in Iterators.flatten((sources, sinks))
        # Register the link in the portlink dictionary
        component.portlink[port] = newlink
    end
    return true
end

tautology(args...) = true

"""
Single rule for connecting ports at the [`TopLevel`](@ref)
"""
struct Offset
    "Offset to add to a source address to reach a destination address"
    offset::Address

    "Name of the source port to start a link at."
    source_port::String

    "Name of the destination port to end a link at."
    dest_port::String

    # Do conversions to correct types.
    function Offset(
        offset::Union{Address,Tuple}, source_port::AbstractString, dest_port::AbstractString
    )
        return new(Address(offset), string(source_port), string(dest_port))
    end
end

# Allow passing of iterators as a constructor.
Offset(A, B, C) = [Offset(a, b, c) for (a, b, c) in zip(A, B, C)]

"""
Global connection rule for connecting ports at the [`TopLevel`](@ref)
"""
struct ConnectionRule
    """
    `Vector{Offset}` - Collection of [`Offset`] rules to be applied to all
    source addresses that pass the filtering stage.
    """
    offsets::Vector{Offset}

    "`Function` - Filter for source addresses. Default: `true`"
    address_filter::Function
    "`Function` - Filter for source components. Default: `true`"
    source_filter::Function
    "`Function` - Filter for destination components. Default: `true`"
    dest_filter::Function
end

# Provide KeyWord alternative.
function ConnectionRule(
    offsets; address_filter = tautology, source_filter = tautology, dest_filter = tautology
)
    return ConnectionRule(offsets, address_filter, source_filter, dest_filter)
end

Base.iterate(c::ConnectionRule) = (c, nothing)
Base.iterate(c::ConnectionRule, ::Any) = nothing

"""
    connection_rule(toplevel, rule::ConnectionRule; metadata = emptymeta())

Apply `rule::ConnectionRule` to `toplevel`. Source addresses will
first be filtered by `rule.address_filter`. Then, for each filtered
address, the [`Component`] at that address will be passed to
`rule.source_filter`.

If the component passes, all elements in `rule.offsets` will be
applied, assuming the component at destination address passes
`rule.dest_filter`. A new link will then be created provided it is
safe to do so.

Method List
-----------
$(METHODLIST)
"""
function connection_rule end

function connection_rule(toplevel::TopLevel, rule::Vector{Offset}; kwargs...)
    return connection_rule(toplevel, ConnectionRule(rule); kwargs...)
end

function connection_rule(toplevel::TopLevel, rule::ConnectionRule; metadata = emptymeta())
    # Count variable for verification - reports the number of links created.
    count = 0

    #for src_address in setdiff(valid_addresses, invalid_addresses)
    for src_address in Iterators.filter(rule.address_filter, addresses(toplevel))
        # Get the source component
        src = toplevel[src_address]
        # Check if the source component fulfills the requirement.
        rule.source_filter(src) || continue

        # Apply all the offsets to the current address
        for offset_rule in rule.offsets
            offset = offset_rule.offset
            src_port = offset_rule.source_port
            dst_port = offset_rule.dest_port

            # Check destination address
            dst_address = src_address + offset
            isaddress(toplevel, dst_address) || continue

            # check destination component
            dst = toplevel[dst_address]
            rule.dest_filter(dst) || continue

            # check port existence
            if !(haskey(src.ports, src_port) && haskey(dst.ports, dst_port))
                continue
            end

            # Build the name for these ports at the top level. If they are
            # free - connect them.
            src_prefix = getname(toplevel, src_address)
            dst_prefix = getname(toplevel, dst_address)

            src_port_path = Path{Port}(src_prefix, src_port)
            dst_port_path = Path{Port}(dst_prefix, dst_port)
            if isfree(toplevel, src_port_path) && isfree(toplevel, dst_port_path)
                # level ports container. If not, initialize them.
                add_link(toplevel, src_port_path, [dst_port_path])
                count += 1
            end
        end
    end

    @debug "Connections made: $count"
    return nothing
end

##############################
#           MUXES
##############################
function build_mux(inputs, outputs; metadata = emptymeta())
    name = "mux_" * string(inputs) * "_" * string(outputs)
    component = Component(name; primitive = "mux", metadata = metadata)
    add_port(component, "in", Input, inputs; metadata = metadata)
    add_port(component, "out", Output, outputs; metadata = metadata)
    return component
end

################################################################################
# Documentation
################################################################################

@doc """
    add_port(component, name, class, [number]; metadata = emptymeta())

Add `number` ports with the given `name` and `class`. Ports names will be the
provided suffix with bracket-vector notation. If `number` is not given, the
port will be added directly.

For example, the function call `add_port(component, "test", "input", 3)` should
add 3 `input` ports to component `c` with names: `test[2], test[1], test[0]`.

If `metadata` is given, it will be assigned to each instantiated port. If
`metadata` is a vector with `length(metadata) == number`, than entries of
`metadata` will be sequentially assigned to each instantiated port.
""" add_port

@doc """
    add_child(component, child::Component, name, [number])

Add a deepcopy child component with the given instance `name` to a component. If
`number` is provided, vectorize instantiation names. Throw error if `name` is
already used as an instance name for a child.
""" add_child

@doc """
    add_link(component, src, dest; metadata = emptymeta(), linkname = "")

Construct a link with the given metadata from the source ports to the
destination ports.

Arguments `src` and `dst` may of type `String`, `Vector{String}`, `PortPath`,
or `Vector{PortPath}`. If keyword argument `linkname` is given, the instantiated
link will be assigned that name. Otherwise, a unique name for the link will be
generated.

An error is raised if:
* Port classes are incorrect for the direction of the link.
* Ports in `src` or `dest` are already assigned to a link.
* A link with the given name already exists in `c`.
""" add_link

@doc """
    build_mux(inputs, outputs; metadata = Dict{String,Any}())

Build a mux with the specified number of inputs and outputs. Inputs and outputs
will be named `in[0], in[1], … , in[inputs-1]` and outputs will be named
`out[0], out[1], … , out[outputs-1]`.

If `metdata` is supplied, the dictionary will be attached to the mux component
itself as well as all ports and links in the mux component.
""" build_mux


text:
# load files and functions
include("set_up.jl")
using PyPlot
using KernelDensity
using Distributions
# scritp

# true parameter valuse
#Κ = 0.3
#Γ = 0.9
#B = 1.
#c = 28.5
#d =  4.
#A = 0.01
#f = 0.
#g = 0.03
#power1 = 1.5
#power2 = 1.8
#sigma =  1.9

################################################################################
###        old data set                                                      ###
################################################################################

# load new data
# load data
Z = convert(Array,readtable("DWPSDE model/Results/sim data small problem/data_usedgp_training_7_par_lunarc_simdata_4_coressimdata.csv"))
Z_data = Z[:,1]

#Z_data = load_data()

PyPlot.figure()
PyPlot.plot(Z_data)

bins = 50
PyPlot.figure()
h1 = PyPlot.plt[:hist](Z_data,bins)

# load data

Z_data = load_data()

PyPlot.figure()
PyPlot.plot(Z_data')

bins = 100
PyPlot.figure()
h1 = PyPlot.plt[:hist](Z_data',bins)

# parameter values to test
Κ = 0.3
Γ = 0.9
B = 1.
c = 28.5
d =  4.
A = 0.01
A_sign = 1
f = 0.
g = 0.03
power1 = 1.5
power2 = 1.8
sigma =  1.9

# new parameter values
#=
B = 1.;
c = 28.5;
d =  4.2;
A = 0.0052;
f = 0.;
g = 0.0610;
power1 = 1.5;
power2 = 1.8;
sigma =  1.9;
=#

# two unknown parameters
theta = log([c,d]) # true values for the unknown parameters
theta_known = [Κ, Γ, A, A_sign, B, f, g, power1, power2, sigma] # set vector with known parameters

# 7 parameters
#theta = log([0.2415 1.0075 28.4738 3.8861 1.5036 1.5500 1.6472]) # true values for the unknown parameters
#theta_known = [A,B,f,g] # set vector with known parameters

(Z_sim, dt, diff_dt, X_thinned) = generate_data(theta, theta_known,1)

PyPlot.figure()
PyPlot.plot(X_thinned)

PyPlot.figure()
PyPlot.plot(Z_sim)

bins = 100
PyPlot.figure()
h1 = PyPlot.plt[:hist](Z_sim,bins)

# compare kernel densities for real data and simulated data
Z_data_m = zeros(2,max(length(Z_data), length(Z_sim)))

for i = 1:length(Z_data)
  Z_data_m[1,i] = Z_data[i]
end

for i = 1:length(Z_sim)
  Z_data_m[2,i] = Z_sim[i]
end

k_real_data = kde(Z_data_m[1,:])
k_sim_data = kde(Z_data_m[2,:])

PyPlot.figure()
PyPlot.plot(k_real_data.x,k_real_data.density, "b")
PyPlot.hold(true)
PyPlot.plot(k_sim_data.x,k_sim_data.density, "b--")

# save data to file
Z_out = zeros(length(Z_sim),1)

for i = 1:length(Z_out)
  Z_out[i] = Z_sim[i]
end

writetable("data_old_new_dt.csv", convert(DataFrame, Z_out))

################################################################################
###          new data set                                                    ###
################################################################################

# load data
file = open("new_data_set.txt")
data = readlines(file)
close(file)

Z_data = zeros(length(data)-1)

idx = 2:length(data)

for i in idx
  try
    Z_data[i-1] = readdlm(IOBuffer(data[i]), Float64)[2]
  catch
    Z_data[i-1] = parse(data[i][end-1-4:end-1])
  end
end

# linear transformation of data to obtain a scaling where it is easier to
# construct the dwp model
Z_data = 50*Z_data

# plot data

PyPlot.figure()
PyPlot.plot(1:length(Z_data), Z_data)

bins = 100
PyPlot.figure()
h1 = PyPlot.plt[:hist](Z_data,bins)

# thinned data

thinning = 100

idx_thinned = 1:thinning:length(Z_data)
Z_data_thinned = zeros(Float64, idx_thinned)

j = 0
for i in idx_thinned
  j = j + 1
  Z_data_thinned[j] = Z_data[i]
end

# plot thinned data
PyPlot.figure()
PyPlot.plot(1:length(Z_data_thinned), Z_data_thinned)

bins = 100
PyPlot.figure()
h1 = PyPlot.plt[:hist](Z_data_thinned,bins)

# simulated data

# parameter values to test
Κ = 0.5;
Γ = 0.9;
B = 1.;
c = 22.5;
d =  13;
A = -0.0025;
f = 0.;
g = 0.;
power1 = 1.3;
power2 = 1.3;
sigma =  2.6;

# two unknown parameters
theta = log([c,d]) # true values for the unknown parameters
theta_known = [Κ, Γ, A, B, f, g, power1, power2, sigma] # set vector with known parameters
#=
thinning_sim = 1

N = Int64(3.5e6/thinning_sim)
dt = 20*1e-9 # 0.005*thinning_sim#20*1e-9
dt_U = 20*1e-9 # 0.005*thinning_sim#20*1e-9
diff_dt = 1

X_sim = zeros(N+1)
U_sim = zeros(N+1)

X_sim[1] = 35
U_sim[1] = 0

dB = zeros(N+1)

for i = 1:length(dB)
  dB[i] = sqrt(dt)*randn()
end

j = 1 # index for U_sim
b_const = sqrt(2*sigma^2 / 2)

#  should also consider the cache memory!
for i = 1:N # numerical integration of the X process
  X_sim[i+1] = X_sim[i] - (A*X_sim[i] - f + (power2*abs(abs(c - B*X_sim[i])^power1/2 - d + g*X_sim[i])^(power2 - 1)*sign(abs(c - B*X_sim[i])^power1/2 - d + g*X_sim[i])*(g - (B*power1*abs(c - B*X_sim[i])^(power1 - 1)*sign(c - B*X_sim[i]))/2))/2)*dt + b_const*dB[i]
  if  mod(i, diff_dt) == 0
    U_sim[j+1] = rand(Normal( U_sim[j]*exp(-Κ*dt_U), sqrt( Γ^2*( 1 - exp( -2*Κ*dt_U ) ) ) ),1)[1]
    j = j + 1
  end
end

X_thinned = X_sim[1:diff_dt:end]
Z_sim = X_thinned + U_sim
=#

Z_sim, X_thinned = generate_data(theta, theta_known)

PyPlot.figure()
PyPlot.plot(X_thinned)

PyPlot.figure()
PyPlot.plot(Z_sim)

bins = 100
PyPlot.figure()
h1 = PyPlot.plt[:hist](Z_sim,bins)

# compare kernel densities for real data and simulated data
Z_data_m = zeros(2,max(length(Z_data_thinned), length(Z_sim)))

for i = 1:length(Z_data_thinned)
  Z_data_m[1,i] = Z_data_thinned[i]
end

for i = 1:length(Z_sim)
  Z_data_m[2,i] = Z_sim[i]
end

k_real_data = kde(Z_data_m[1,:])
k_sim_data = kde(Z_data_m[2,:])

PyPlot.figure()
PyPlot.plot(k_real_data.x,k_real_data.density, "b")
PyPlot.hold(true)
PyPlot.plot(k_sim_data.x,k_sim_data.density, "b--")

# save data to file
Z_out = zeros(length(Z_sim),1)

for i = 1:length(Z_out)
  Z_out[i] = Z_sim[i]
end

writetable("data_new.csv", convert(DataFrame, Z_out))

max_stars_repo_path: DWPSDE model/generate_data_script.jl

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_path: DWPSDE model/generate_data_script.jl

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_path: DWPSDE model/generate_data_script.jl

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
"""
    coherancy(model::TimeSeriesModel, ω)

Compute the coherancy of a given model at frequency `ω`.
"""
function coherancy(model::TimeSeriesModel{D}, ω) where {D}
    S = sdf(model, ω)
    return [S[i,j]/sqrt(S[i,i]*S[j,j]) for i = 1:D, j = 1:D]
end
coherancy(::TimeSeriesModel{1}, ω) = error("coherancy not applicable for univariate models.")

"""
    coherance(model::TimeSeriesModel, ω)

Compute the coherance of a given model at frequency `ω`.
"""
coherance(model::TimeSeriesModel, ω) = abs.(coherancy(model, ω))

"""
    groupdelay(model::TimeSeriesModel, ω)

Compute the groupdelay of a given model at frequency `ω`.
"""
groupdelay(model::TimeSeriesModel, ω) = angle.(coherancy(model, ω))


text:
using Documenter, IterativeSolvers

makedocs(
	modules = [IterativeSolvers],
	format = :html,
	doctest = false,
	clean = true,
	sitename = "IterativeSolvers.jl",
	pages = [
		"Home" => "index.md",
		"Getting started" => "getting_started.md",
		"Preconditioning" => "preconditioning.md",
		"Linear systems" => [
			"Conjugate Gradients" => "linear_systems/cg.md",
			"Chebyshev iteration" => "linear_systems/chebyshev.md",
			"MINRES" => "linear_systems/minres.md",
			"BiCGStab(l)" => "linear_systems/bicgstabl.md",
			"IDR(s)" => "linear_systems/idrs.md",
			"Restarted GMRES" => "linear_systems/gmres.md",
			"LSMR" => "linear_systems/lsmr.md",
			"LSQR" => "linear_systems/lsqr.md",
			"Stationary methods" => "linear_systems/stationary.md"
		],
		"Eigenproblems" => [
			"Power method" => "eigenproblems/power_method.md",
			"LOBPCG" => "eigenproblems/lobpcg.md"
		],
		"SVDL" => "svd/svdl.md",
		"The iterator approach" => "iterators.md",
		"About" => [
			"Contributing" => "about/CONTRIBUTING.md",
			"License" => "about/license.md",
	    ]
	]
)

deploydocs(
	repo = "github.com/JuliaMath/IterativeSolvers.jl.git",
	target = "build",
	osname = "linux",
	julia  = "0.7",
	deps = nothing,
	make = nothing,
)


text:
function I(n)
    mat = zeros(n, n)
    for i in 1:n
        mat[i, i] = 1.0
    end
    return mat
end

mean(x) = sum(x)/length(x)

function var(x)
    m = mean(x)
    xm = x .- m
    return sum((xm .* xm) / (length(x) - 1))
end

std(x) = var(x) |> sqrt

geomean(x::Array{<: Number, 1}) = exp(sum(log.(x))/length(x))

function cor(x, y)
    xmd= x .- mean(x)
    ymd =y .- mean(y)
    return (sum(xmd .* ymd) / sqrt(var(x) * var(y))) / (length(x) - 1.0)
end

function cor(x::Matrix{Float64})
    _, p = size(x)
    cormat = ones(p, p)
    for i in 1:p
        for j in i:p
            if i != j
                cr = cor(x[:,i], x[:,j])
                cormat[i, j] = cr
                cormat[j, i] = cr
            end
        end
    end
    return cormat
end

function euclidean(v1::Array{T1,1}, v2::Array{T2,1})::Float64 where {T1 <: Number,T2 <: Number}
    (v1 .- v2).^2.0 |> sum |> sqrt
end

function euclidean(v1::Array{T,1})::Float64 where T <: Number
    v2 = zeros(Float64, length(v1))
    return euclidean(v1, v2)
end

function normalize(v1::Array{T,1})::Array{Float64,1} where T <: Number
    return v1 ./ euclidean(v1)
end

function normalize(data::DataFrame)::DataFrame
    df = similar(data)
    _, p = size(df)
    for i in 1:p
        df[:, i] = normalize(data[:, i])
    end
    return df
end

function apply_columns(fs::Array{Function,1}, data::Union{DataFrame,Array{T,2}} where T <: Number)
    _, m = size(data)
    return [fs[i](data[:,i]) for i in 1:m]
end

function apply_columns(f::Function, data::Union{DataFrame,Array{T,2}} where T <: Number)
    return [f(c) for c in eachcol(data)]
end

function colmins(data::DataFrame)
    return apply_columns(minimum, data)
end

function colmaxs(data::DataFrame)
    return apply_columns(maximum, data)
end

function colsums(data)::Array{Number,1}
    return apply_columns(sum, data)
end

function colmeans(data)::Array{Number,1}
    return apply_columns(mean, data)
end

function apply_rows(f::Function, data::Union{DataFrame,Array{T,2}} where T <: Number)
    return [f(c) for c in eachrow(data)]
end

function rowmins(data::DataFrame)::Array{Float64,1}
    return apply_rows(minimum, data)
end

function rowmaxs(data::DataFrame)::Array{Float64,1}
    return apply_rows(maximum, data)
end

function rowsums(data)::Array{Float64,1}
    return apply_rows(sum, data)
end

function rowmeans(data)::Array{Float64,1}
    return apply_rows(mean, data)
end

function unitize(v::Array{T,1})::Array{Float64,1} where {T <: Number}
    return v ./ sum(v)
end

function Base.:*(w::Array{Float64,1}, data::DataFrame)::DataFrame
    newdf = copy(data)
    _, p = size(newdf)
    for i in 1:p
        newdf[!, i] = convert.(Float64, newdf[!,i])
        newdf[:, i] = w[i] .* data[:, i]
    end
    return newdf
end

function Base.:-(r1::DataFrameRow, r2::DataFrameRow)::Array{Number,1}
    # v1 = convert(Array{Float64,1}, r1)
    v1 = Vector{Float64}(r1)
    # v2 = convert(Array{Float64,1}, r2)
    v2 = Vector{Float64}(r2)
    return v1 .- v2
end

function Base.:-(r1::Array{T,1}, r2::DataFrameRow)::Array{T,1} where T <: Number
    # v2 = convert(Array{Float64,1}, r2)
    v2 = Vector{Float64}(r2)
    return r1 .- v2
end

function Base.:-(r1::DataFrameRow, r2::Array{T,1})::Array{T,1} where T <: Number
    v1 = convert(Array{Float64,1}, r1)
    return v1 .- r2
end

"""
    makeDecisionMatrix(mat; names)

    Create a DataFrame using a decision matrix optionally using column names.

# Arguments:
 - `mat::Array{T,2}`: Matrix of numbers in any set
 - `names::Union{Nothing,Array{String,1}}`: Column names. Default is nothing. If column names are not given, they are labelled as Crt 1, Crt 2, ..., etc.

"""
function makeDecisionMatrix(mat::Array{T,2}; names::Union{Nothing,Array{String,1}}=nothing)::DataFrame where {T <: Number}
    _, m = size(mat)
    df = DataFrame()
    for i in 1:m
        if names isa Nothing
            name = string("Crt", i)
        else
            name = names[i]
        end
        df[:,Symbol(name)] = convert(Array{Float64,1}, mat[:, i])
    end
    return df
end

function Base.minimum(df::DataFrame)
    df |> Matrix |> minimum
end

function reverseminmax(fns::Array{Function,1})::Array{Function,1}
    newfs = map(x -> if x == minimum maximum else minimum end, fns)
    return newfs
end

function makeminmax(fns::Array{K,1} where K)::Array{Function,1}
    return convert(Array{Function,1}, fns)
end

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
export fields

add_verbose_scope(:Fields)
add_assert_scope(:Fields)

add_verbose_scope(:FieldsNonFancy)

mutable struct FieldsTower
  field::AnticNumberField
  generators_of_automorphisms::Vector{NfToNfMor}
  subfields::Vector{NfToNfMor}
  ramified_primes::Vector{fmpz}

  #Solvable embedding problems for the extension
  #They are here to improve the conductor computation
  has_info::Bool
  imgs_autos::Vector{Main.ForeignGAP.MPtr}
  auts_for_conductors::Vector{Main.ForeignGAP.MPtr}
  proj_ext::Main.ForeignGAP.MPtr

  function FieldsTower(K::AnticNumberField, auts::Vector{NfToNfMor}, subfields::Vector{NfToNfMor})
    z = new()
    z.field = K
    z.generators_of_automorphisms = auts
    z.subfields = subfields
    z.has_info = false
    return z
  end

end

include("./brauer.jl")
include("./merge.jl")
include("./abelian_layer.jl")
include("./read_write.jl")
include("./conductors.jl")

Generic.degree(F::FieldsTower) = degree(F.field)
Hecke.maximal_order(F::FieldsTower) = maximal_order(F.field)

function ramified_primes(F::FieldsTower)
  if !isdefined(F, :ramified_primes)
    f = factor(discriminant(maximal_order(F.field)))
    F.ramified_primes = collect(keys(f.fac))
  end
  return F.ramified_primes
end

###############################################################################
#
#  From automorphisms to permutation groups
#
###############################################################################

function Base.push!(G::AbstractAlgebra.Generic.geobucket{T}, p::T) where {T <: AbstractAlgebra.MPolyElem}
   R = parent(p)
   i = max(1, ndigits(length(p), base=4))
   l = length(G.buckets)
   if length(G.buckets) < i
     resize!(G.buckets, i)
     for j in (l + 1):i
       G.buckets[j] = zero(R)
     end
   end
   G.buckets[i] = addeq!(G.buckets[i], p)
   while i <= G.len
      if length(G.buckets[i]) >= 4^i
         G.buckets[i + 1] = addeq!(G.buckets[i + 1], G.buckets[i])
         G.buckets[i] = R()
         i += 1
      end
      break
   end
   if i == G.len + 1
      Base.push!(G.buckets, R())
      G.len += 1
   end
end

function permutation_group(G::Vector{Hecke.NfRelNSToNfRelNSMor{nf_elem}})
  permutations = permutation_group1(G)
  return _perm_to_gap_grp(permutations)
end

function permutations(G::Array{Hecke.NfToNfMor, 1})
  K = domain(G[1])
  n = length(G)
  dK = degree(K)
  d = numerator(discriminant(K.pol))
  p = 11
  while mod(d, p) == 0
    p = next_prime(p)
  end
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)
  fmod = Rx(K.pol)

  pols = gfp_poly[x]
  gpol = Rx(G[1].prim_img)
  if gpol != x
    push!(pols, gpol)
    gpol = compose_mod(gpol, pols[2], fmod)
    while gpol != x
      push!(pols, gpol)
      gpol = compose_mod(gpol, pols[2], fmod)
    end
  end
  order = length(pols)

  for i in 2:n
    pi = Rx(G[i].prim_img)
    if !(pi in pols)
      previous_order = order
      order = order + 1
      push!(pols, pi)
      for j in 2:previous_order
        order = order + 1
        push!(pols, compose_mod(pols[j], pi, fmod))
      end
      if order == n
        break
      end
      rep_pos = previous_order + 1
      while rep_pos <= order
        for k in 1:i
          po = Rx(G[k].prim_img)
          att = compose_mod(pols[rep_pos], po, fmod)
          if !(att in pols)
            order = order + 1
            push!(pols, att)
            for j in 2:previous_order
              order = order + 1
              push!(pols, compose_mod(pols[j], att, fmod))
            end
            if order == dK
              break
            end
          end
          if order == dK
            break
          end
        end
        if order == dK
          break
        end
        rep_pos = rep_pos + previous_order
      end
    end
  end
  #Now, I have the images mod p
  Dcreation = Vector{Tuple{gfp_poly, Int}}(undef, length(pols))
  for i = 1:length(pols)
    Dcreation[i] = (pols[i], i)
  end
  permutations = Array{Array{Int, 1},1}(undef, n)
  for i = 1:n
    permutations[i] = Vector{Int}(undef, dK)
  end
  gen_pols = [Rx(x.prim_img) for x in G]
  D = Dict{gfp_poly, Int}(Dcreation)
  for s = 1:n
    for i = 1:length(pols)
      permutations[s][i] = D[Hecke.compose_mod(gen_pols[s], pols[i], fmod)]
    end
  end
  return permutations
end

function permutation_group(G::Array{Hecke.NfToNfMor, 1})
  return _perm_to_gap_grp(permutations(G))
end

function _from_autos_to_perm(G::Array{Hecke.NfToNfMor,1})

  K = domain(G[1])
  @assert degree(K) == length(G)
  n = length(G)
  #First, find a good prime
  p = 2
  d = numerator(discriminant(K.pol))
  while mod(d, p) == 0
    p = next_prime(p)
  end
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)
  fmod = Rx(K.pol)
  pols = Vector{Tuple{gfp_poly, Int}}(undef, n)
  for i = 1:n
    pols[i] = (Rx(G[i].prim_img), i)
  end
  D = Dict{gfp_poly, Int}(pols)
  permutations = Array{Array{Int, 1},1}(undef, n)
  for s = 1:n
    perm = Array{Int, 1}(undef, n)
    for i = 1:n
      perm[i] = D[Hecke.compose_mod(pols[i][1], pols[s][1], fmod)]
    end
    permutations[s] = perm
  end
  return permutations

end

function _perm_to_gap_grp(perm::Array{Array{Int, 1},1})
  g = Main.ForeignGAP.MPtr[]
  for x in perm
    z = _perm_to_gap_perm(x)
    push!(g, z)
  end
  return GAP.Globals.Group(GAP.julia_to_gap(g))
end

function _perm_to_gap_perm(x::Array{Int, 1})
  z = GAP.Globals.PermList(GAP.julia_to_gap(x))
  return z
end

function IdGroup(autos::Array{NfToNfMor, 1})
  G = permutation_group(autos)
  return GAP.Globals.IdGroup(G)
end

################################################################################
#
#  final computation of the maximal order and automorphisms
#
################################################################################

function _from_relative_to_abs_with_embedding(L::Hecke.NfRelNS{T}, autL::Array{Hecke.NfRelNSToNfRelNSMor{T}, 1}, use_simplify::Bool = true) where T

  S, mS = simple_extension(L)
  K, mK, MK = absolute_field(S, false)

  #First, we compute the maximal order of the absolute field.
  #We start from the maximal orders of the relative extension and of the base field.
  #FALSE: Since the computation of the relative maximal order is slow, I prefer to bring to the absolute field the elements
  # generating the equation order.
  pols = L.pol
  gL = gens(L)
  B = Array{nf_elem, 1}(undef, degree(K))
  B[1] = K(1)
  ind = total_degree(pols[1])
  genjj = mK\(mS\gL[1])
  for i = 2:ind
    B[i] = genjj*B[i-1]
  end
  for jj = 2:length(pols)
    genjj = mK\(mS\gL[jj])
    el = deepcopy(genjj)
    new_deg = total_degree(pols[jj])
    for i = 2:new_deg
      for j = 1:ind
        B[(i-1)* ind + j] = B[j]* el
      end
      mul!(el, el, genjj)
    end
    ind *= new_deg
  end

  #Now, I add the elements of the maximal order
  OB = maximal_order(base_field(S))
  for i = 1:degree(OB)
    el = MK(OB.basis_nf[i])
    for j = 1:ind
      B[(i-1)* ind + j] = B[j] * el
    end
  end
  @vprint :Fields 2 "Product basis computed\n"
  #Now, we compute the maximal order. Then we simplify.
  #We simplify only if the degree of the field is lower than 30

  BasisMat = basis_matrix(B, FakeFmpqMat)
  @vtime :Fields 3 Hecke.hnf_modular_eldiv!(BasisMat.num, BasisMat.den, :lowerleft)
  NewBMat = FakeFmpqMat(BasisMat.num, BasisMat.den)
  @vtime :Fields 3 Ostart = NfAbsOrd(K, NewBMat)
  Ostart.index = divexact(NewBMat.den^degree(K), prod_diagonal(NewBMat.num))
  Ostart.gen_index = fmpq(Ostart.index)
  Ostart.disc = divexact(numerator(discriminant(K)), Ostart.index^2)
  ram_primes_rel = numerator(norm(discriminant(L)))
  ram_primes_down = Hecke.ramified_primes(OB)
  for p in ram_primes_down
    if isone(gcd(p, ram_primes_rel))
      push!(Ostart.primesofmaximality, p)
    end
  end
  @vtime :Fields 3 O1 = MaximalOrder(Ostart)
  O1.ismaximal = 1
  Hecke._set_maximal_order_of_nf(K, O1)
  if use_simplify
    @vtime :Fields 3 Ks, mKs = Hecke.simplify(K)
    #Now, we have to construct the maximal order of this field.
    #I compute the inverse of mKs
    @vtime :Fields 3 mKsI = find_inverse(mKs)
    if isdefined(O1, :lllO)
      lO = O1.lllO::NfOrd
      O2 = NfOrd(nf_elem[mKsI(x) for x in basis(lO, K, copy = false)], false)
      #O2.lllO = O2
    else
      O2 = NfOrd(nf_elem[mKsI(x) for x in basis(O1, K, copy = false)], false)
    end
    O2.ismaximal = 1
    @assert isdefined(O1, :disc)
    O2.disc = O1.disc
    O2.index = root(divexact(numerator(discriminant(Ks)), O2.disc), 2)
    @vtime :Fields 3 OLLL = lll(O2)
    Hecke._set_maximal_order_of_nf(Ks, OLLL)
  else
    Ks = K
    mKs = id_hom(K)
    mKsI = id_hom(K)
  end
  #I want also the embedding of the old field into the new one.
  #It is enough to find the image of the primitive element.
  k = base_field(S)
  a = MK(gen(k))
  embed = NfToNfMor(k, Ks, mKsI(a))
  #@assert iszero(k.pol(img_a))
  @vprint :Fields 3 "MaximalOrder Computed. Now Automorphisms\n"
  #Now, the automorphisms.
  # I need both generators and the whole group.
  autos = Array{NfToNfMor, 1}(undef, length(autL))
  el = mS(mK(mKs.prim_img))
  el1 = mS(mK(gen(K)))
  for i=1:length(autL)
    #@assert iszero(K.pol(mK(mS\(autL[i](el1)))))
    x = mKsI(mK\(mS\(autL[i](el))))
    #@assert Ks.pol(y) == 0
    autos[i] = Hecke.NfToNfMor(Ks, Ks, x)
  end
  #And the generators are done. Now the closure
  #@vtime :Fields 3 Hecke._set_automorphisms_nf(Ks, closure(autos, degree(Ks)))
  #No! I set the automorphisms only if I need them
  #Hecke._set_automorphisms_nf(Ks, autsKs)
  @vprint :Fields 2 "Finished\n"
    #@assert codomain(embed) == Ks
  return Ks, autos, embed
end

function find_inverse(f::NfToNfMor)
  v = Vector{nf_elem}(undef, degree(domain(f)))
  v[1] = one(codomain(f))
  v[2] = f.prim_img
  for i = 3:degree(domain(f))
    v[i] = v[2]*v[i-1]
  end
  M = basis_matrix(v)
  w = zero_matrix(FlintQQ, 1, nrows(M))
  w[1, 2] = 1
  fl, s = can_solve(M, w, side = :left)
  @assert fl
  s1 = FakeFmpqMat(s)
  pre_img = elem_from_mat_row(domain(f), s1.num, 1, s1.den)
  return hom(codomain(f), domain(f), pre_img)
end

function find_inverse2(f::NfToNfMor)
  K = codomain(f)
  arr_prim_img = Array{nf_elem, 1}(undef, degree(K))
  arr_prim_img[1] = K(1)
  arr_prim_img[2] = f.prim_img
  for i = 3:degree(K)
    arr_prim_img[i] = arr_prim_img[i-1]*f.prim_img
  end
  M1 = inv(basis_matrix(arr_prim_img, FakeFmpqMat))
  return M1
end

###############################################################################
#
#  Split Extensions
#
###############################################################################

function _split_extension(G::Array{Hecke.NfToNfMor, 1}, mats::Array{Hecke.GrpAbFinGenMap, 1})

  gtype = map(Int, domain(mats[1]).snf)
  G1 = permutation_group(G)
  gensG1 = GAP.Globals.GeneratorsOfGroup(G1)
  A = GAP.Globals.AbelianGroup(GAP.julia_to_gap(gtype))
  gens = GAP.Globals.GeneratorsOfGroup(A)
  auts = Array{Main.ForeignGAP.MPtr, 1}(undef, length(mats))
  for i = 1:length(mats)
    images = Array{Main.ForeignGAP.MPtr, 1}(undef, length(gtype))
    for j = 1:length(gtype)
      g = GAP.Globals.Identity(A)
      for k = 1:length(gtype)
        if !iszero(mats[i].map[j,k])
          g *= gens[k]^Int(mats[i].map[j,k])
        end
      end
      images[j] = g
    end
    auts[i] = GAP.Globals.GroupHomomorphismByImages(A, A, gens, GAP.julia_to_gap(images))
  end
  AutGrp = GAP.Globals.Group(GAP.julia_to_gap(auts))
  mp = GAP.Globals.GroupHomomorphismByImages(G1, AutGrp, gensG1, GAP.julia_to_gap(auts))
  return GAP.Globals.SplitExtension(G1, mp, A)

end

###############################################################################
#
#  Check group extension
#
###############################################################################

function check_group_extension(TargetGroup::Main.ForeignGAP.MPtr, autos::Array{NfToNfMor, 1}, res_act::Array{GrpAbFinGenMap, 1})

  GS = domain(res_act[1])
  expo = Int(GS.snf[end])
  K = domain(autos[1])
  d = degree(K)
  com, uncom = ppio(expo, d)

  if com == 1
    # I only need to check the split extension, since the second cohomology group is
    # trivial, regardless of the action
    H = _split_extension(autos, res_act)
    return GAP.Globals.IdGroup(H) == TargetGroup
  end

  if uncom == 1
    #Need a cohomological check. Only useful in the prime power case.
    return true
  end

  # I check the split extension related to only uncom
  #Now, I have to check if the split extension is isomorphic to IdH
  Qn, mQn = quo(GS, uncom, false)
  S1, mS1 = snf(Qn)
  new_res_act = Array{GrpAbFinGenMap, 1}(undef, length(res_act))
  for i = 1:length(res_act)
    Mat = mS1.map*mQn.imap*res_act[i].map*mQn.map*mS1.imap
    Hecke.reduce_mod_snf!(Mat, S1.snf)
    new_res_act[i] = hom(S1, S1, Mat)
  end
  H = _split_extension(autos, new_res_act)
  return GAP.Globals.IdGroup(H) == TargetGroup

end

###############################################################################
#
#  Interface to find Fields
#
###############################################################################

function field_extensions(list::Vector{FieldsTower}, bound::fmpz, IsoE1::Main.ForeignGAP.MPtr, l::Array{Int, 1}, only_real::Bool, simplify::Bool)

  grp_to_be_checked = Dict{Int, Main.ForeignGAP.MPtr}()
  d = degree(list[1])
  n = lcm(l)
  com, uncom = ppio(n, d)
  fac = factor(n)
  for (p, v) in fac
    grp_to_be_checked[Int(p)] = _construct_grp(IsoE1, Int(p)^v)
  end
  if uncom != 1
    IsoCheck = _construct_grp(IsoE1, uncom)
  else
    IsoCheck = IsoE1
  end
  final_list = FieldsTower[]
  for (j, x) in enumerate(list)
    @vprint :Fields 1 "Field $(j)/$(length(list)): $(x.field.pol)"
    @vprint :FieldsNonFancy 1 "Field $(j)/$(length(list)): $(x.field.pol)\n"
    append!(final_list, field_extensions(x, bound, IsoCheck, l, only_real, grp_to_be_checked, IsoE1, simplify))
  end
  return final_list

end

function field_extensions(x::FieldsTower, bound::fmpz, IsoE1::Main.ForeignGAP.MPtr, l::Array{Int, 1}, only_real::Bool, grp_to_be_checked::Dict{Int, Main.ForeignGAP.MPtr}, IsoG::Main.ForeignGAP.MPtr, simplify::Bool)

  list_cfields = _abelian_normal_extensions(x, l, bound, IsoE1, only_real, IsoG)
  if isempty(list_cfields)
    @vprint :Fields 1 "\e[1F$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())Number of new fields found: 0\n\n"
    @vprint :FieldsNonFancy 1 "Number of new fields found: 0\n"
    return Vector{FieldsTower}()
  end
  list = from_class_fields_to_fields(list_cfields, x.generators_of_automorphisms, grp_to_be_checked)
  @vprint :Fields 1 "Computing maximal orders"
  @vprint :FieldsNonFancy 1 "Computing maximal orders\n"
  final_list = Vector{FieldsTower}(undef, length(list))
  for j = 1:length(list)
    fld, autos, embed = _from_relative_to_abs_with_embedding(list[j][1], list[j][2], simplify)
    previous_fields = Array{NfToNfMor, 1}(undef, length(x.subfields)+1)
    for s = 1:length(x.subfields)
      previous_fields[s] = x.subfields[s]
    end
    previous_fields[end] = embed
    final_list[j] = FieldsTower(fld, autos, previous_fields)
  end
  @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
  @vprint :Fields 1 "Number of new fields found: $(length(final_list))\n\n"
  @vprint :FieldsNonFancy 1 "Number of new fields found: $(length(final_list))\n"
  return final_list

end

###############################################################################
#
#  Interface
#
###############################################################################

function fields_direct_product(g1, g2, red, redfirst, absolute_bound; only_real = false)
  b1 = root(absolute_bound, g2[1])
  b2 = root(absolute_bound, g1[1])
  l2 = fields(g2[1], g2[2], b2, only_real = only_real)
  if isempty(l2)
    return FieldsTower[]
  end
  if g1 == g2
    return _merge(l2, l2, absolute_bound, red, redfirst)
  end
  l1 = fields(g1[1], g1[2], b1, only_real = only_real)
  if isempty(l1)
    return FieldsTower[]
  end
  return _merge(l1, l2, absolute_bound, red, redfirst)
end

function fields(a::Int, b::Int, absolute_bound::fmpz; using_direct_product::Bool = true, only_real::Bool = false, simplify::Bool = true)
  @assert absolute_bound > 0
  if a == 1
    @assert b == 1
    Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
    K, a = NumberField(x-1, cached = false)
    g = NfToNfMor(K, K, K(1))
    return FieldsTower[FieldsTower(K, NfToNfMor[g], Array{NfToNfMor, 1}())]
  end
  G = GAP.Globals.SmallGroup(a, b)
  @assert GAP.Globals.IsSolvable(G)
  if using_direct_product
    g1, g2, red, redfirst = direct_product_decomposition(G, (a, b))
    if g2 != (1, 1)
      return fields_direct_product(g1, g2, red, redfirst, absolute_bound; only_real = false)
    end
  end
  @vprint :Fields 1 "Doing Group ($a, $b) with bound $absolute_bound\n"
  @vprint :FieldsNonFancy 1 "Doing Group ($a, $b) with bound $absolute_bound\n"
  L = GAP.Globals.DerivedSeries(G)
  lvl = _real_level(L)
  #First step by hand
  l = GAP.gap_to_julia(Vector{Int64}, GAP.Globals.AbelianInvariants(G))
  @vprint :Fields 1 "contructing abelian extensions with invariants $l \n"
  @vprint :FieldsNonFancy 1 "contructing abelian extensions with invariants $l \n"
  #@vtime :Fields 2
  list = abelian_extensionsQQ(l, root(absolute_bound, div(a,prod(l))), (lvl > 1 || only_real))
  if isempty(list)
    return list
  end
  @vprint :Fields 1 "First step completed\n \n"
  @vprint :FieldsNonFancy 1 "First step completed\n \n"
  for i = 2:length(L)-1

    E1 = GAP.Globals.FactorGroup(L[1], L[i+1])
    l = GAP.gap_to_julia(Vector{Int64}, GAP.Globals.AbelianInvariants(L[i]))
    @vprint :Fields 1 "contructing abelian extensions with invariants $l \n"
    @vprint :FieldsNonFancy 1 "contructing abelian extensions with invariants $l \n"
    o = divexact(a, GAP.Globals.Size(E1))
    bound = root(absolute_bound, o)
    IsoE1 = GAP.Globals.IdGroup(E1)
    @vprint :Fields 1 "Number of fields at the $i -th step: $(length(list)) \n"
    @vprint :FieldsNonFancy 1 "Number of fields at the $i -th step: $(length(list)) \n"
    lG = snf(DiagonalGroup(l))[1]
    invariants = map(Int, lG.snf)
    onlyreal = (lvl > i || only_real)
    @vprint :Fields 1 "Computing obstructions\n"
    @vprint :FieldsNonFancy 1 "Computing obstructions\n"
    #@vtime :Fields 1
    list = check_Brauer_obstruction(list, L, i, invariants)
    @vprint :Fields 1 "Fields to check: $(length(list))\n\n"
    @vprint :FieldsNonFancy 1 "Fields to check: $(length(list))\n\n"
    if isempty(list)
      return FieldsTower[]
    end
    simplify_i_th = true
    if i == length(L)-1 && simplify == false
      simplify_i_th = false
    end
    list = field_extensions(list, bound, IsoE1, invariants, onlyreal, simplify)
    @vprint :Fields 1 "Step $i completed\n"
    @vprint :FieldsNonFancy 1 "Step $i completed\n"
    if isempty(list)
      return FieldsTower[]
    end
  end
  return list

end

function fields(a::Int, b::Int, list::Vector{FieldsTower}, absolute_bound::fmpz; only_real::Bool = false, simplify::Bool = true)
  G = GAP.Globals.SmallGroup(a, b)
  return fields(list, G, absolute_bound, only_real = only_real, simplify = simplify)
end

function fields(list::Vector{FieldsTower}, G, absolute_bound::fmpz; only_real::Bool = false, simplify::Bool = true)
  L = GAP.Globals.DerivedSeries(G)
  lvl = _real_level(L)
  first = true
  for i = 2:length(L)-1
    G1 = GAP.Globals.FactorGroup(L[1], L[i])
    if first && GAP.Globals.Size(G1) != degree(list[1].field)
      continue
    end
    first = false
    E1 = GAP.Globals.FactorGroup(L[1], L[i+1])
    H1 = GAP.Globals.FactorGroup(L[i], L[i+1])
    l = GAP.gap_to_julia(Vector{Int64}, GAP.Globals.AbelianInvariants(H1))
    @vprint :Fields 1 "contructing abelian extensions with invariants $l \n"
    @vprint :FieldsNonFancy 1 "contructing abelian extensions with invariants $l \n"
    o = divexact(GAP.Globals.Size(G), GAP.Globals.Size(E1))
    bound = root(absolute_bound, o)
    IsoE1 = GAP.Globals.IdGroup(E1)
    @vprint :Fields 1 "Number of fields at the $i -th step: $(length(list)) \n"
    @vprint :FieldsNonFancy 1 "Number of fields at the $i -th step: $(length(list)) \n"
    lG = snf(DiagonalGroup(l))[1]
    invariants = map(Int, lG.snf)
    onlyreal = (lvl > i || only_real)
    #First, I search for obstruction.
    if iscyclic(lG)
      @vprint :Fields 1 "Computing obstructions\n"
      @vprint :FieldsNonFancy 1 "Computing obstructions\n"
      #@vtime :Fields 1
      list = check_Brauer_obstruction(list, L, i, invariants[1])
      @vprint :Fields 1 "Fields to check: $(length(list))\n\n"
      @vprint :FieldsNonFancy 1 "Fields to check: $(length(list))\n\n"
    end
    if isempty(list)
      return FieldsTower[]
    end
    list = field_extensions(list, bound, IsoE1, invariants, onlyreal, simplify)
    @vprint :Fields 1 "Step $i completed\n"
    @vprint :FieldsNonFancy 1 "Step $i completed\n"
    if isempty(list)
      return FieldsTower[]
    end
  end
  return list
end


text:
module ContractedMachines

#-- Required modules
# Wiring diagrams
using AlgebraicDynamics
using AlgebraicDynamics.DWDDynam
using Catlab.WiringDiagrams
import AlgebraicDynamics: oapply

# ode solver
using DifferentialEquations
import DifferentialEquations: ODEProblem

# Contract display
using IntervalSets
using PrettyTables
using Printf

# Static contracts
include("StaticContracts.jl")
using .StaticContracts	# Must have period as submodule

#-- Information accesible by using module
export ContractedMachine, oapply, check_contract, ODEProblem

#-- Datatypes used by contract machine
const ContractOutputTable = Dict{KEY, NamedTuple{ (:input, :output),
                                                  Tuple{ Vector{ Pair{Symbol, Bool} },
                                                         Vector{ Pair{Symbol, Bool} } } }
                                 } where KEY <: Union{Symbol, Any}

const ContractTimeTable = Dict{Any, Any}

const ContractOutputBox = NamedTuple{ (:input, :output),
                                      Tuple{ Vector{Bool}, Vector{Bool} } }

#-- Type to display pretty table
struct ContractTable
    table::Union{ContractOutputTable, ContractTimeTable, ContractOutputBox}
end

#-- Type inheritance from AlgebraicDynamics
struct AbstractContractMachine{T, I, S} <: AbstractMachine{T, I, S}
    contract::S
    machine::I
    fcontract::Function
end

#-- Contracts are defined via intervals
const ContractedMachine{T} = AbstractContractMachine{ T, AbstractMachine{T}, StaticContract{T} }

function ContractedMachine{T}(static_contract::StaticContract{T}, machine::AbstractMachine{T};
                            fcontract = nothing) where T
    # contract evaluation function
    if fcontract == nothing
        # Only define function in none is provided
        function fbox(u::AbstractVector, xs::AbstractVector, p=nothing, t=0)
            # evaluate functional inputs
            x = map(x -> typeof(x) <: Function ? x(t) : x, xs)

            # check whether contracts at input ports are satisfied
            Rin = map( (xin,cont) -> xin in cont, x, static_contract.input )

            # check whether contracts at output ports are satisfied
            Rout = map( (rout,cont) -> rout in cont, DWDDynam.readout(machine, u, p, t), static_contract.output )
            return ContractTable( (input = Rin, output = Rout) )
        end

    # make a new machine satisfying these restrictions
        return ContractedMachine{T}(static_contract, machine, fbox)
    else
        return ContractedMachine{T}(static_contract, machine, fcontract)
    end
end

function ContractedMachine{T}(cinput::Vector, nstates::Int, coutput::Vector,
                            dynamics::Function, readout::Function;
                            fcontract = nothing, mtype = :continuous) where T
    # contracts
    static_contract = StaticContract{T}(cinput, coutput)

    # each element in the vectors is a contract for the port on a box
    ninputs = length(cinput)
    noutputs = length(coutput)

    # select a machine type
    if mtype == :continuous
        machine = ContinuousMachine{T}(ninputs, nstates, noutputs, dynamics, readout)
    elseif mtype :discrete
        machine = DiscreteMachine{T}(ninputs, nstates, noutputs, dynamics, readout)
    end
    # make a machine using inner constructor
    return ContractedMachine{T}(static_contract, machine, fcontract=fcontract )
end

#-- Getter functions --
function eval_dynamics(f::ContractedMachine{T}, u, xs, p=nothing, t=0) where T<:Real
    return DWDDynam.eval_dynamics(f.machine, u, xs, p, t)
end

function readout(f::ContractedMachine, u, p=nothing, t=0)
    return DWDDynam.readout(f.machine, u, p, t)
end

#-- Compose multiple contract machines --

function oapply(d::WiringDiagram, ms::Vector{ContractedMachine{T}}) where T<:Real
    # compose contracts
    static_contract = StaticContracts.oapply(d, map(m -> m.contract, ms))

    # Get the initial index of the state vector of each each box.
    # The composition concatenates these vectors into a single column.
    nstate = 0
    index = Array{ UnitRange{Int} }(undef, nboxes(d))

    for i in 1:nboxes(d)
        nstate += nstates(ms[i].machine)                           # initial index of vector
        index[i] = nstate - nstates(ms[i].machine) + 1 : nstate    # store vector indeces
    end

    # get the name of each box
    name = map( box -> box.value, boxes(d) )

    #---- Evaluate contracts
    function fcontract(u::AbstractVector, x::AbstractVector, p=nothing, t=0)
        # get the output of each box
        rout = map( id -> ContractedMachines.readout(ms[id], u[index[id]], p, t ), 1:nboxes(d) )

        # evaluate the contract function
        fout = Array{Dict}(undef, nboxes(d))

        for id in 1:nboxes(d)                       # check all boxes in the diagram
            # collect the inputs of each box
            xin = zeros( ninputs(ms[id].machine) )

            for w in in_wires(d, id)                # check all wires going into a box
                if w.source.box != input_id(d)      # internal inputs use are the readout of another box
                    xin[w.target.port] += rout[w.source.box][w.source.port]
                else                                # external inputs make use of a given vector
                    x_source = x[w.source.port]

                    # check for function input
                    if typeof(x_source) <: Function
                        xin[w.target.port] += x_source(t)
                    else
                        xin[w.target.port] += x_source
                    end
                end
            end

            # evaluate the contract function of each box
            foutput = ms[id].fcontract( u[index[id]], xin, p, t ).table

            # for atomic boxes, assign the box name to the inputs and outputs
            if typeof(foutput) <: NamedTuple
                fout[id] = Dict(name[id] => (input = input_ports(d, id) .=> foutput.input,
                                             output = output_ports(d, id) .=> foutput.output))
            else    # for nested boxes, append the name of the parent box to the child boxes
                fout[id] = Dict( (name[id] .=> keys(foutput)) .=> values(foutput) )
            end
        end
        # return a contract table with a box directory assigned to inputs and outputs
        return ContractTable( merge(fout...) )
    end

    # compose the dynamics functions
    machine = DWDDynam.oapply(d, map(m -> m.machine, ms))

    # return the composed machine
    return ContractedMachine{T}(static_contract, machine, fcontract=fcontract)
end

# Identify during which time intervals a signal is zero [false == contract is violated]
function failureInterval(arr::AbstractVector, time=nothing, output="time")
    # check each element in the array
    index = Array{Tuple}(undef, 0)
    start = 1

    for i in 1:length(arr)             # no previous element
        if i != 1
            if arr[i] > arr[i-1]       # rising signal
                push!( index, (start, i-1) )

            elseif arr[i] < arr[i-1]   # falling signal
                start = i
            end

            if i == length(arr) && arr[end] <= 0   # end of interval
                push!( index, (start, i) )
            end
        end
    end

    # map indeces to given array
    if time != nothing && output == "time"
        index = map( set -> (time[set[1]], time[set[2]]), index)
    end
    # output 2d array of start and stop times
    return index
end

# Assign contract failure times to each wire of each box
function check_contract(sol::T1, machine::ContractedMachine{T2},
                        x0::AbstractVector, p=nothing, t=0;
                        output="time" ) where {T1<:ODESolution, T2<:Real}

    # evalutate the contract function throughout time interval
    fout = map( t -> machine.fcontract(sol(t), x0, p, t), sol.t )

    # Store time at which each contract fails for each box directory:
    # Initial evaluation of contracts. Used to get the names and wires of each box.
    set0 = fout[begin].table
    dict = Dict()

    # go through all directories
    for key in keys(set0)
        # function to map indeces of contract outputs
        mapOut = index -> map(i -> set0[key][index][i].first => begin
                                  # array of the contract state for each time step
                                  contract = map(set -> set.table[key][index][i].second, fout);
                                  # find the duration of 0's in the array
                                  failureInterval(contract, sol.t, output);
                              end,
                              1:length(set0[key][index]) )
        # store the times at the directory
        dict[key] = ( input=mapOut(:input), output=mapOut(:output) )
    end

    # output a data structure like that of fcontract but with failure intervals.
    return ContractTable(dict)
end

#-- Pretty printing of output

# display contract machines as product of intervals
function Base.show(io::IO, vf::ContractedMachine)
    display(vf.contract)
    print("[Machine]") 		# Indicate it's a machine
end

# display contract tables (output of fcontract and eval_contract) as a pretty table
function Base.show(io::IO, data::ContractTable)
    # get the data structure
    dict = data.table

    # select the appriate table
    if typeof(dict) <: ContractOutputBox        # atomic box
        # assign a contract state to each port
        mapOut = index -> join( string.(1:length(dict[index])) .* " : " .* string.(dict[index]), "\n" )

        # do not assign a box as there is no diagram
        output = hcat( [mapOut(:input)], [mapOut(:output)] )
        header = ( ["input", "output"], ["port: contract", "port: contract"] )

    else # composed boxes
        box = string.(keys(dict))

        if typeof(dict) <: ContractOutputTable  # output of fcontract
            mapOut = index -> map( x -> join(
                map(pair -> string(pair.first) * " : " * string(pair.second), x[index]), "\n"),
                                   values(dict) )

            header =  ( ["box", "input", "output"], ["directory", "wire: contract", "wire: contract"] )

        else    # output of eval_contract
            mapOut = index::Symbol -> map( value -> begin
                                           # make array with strings of time intervals for each box
                                           # Note: need way to change formatting for intergers
                                           out = map( pair ->
                                                      map( set -> @sprintf("%s : %f , %f", pair.first, set...), pair.second ),
                                                 value[index] );

                                           # combine strings into large string seperated by line jumps
                                           join( vcat(out...), "\n")
                                           end,
                                           values(dict) )

            header = ( ["box", "input", "output"], ["directory", "wire: failure interval", "wire: failure interval"] )
        end
        # both types have the same structure
        output = hcat(box, mapOut(:input), mapOut(:output))
    end
    # display a pretty table
    pretty_table(output, header = header, backend=:html, linebreaks=true, tf=tf_html_minimalist, alignment=:c)
end

#-- helper functions

# compose machines given the name of each box
function oapply(d::WiringDiagram, ms::Dict{Symbol, ContractedMachine{T}}) where T<:Real
    oapply(d, map(box -> ms[box.value], boxes(d)) )
end

# solve the dynamics of the contract machines
function ODEProblem(m::ContractedMachine{T}, u0::AbstractVector, xs::AbstractVector, tspan::Tuple, p=nothing) where T<:Real
    ODEProblem(m.machine, u0, xs, tspan, p)
end

end # module


text:
# Model for dual-clocked fifo.
import Base: write, read

abstract type AbstractFifo end

##################
# Fifo Interface #
##################

"""
    isreadready(f::AbstractFifo)

Return `true` if `read(f)` may be called and will return valid read data.
"""
isreadready(f::AbstractFifo) = false

"""
    peek(f::AbstractFifo)

Return the top element of the read side of the fifo. Must not change that read
element. That is, `peek()` must be save to call multiple times in a row.
"""
peek(f::AbstractFifo) = nothing

"""
    read(f::AbstractFifo)

Return the top element of the read side of the fifo. Mark the fifo so it will
be updated to the next element on the next call to `readupdate!(f)`.
"""
read(f::AbstractFifo) = nothing

"""
    readupdate!(f::AbstractFifo)

Update the read side of `f`. Usually analogous to "clocking" the read side of
the fifo.
"""
readupdate!(f::AbstractFifo) = nothing

"""
    iswriteready(f::AbstractFifo)

Return `true` if `write(f)` followed by `writeupdate!` will be successful.
"""
iswriteready(f::AbstractFifo) = false

"""
    write(f::AbstractFifo, data)

Send data to `f` to store on the next call to `writeupdate!(f)`.
"""
write(f::AbstractFifo, data) = nothing

"""
    writeupdate!(f::AbstractFifo)

Update the write side of `f`. Usually analogous to "clocking" the write side
of the fifo.
"""
writeupdate!(f::AbstractFifo, data) = nothing

# A test fifo for testing various methods involving cores.
@with_kw mutable struct EmptyFullFifo <: AbstractFifo
    empty :: Bool = false
    full  :: Bool = false
end

isreadready(t::EmptyFullFifo) = !t.empty
iswriteready(t::EmptyFullFifo) = !t.full


text:
# Evolvability using a dictionary to keep track of the phenotypes that contribute to evolution evolvability.
# Additional objectives are to simplify the logical structure form geno_complexity() and to
#   use neutral_evolution() and lambda_evolution() instead of mut_evolve()
export run_evo_dict, print_matrix

# ncircuits is the number of circuits (chromosomes) per goal used to compute evolvability
# ntries is the number of evolution attempts done in trying to evolve each chromosome
# maxsteps is the maximum number of steps for each evolution
# mutrate is the mutation rate if lambda_evolution is used.
# mutrate<=0 means use neutral_evolution
# As of 7/3/21 does not work for cartesian==false because mutata_all() is not definined in LinCircuit.jl, and
# it looks like the definition of circuit_int() in LinCircuit.jl overrides the definition in Chromosome.jl
#   (so I turned off the export in LinCircuit.jl).
function run_evo_dict( p::Parameters, gl::GoalList, ncircuits::Int64, numtries::Int64, maxsteps::Int64, mutrate::Float64=0.0;
      cartesian::Bool=true, csvfile::String="" )
  df = DataFrame()
  #df.goal = Goal[]
  df.goal = MyInt[]
  df.numinputs = Int64[]
  df.numoutputs = Int64[]
  df.numgates = Int64[]
  df.numlevsback = Int64[]
  df.ncircuits = Int64[]
  df.numtries = Int64[]
  df.maxsteps = Int64[]
  df.mutrate = Float64[]
  df.successes = Int64[]
  df.tries = Int64[]
  df.totalsteps = Int64[]
  df.robust_sum = Int64[]
  df.count_circs = Int64[]
  df.evo_count = Int64[]
  #println("size(df): ",size(df))
  results = pmap( g->evo_dict( p, g, ncircuits, numtries, maxsteps, mutrate, cartesian=cartesian ), gl )
  for res in results
    result_list = res[2:end]  # All results from evo_dict() except goal. Last element is evdict
    #println("result_list: ",result_list)
    param_list = [ res[1][1], p.numinputs, p.numoutputs, p.numinteriors, p.numlevelsback, ncircuits, numtries, maxsteps, mutrate ]
    df_list = vcat( result_list[1:(end-1)], length(result_list[end]))  # length(result_list[end]) is evo_count
    #println("df_list: ",df_list)
    dfrow = vcat(param_list,df_list)
    #println("length(dfrow): ",length(dfrow))
    #println("dfrow: ",dfrow)
    push!(df,vcat(param_list,df_list))
  end
  evdict = Dict{MyInt, Int64 }[]
  evset = Set{MyInt}[]
  evlist = Vector{Int64}[]
  pheno_matrix = zeros(Int64,length(gl),length(gl))
  gli = map( x->x[1], gl )
  println("gli: ",gli)
  for res in results
    push!(evdict,res[end])
    push!(evset,Set(keys(res[end])))
    evo_list = sort([Int64(k) for k in keys(res[end])])
    #println("evo_list: ",evo_list)
    push!(evlist,evo_list)
    goal_index = findfirst( x->x==res[1], gl )
    update_pheno_matrix( res[1], gli, res[end], pheno_matrix )
  end
  #println("evlist: ",evlist)
  insertcols!(df,size(df)[2]+1,:evo_phenos=>evlist)
  intsdf = DataFrame()
  for g in gl
    insertcols!(intsdf,size(intsdf)[2]+1,Symbol(@sprintf("0x00%x",g[1]))=>Int64[])
  end
  for evs1 in evset
    intersects = [ length(intersect(evs1,evs2)) for evs2 in evset ]
    println("intersects: ",intersects)
    push!(intsdf,intersects)
  end
  print_matrix(pheno_matrix)
  phdf_hex = matrix_to_dataframe( pheno_matrix, gli, hex=true )
  phdf_dec = matrix_to_dataframe( pheno_matrix, gli, hex=false )
  if length(csvfile) > 0
    hostname = chomp(open("/etc/hostname") do f read(f,String) end)
    open( csvfile[1:end-4]*"_phmatrix_hex.csv", "w" ) do f
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      print_parameters(f,p,comment=true)
      println(f,"# ncircuits: ",ncircuits)
      CSV.write( f, phdf_hex, append=true, writeheader=true )
    end
    open( csvfile[1:end-4]*"_phmatrix_dec.csv", "w" ) do f
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      print_parameters(f,p,comment=true)
      println(f,"# ncircuits: ",ncircuits)
      CSV.write( f, phdf_dec, append=true, writeheader=true )
    end
    open( csvfile[1:end-4]*"_ints.csv", "w" ) do f
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      println(f,"# numinputs: ",p.numinputs)
      println(f,"# numoutputs: ",p.numoutputs)
      println(f,"# numinteriors: ",p.numinteriors)
      println(f,"# numlevelsback: ",p.numlevelsback)
      println(f,"# ncircuits: ",ncircuits)
      println(f,"# numtries: ",numtries)
      println(f,"# maxsteps: ",maxsteps)
      CSV.write(f, intsdf, append=true, writeheader=true )
    end
    open( csvfile, "w" ) do f
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      println(f,"# ncircuits: ",ncircuits)
      println(f,"# numtries: ",numtries)
      println(f,"# maxsteps: ",maxsteps)
      CSV.write(f, df, append=true, writeheader=true )
    end
  end
  (intsdf,phdf_hex,df)
end

# Same as run_evo_dict, only for just one goal.
# if mutrate <= 0.0 use neutral_evolution, otherwise use lambda_evolution
function evo_dict( p::Parameters, g::Goal, ncircuits::Int64, numtries::Int64, maxsteps::Int64, mutrate::Float64;
    cartesian::Bool=true )
  funcs = cartesian ? default_funcs(p.numinputs) : lin_funcs(p.numinputs)
  # evdict counts the number of occurences of the phenotypes that are mutationally adjacent to g
  evdict = Dict{MyInt, Int64 }()
  # circuit_dict is indexed by the integer representation of the circuit.
  circuit_dict = cartesian ? Dict{Int64,Int64}() : Dict{Vector{Int64},Int64}()
  robust_sum = 0
  total_steps = 0
  successes = 0
  t = 0
  while t < numtries && successes < ncircuits
    t += 1
    if cartesian
      c = random_chromosome( p, funcs )
    else
      c = rand_lcircuit(p,funcs)
    end
    if mutrate <= 0.0
      (new_c,steps) = neutral_evolution( c, g, maxsteps )
    else
      (new_c,steps) = lambda_evolution( c, g, maxsteps, mutrate )
    end
    total_steps += steps
    if steps < maxsteps
      if haskey(circuit_dict,circuit_int(new_c))
        circuit_dict[circuit_int(new_c)] += 1
      else
        circuit_dict[circuit_int(new_c)] = 1
      end
      successes += 1
      ov = output_values( new_c )
      @assert g == ov
      #println("g: ",g,"  ov: ",ov)
      neighboring_phenos = mutate_all( new_c, funcs, output_outputs=:true )
      #println("neighboring_phenos: ",neighboring_phenos[1:5])
      robust_outputs = filter( x->x==ov, neighboring_phenos )
      robust_sum += length(robust_outputs)
      for ph in neighboring_phenos
        if haskey(evdict,ph[1])
          evdict[ph[1]] += 1
        else
          evdict[ph[1]] = 1
        end
      end
    end  # if
  end  # while
  count_circs = length(keys(circuit_dict))
  return [ g, successes, t, total_steps, robust_sum, count_circs, evdict ]
end

function update_pheno_matrix( g::Goal, gli::Vector{MyInt}, evdict::Dict{MyInt,Int64}, pheno_matrix::Array{Int64,2} )
  #println("update_pheno_matrix: goal: ",g,"  gli: ",gli, "  evdict: ",evdict)
  println("update_pheno_matrix: goal: ",g,"  gli: ",gli )
  goal_index = findfirst( x->x==g[1], gli )
  for k in keys(evdict)
    ph_index = findfirst( x->x==k, gli )
    if ph_index != nothing
      #println("k: ",k,"  goal_index: ",goal_index,"  ph_index: ",ph_index)
      pheno_matrix[goal_index,ph_index] = evdict[k]
    end
  end
  #print_matrix(pheno_matrix)
end

function print_matrix( f::IO, matrix::Array{Int64,2} )
  (rows,columns) = size(matrix)
  for r = 1:rows
    for c = 1:columns
      @printf(f,"%4d ",matrix[r,c])
    end
    println(f)
  end
end

function print_matrix( matrix::Array{Int64,2} )
  print_matrix( Base.stdout, matrix )
end

# Convert matrix to a dataframe with a goal column and with column names taken from gli.
function matrix_to_dataframe( pheno_matrix::Array{Int64,2}, gli::Vector{MyInt}; hex::Bool=true)
  (rows,columns) = size(pheno_matrix)
  @assert rows==length(gli)
  df = DataFrame()
  if hex
    rnames = map(x->@sprintf("0x%x",x), gli )
  else
    rnames = map(x->@sprintf("%d",x), gli )
  end
  df.goal = rnames
  for c = 1:columns
    insertcols!(df,size(df)[2]+1,Symbol(rnames[c])=>pheno_matrix[:,c])
  end
  df
end


text:
"""
    KMarkovEnvironment{OV, M<:POMDP, S, R<:AbstractRNG} <: AbstractEnvironment{OV}
A k-markov wrapper for MDPs and POMDPs, given a MDP or POMDP create an AbstractEnvironment where s_t = (o_t, ..., o_t-k)
The K-Markov observation is represented by a vector of k observations.
"""
mutable struct KMarkovEnvironment{OV, M<:POMDP, S, R<:AbstractRNG} <: AbstractEnvironment{OV}
    problem::M
    k::Int64
    state::S
    obs::Vector{OV}
    rng::R
end
function KMarkovEnvironment(problem::M,
                            ov::Type{<:AbstractArray} = obsvector_type(problem);
                            k::Int64=1,
                            rng::AbstractRNG=MersenneTwister(0)
                            ) where {M<:POMDP}
    # determine size of obs vector
    s = rand(rng, initialstate(problem))
    o = rand(rng, initialobs(problem, s))
    obs = convert_o(ov, o, problem)
    # init vector of obs
    obsvec = fill(zeros(eltype(ov), size(obs)...), k)
    return KMarkovEnvironment{ov, M, typeof(s), typeof(rng)}(
        problem,
        k,
        rand(rng, initialstate(problem)),
        obsvec,
        rng)
end

"""
    reset!(env::KMarkovEnvironment{OV})
Reset an POMDP environment by sampling an initial state,
generating an observation and returning it.
"""
function reset!(env::KMarkovEnvironment{OV}) where OV
    s = rand(env.rng, initialstate(env.problem))
    env.state = s
    o = rand(env.rng, initialobs(env.problem, s))
    obs = convert_o(OV, o, env.problem)
    fill!(env.obs, obs)
    return env.obs
end

"""
    step!(env::POMDPEnvironment{OV}, a::A)
Take in an POMDP environment, and an action to execute, and
step the environment forward. Return the observation, reward,
terminal flag and info
"""
function step!(env::KMarkovEnvironment{OV}, a::A) where {OV, A}
    s, o, r, info = @gen(:sp, :o, :r, :info)(env.problem, env.state, a, env.rng)
    env.state = s
    t = isterminal(env.problem, s)
    obs = convert_o(OV, o, env.problem)
    # shift the old observation to lower indices
    for i=1:env.k-1
        env.obs[i] = env.obs[i + 1]
    end
    env.obs[env.k] = obs
    return env.obs, r, t, info
end

"""
    actions(env::KMarkovEnvironment)
Return an action object that can be sampled with rand.
"""
function POMDPs.actions(env::KMarkovEnvironment)
    return actions(env.problem)
end

"""
    sample_action(env::Union{POMDPEnvironment, MDPEnvironment})
Sample an action from the action space of the environment.
"""
function sample_action(env::KMarkovEnvironment)
    return rand(env.rng, actions(env))
end

"""
    obs_dimensions(env::KMarkovEnvironment{OV})
returns the size of the observation vector.
The object return by `step!` and `reset!` is a vector of k observation vector of size `obs_dimensions(env)`
It generates an initial state, converts it to an array and returns its size.
"""
function obs_dimensions(env::KMarkovEnvironment{OV}) where OV
    s = rand(env.rng, initialstate(env.problem))
    o = rand(env.rng, initialobs(env.problem, s))
    obs_dim = size(convert_o(OV, o, env.problem))
    return (obs_dim..., env.k)
end


text:
function say_name(name)
  println("$name")
end

say_name("mei")
# ワンライナーでかく
say_name_oneline(name) = println("$name")
say_name_oneline("mei_one")

# 複数行
say_name_lines(name) = begin
    output = name * ": 複数行で書くよ"
    println(output)
end

say_name_lines("mei")

# lambda式

say_name_lambda = name -> println("$name da")
say_name_lambda("mae")

# juliaでは最終で返す値がかえる return はなくてもいい

function sample(x,y)
    x + y
end

println(sample(3,2))

# operatorを利用したやり方

println(+(1, 2, 3, 4))
println(*(1, 2, 3, 4))


text:
using LaserFileSplitter

svgfile = "data/exp_raw/DeathMoonDevil.svg"
outfolder = "data/processed/chunked_deathmoondevil3"
max_num_chunks = 10
save_svg_in_chunks(svgfile, outfolder, max_num_chunks)

max_stars_repo_name: zoemcc/LaserFileSplitter

max_issues_repo_name: zoemcc/LaserFileSplitter

max_forks_repo_name: zoemcc/LaserFileSplitter


text:
module Bench

include("args.jl")
#
import Bench.Args: generate, smallest, tag

import BenchmarkTools: BenchmarkGroup, @benchmarkable

function valid(f, args)
    try
        f(args...)
        return true
    catch
        return false
    end
end

function push!(suite, ns, f, args)
    key = Tuple( (f, tag.(args)...) )
    suite[key] = BenchmarkGroup(collect(key))
    for n in ns
        arg = generate(n, args)
        suite[key][n] = @benchmarkable $f( $arg... )
    end
end

function productsuite(ns, fs, arglist...)
    suite = BenchmarkGroup()
    for f in fs
        for args in Iterators.product(arglist...)
            push!(suite, ns, f, args)
        end
    end
    return suite
end

end # module


text:
include("../highs_common.jl")

script = build_script(shared_libs = "ON")

products = [
    LibraryProduct("libhighs", :libhighs),
    ExecutableProduct("highs", :highs),
]

platforms = expand_cxxstring_abis(platforms)

dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
]

build_tarballs(
    ARGS,
    name,
    version,
    sources,
    script,
    platforms,
    products,
    dependencies;
    preferred_gcc_version = v"6",
    julia_compat = "1.6",
)

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
#=
Random Forests in the "Predict, then Optimize" (PO) setting. Uses package DecisionTree.
=#

using Parameters, DecisionTree

## Type for Random Forest algorithm parameters with default values
@with_kw struct rf_parms
    num_trees::Int = 100
    num_features_per_split::Int = -1
end

"""
    train_random_forests_po(X::Matrix{Float64}, c::Matrix{Float64};
        num_trees = 100, num_features_per_split = -1)
Train d different Random Forests regression models, each to predict a different
component of the cost vector as a function of the features. Returns a list of the
d different RF models each trained from the `DecisionTree` package.
# Arguments
- `X`: p x n training set feature matrix
- `c`: d x n training set matrix of cost vectors
- `num_trees` is the total number of trees in each RF model
- `num_features_per_split` is the number of randomly sampled features at each split
"""
function train_random_forests_po(X::Matrix{Float64}, c::Matrix{Float64};
    rf_alg_parms::rf_parms = rf_parms())

    @unpack num_trees, num_features_per_split = rf_alg_parms

    (p, n) = size(X)
    (d, n2) = size(c)
    if n != n2
        error("Dimensions of the input are mismatched.")
    end

    X_t = copy(transpose(X))

    # If num_features_per_split is not specified, use default for regression
    if num_features_per_split < 1
        num_features_per_split = ceil(Int, p/3)
    end

    rf_model_list = DecisionTree.Ensemble[]

    # train one model for each component
    for j = 1:d
        c_vec = c[j, :]
        rf_model = build_forest(c_vec, X_t, num_features_per_split, num_trees)
        push!(rf_model_list, rf_model)
    end

    return rf_model_list
end

"""
    predict_random_forests_po(rf_model_list, X_new)
Given a list of random forests models trained previously by the `run_random_forests_po`
function and a p x n matrix of new features values, return a d x n matrix of predictions.
# Arguments
- `rf_model_list` is the length d list of RF model objects, each to predict a different component of the cost vector based on a feature vector of length p.
- `X_new`: p x n feature matrix. Each column is one of n observations.
"""
function predict_random_forests_po(rf_model_list, X_new)
    (p, n) = size(X_new)
    d = length(rf_model_list)

    X_new_t = copy(transpose(X_new))

    preds = zeros(d, n)
    for j = 1:d
        preds[j, :] = apply_forest(rf_model_list[j], X_new_t)
    end

    return preds
end

max_stars_repo_name: katerinagiannoutsou/dso-coding-project

max_issues_repo_name: katerinagiannoutsou/dso-coding-project

max_forks_repo_name: katerinagiannoutsou/dso-coding-project


text:
################################################################################

"Reset hidden Markov model object."
function reset!(self::HMM)
  self.model = [χ .= [0; repeat(-1:-1, length(self.model[1]) - 1)] for χ = self.model]
end

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