text:
abstract type MultiFieldStyle end

struct ConsecutiveMultiFieldStyle <: MultiFieldStyle end

struct StridedMultiFieldStyle <: MultiFieldStyle end

"""
    struct MultiFieldFESpace{S<:MultiFieldStyle,B} <: FESpace
      spaces::Vector{<:SingleFieldFESpace}
      multi_field_style::S
      constraint_style::Val{B}
    end
"""
struct MultiFieldFESpace{MS<:MultiFieldStyle,CS<:ConstraintStyle,V} <: FESpace
  vector_type::Type{V}
  spaces::Vector{<:SingleFieldFESpace}
  multi_field_style::MS
  constraint_style::CS
  function MultiFieldFESpace(
    ::Type{V},
    spaces::Vector{<:SingleFieldFESpace},
    multi_field_style::MultiFieldStyle) where V
    @assert length(spaces) > 0

    MS = typeof(multi_field_style)
    if any( map(has_constraints,spaces) )
      constraint_style = Constrained()
    else
      constraint_style = UnConstrained()
    end
    CS = typeof(constraint_style)
    new{MS,CS,V}(V,spaces,multi_field_style,constraint_style)
  end
end

"""
    MultiFieldFESpace(spaces::Vector{<:SingleFieldFESpace})
"""
function MultiFieldFESpace(spaces::Vector{<:SingleFieldFESpace})
  Ts = map(get_dof_value_type,spaces)
  T = typeof(*(map(zero,Ts)...))
  MultiFieldFESpace(Vector{T},spaces,ConsecutiveMultiFieldStyle())
end

function MultiFieldFESpace(::Type{V},spaces::Vector{<:SingleFieldFESpace}) where V
  MultiFieldFESpace(V,spaces,ConsecutiveMultiFieldStyle())
end

MultiFieldStyle(::Type{MultiFieldFESpace{S,B,V}}) where {S,B,V} = S()
MultiFieldStyle(f::MultiFieldFESpace) = MultiFieldStyle(typeof(f))

# Implementation of FESpace

function FESpaces.get_triangulation(f::MultiFieldFESpace)
  s1 = first(f.spaces)
  trian = get_triangulation(s1)
  @check all(map(i->trian===get_triangulation(i),f.spaces))
  trian
end

function FESpaces.num_free_dofs(f::MultiFieldFESpace)
  n = 0
  for U in f.spaces
    n += num_free_dofs(U)
  end
  n
end

function FESpaces.get_free_dof_ids(f::MultiFieldFESpace)
  get_free_dof_ids(f,MultiFieldStyle(f))
end

function FESpaces.get_free_dof_ids(f::MultiFieldFESpace,::MultiFieldStyle)
  @abstractmethod
end

function FESpaces.get_free_dof_ids(f::MultiFieldFESpace,::ConsecutiveMultiFieldStyle)
  block_num_dofs = Int[]
  for U in f.spaces
    push!(block_num_dofs,num_free_dofs(U))
  end
  BlockArrays.blockedrange(block_num_dofs)
end

FESpaces.get_dof_value_type(f::MultiFieldFESpace{MS,CS,V}) where {MS,CS,V} = eltype(V)

FESpaces.get_vector_type(f::MultiFieldFESpace) = f.vector_type

FESpaces.ConstraintStyle(::Type{MultiFieldFESpace{S,B,V}}) where {S,B,V} = B()

struct MultiFieldFEBasisComponent{B} <: FEBasis
  cell_basis::AbstractArray
  single_field::B
  fieldid::Int
  nfields::Int
  function MultiFieldFEBasisComponent(
    single_field::SingleFieldFEBasis,fieldid::Integer,nfields::Integer)
    function block_dofs(cell_bs,::TestBasis,fieldid,nfields)
      cell_basis = lazy_map(BlockMap(nfields,fieldid),cell_bs)
    end
    function block_dofs(cell_bs,::TrialBasis,fieldid,nfields)
      cell_basis = lazy_map(BlockMap((1,nfields),fieldid),cell_bs)
    end
    B = typeof(single_field)
    cell_bs = get_data(single_field)
    bsty = BasisStyle(single_field)
    cell_basis = block_dofs(cell_bs,bsty,fieldid,nfields)
    new{B}(cell_basis,single_field,fieldid,nfields)
  end
end

CellData.get_data(f::MultiFieldFEBasisComponent) = f.cell_basis
CellData.get_triangulation(f::MultiFieldFEBasisComponent) = get_triangulation(f.single_field)
FESpaces.BasisStyle(::Type{<:MultiFieldFEBasisComponent{B}}) where B = BasisStyle(B)
CellData.DomainStyle(::Type{<:MultiFieldFEBasisComponent{B}}) where B = DomainStyle(B)
function FESpaces.CellData.similar_cell_field(
  f::MultiFieldFEBasisComponent,cell_data,trian,ds::DomainStyle)
  @notimplemented
end
function FESpaces.similar_fe_basis(
  f::MultiFieldFEBasisComponent,cell_data,trian,bs::BasisStyle,ds::DomainStyle)
  @notimplemented
end

for fun in (:gradient,:DIV,:∇∇)
  @eval begin
    function $fun(f::MultiFieldFEBasisComponent)
      g = $fun(f.single_field)
      MultiFieldFEBasisComponent(g,f.fieldid,f.nfields)
    end
  end
end

function CellData.change_domain(
  a::MultiFieldFEBasisComponent,
  tdomain::DomainStyle)
  sf = change_domain(a.single_field,tdomain)
  MultiFieldFEBasisComponent(sf,a.fieldid,a.nfields)
end

function CellData.change_domain(
  a::MultiFieldFEBasisComponent,
  ttrian::Triangulation,
  tdomain::DomainStyle)
  sf = change_domain(a.single_field,ttrian,tdomain)
  MultiFieldFEBasisComponent(sf,a.fieldid,a.nfields)
end

function MultiFieldFEBasisComponent(
  single_field::CellFieldAt{S,<:SingleFieldFEBasis},
  fieldid::Integer,
  nfields::Integer) where S

  sf = single_field.parent
  mf = MultiFieldFEBasisComponent(sf,fieldid,nfields)
  CellFieldAt{S}(mf)
end

function CellData.change_domain(
  a::CellFieldAt{S,<:MultiFieldFEBasisComponent},
  tdomain::DomainStyle) where S
  mf = a.parent
  sfin = CellFieldAt{S}(mf.single_field)
  sfout = change_domain(sfin,tdomain)
  MultiFieldFEBasisComponent(sfout,mf.fieldid,mf.nfields)
end

function CellData.change_domain(
  a::CellFieldAt{S,<:MultiFieldFEBasisComponent},
  ttrian::Triangulation,
  tdomain::DomainStyle) where S
  mf = a.parent
  sfin = CellFieldAt{S}(mf.single_field)
  sfout = change_domain(sfin,ttrian,tdomain)
  MultiFieldFEBasisComponent(sfout,mf.fieldid,mf.nfields)
end

function FESpaces.get_fe_basis(f::MultiFieldFESpace)
  nfields = length(f.spaces)
  all_febases = MultiFieldFEBasisComponent[]
  for field_i in 1:nfields
    dv_i = get_fe_basis(f.spaces[field_i])
    @assert BasisStyle(dv_i) == TestBasis()
    dv_i_b = MultiFieldFEBasisComponent(dv_i,field_i,nfields)
    push!(all_febases,dv_i_b)
  end
  MultiFieldCellField(all_febases)
end

function FESpaces.get_trial_fe_basis(f::MultiFieldFESpace)
  nfields = length(f.spaces)
  all_febases = MultiFieldFEBasisComponent[]
  for field_i in 1:nfields
    du_i = get_trial_fe_basis(f.spaces[field_i])
    @assert BasisStyle(du_i) == TrialBasis()
    du_i_b = MultiFieldFEBasisComponent(du_i,field_i,nfields)
    push!(all_febases,du_i_b)
  end
  MultiFieldCellField(all_febases)
end

function FESpaces.FEFunction(fe::MultiFieldFESpace, free_values)
  blocks = map(1:length(fe.spaces)) do i
    free_values_i = restrict_to_field(fe,free_values,i)
    FEFunction(fe.spaces[i],free_values_i)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end

function FESpaces.EvaluationFunction(fe::MultiFieldFESpace, free_values)
  blocks = map(1:length(fe.spaces)) do i
    free_values_i = restrict_to_field(fe,free_values,i)
    EvaluationFunction(fe.spaces[i],free_values_i)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end

function CellData.CellField(fe::MultiFieldFESpace,cell_values)
  single_fields = map(1:length(fe.spaces)) do i
    cell_values_field = lazy_map(a->a.array[i],cell_values)
    CellField(fe.spaces[i],cell_values_field)
  end
  MultiFieldCellField(single_fields)
end

"""
    restrict_to_field(f::MultiFieldFESpace,free_values::AbstractVector,field::Integer)
"""
function restrict_to_field(f::MultiFieldFESpace,free_values::AbstractVector,field::Integer)
  _restrict_to_field(f,MultiFieldStyle(f),free_values,field)
end

function  _restrict_to_field(f,::MultiFieldStyle,free_values,field)
  @notimplemented
end

function  _restrict_to_field(f,::ConsecutiveMultiFieldStyle,free_values,field)
  offsets = compute_field_offsets(f)
  U = f.spaces
  pini = offsets[field] + 1
  pend = offsets[field] + num_free_dofs(U[field])
  SubVector(free_values,pini,pend)
end

"""
    compute_field_offsets(f::MultiFieldFESpace)
"""
function compute_field_offsets(f::MultiFieldFESpace)
  @assert MultiFieldStyle(f) == ConsecutiveMultiFieldStyle()
  U = f.spaces
  n = length(U)
  offsets = zeros(Int,n)
  for i in 1:(n-1)
    Ui = U[i]
    offsets[i+1] = offsets[i] + num_free_dofs(Ui)
  end
  offsets
end

function FESpaces.get_cell_isconstrained(f::MultiFieldFESpace)
  msg = """\n
  This method does not make sense for multi-field
  since each field can be defined on a different triangulation.
  Pass a triangulation in the second argument to get
  the constrain flag for the corresponding cells.
  """
  trians = map(get_triangulation,f.spaces)
  trian = first(trians)
  @check all(map(t->t===trian,trians)) msg
  get_cell_isconstrained(f,trian)
end

function FESpaces.get_cell_isconstrained(f::MultiFieldFESpace,trian::Triangulation)
  data = map(f.spaces) do space
    trian_i = get_triangulation(space)
    if is_change_possible(trian_i,trian)
      get_cell_isconstrained(space,trian)
    else
      Fill(false,num_cells(trian))
    end
  end
  lazy_map( (args...) -> +(args...)>0,  data...)
end

#function FESpaces.get_cell_isconstrained(f::MultiFieldFESpace,trian::SkeletonTriangulation)
#  plus = get_cell_isconstrained(f,trian.plus)
#  minus = get_cell_isconstrained(f,trian.minus)
#  lazy_map((l,r)-> l||r,plus,minus)
#end

function FESpaces.get_cell_is_dirichlet(f::MultiFieldFESpace)
  msg = """\n
  This method does not make sense for multi-field
  since each field can be defined on a different triangulation.
  Pass a triangulation in the second argument to get
  the constrain flag for the corresponding cells.
  """
  trians = map(get_triangulation,f.spaces)
  trian = first(trians)
  @check all(map(t->t===trian,trians)) msg
  get_cell_is_dirichlet(f,trian)
end

function FESpaces.get_cell_is_dirichlet(f::MultiFieldFESpace,trian::Triangulation)
  data = map(f.spaces) do space
    trian_i = get_triangulation(space)
    if is_change_possible(trian_i,trian)
      get_cell_is_dirichlet(space,trian)
    else
      Fill(false,num_cells(trian))
    end
  end
  lazy_map( (args...) -> +(args...)>0,  data...)
end

#function FESpaces.get_cell_is_dirichlet(f::MultiFieldFESpace,trian::SkeletonTriangulation)
#  plus = get_cell_is_dirichlet(f,trian.plus)
#  minus = get_cell_is_dirichlet(f,trian.minus)
#  lazy_map((l,r)-> l||r,plus,minus)
#end

function FESpaces.get_cell_constraints(f::MultiFieldFESpace)
  msg = """\n
  This method does not make sense for multi-field
  since each field can be defined on a different triangulation.
  Pass a triangulation in the second argument to get
  the constrains for the corresponding cells.
  """
  trians = map(get_triangulation,f.spaces)
  trian = first(trians)
  @check all(map(t->t===trian,trians)) msg
  get_cell_constraints(f,trian)
end

function FESpaces.get_cell_constraints(f::MultiFieldFESpace,trian::Triangulation)
  nfields = length(f.spaces)
  blockmask = [ is_change_possible(get_triangulation(Vi),trian) for Vi in f.spaces ]
  active_block_ids = findall(blockmask)
  active_block_data = Any[ get_cell_constraints(f.spaces[i],trian) for i in active_block_ids ]
  blockshape = (nfields,nfields)
  blockindices = [(i,i) for i in active_block_ids]
  lazy_map(BlockMap(blockshape,blockindices),active_block_data...)
end

#function FESpaces.get_cell_constraints(f::MultiFieldFESpace,trian::SkeletonTriangulation)
#  cell_values_plus = get_cell_constraints(f,trian.plus)
#  cell_values_minus = get_cell_constraints(f,trian.minus)
#  lazy_map(BlockMap((2,2),[(1,1),(2,2)]),cell_values_plus,cell_values_minus)
#end

function FESpaces.get_cell_dof_ids(f::MultiFieldFESpace)
  msg = """\n
  This method does not make sense for multi-field
  since each field can be defined on a different triangulation.
  Pass a triangulation in the second argument to get the DOF ids
  on top of the corresponding cells.
  """
  trians = map(get_triangulation,f.spaces)
  trian = first(trians)
  @check all(map(t->t===trian,trians)) msg
  get_cell_dof_ids(f,trian)
end

function FESpaces.get_cell_dof_ids(f::MultiFieldFESpace,trian::Triangulation)
  get_cell_dof_ids(f,trian,MultiFieldStyle(f))
end

#function FESpaces.get_cell_dof_ids(f::MultiFieldFESpace,trian::SkeletonTriangulation)
#  cell_values_plus = get_cell_dof_ids(f,trian.plus)
#  cell_values_minus = get_cell_dof_ids(f,trian.minus)
#  lazy_map(BlockMap(2,[1,2]),cell_values_plus,cell_values_minus)
#end

function FESpaces.get_cell_dof_ids(f::MultiFieldFESpace,::Triangulation,::MultiFieldStyle)
  @notimplemented
end

function FESpaces.get_cell_dof_ids(f::MultiFieldFESpace,trian::Triangulation,::ConsecutiveMultiFieldStyle)
  offsets = compute_field_offsets(f)
  nfields = length(f.spaces)
  blockmask = [ is_change_possible(get_triangulation(Vi),trian) for Vi in f.spaces ]
  active_block_ids = findall(blockmask)
  active_block_data = Any[]
  for i in active_block_ids
    cell_dofs_i = get_cell_dof_ids(f.spaces[i],trian)
    if i == 1
      push!(active_block_data,cell_dofs_i)
    else
      offset = Int32(offsets[i])
      o = Fill(offset,length(cell_dofs_i))
      cell_dofs_i_b = lazy_map(Broadcasting(_sum_if_first_positive),cell_dofs_i,o)
      push!(active_block_data,cell_dofs_i_b)
    end
  end
  lazy_map(BlockMap(nfields,active_block_ids),active_block_data...)
end

function _sum_if_first_positive(a,b)
  if a>0
    return a+b
  else
    return a
  end
end

function Arrays.return_value(k::Broadcasting{typeof(_sum_if_first_positive)},dofs::VectorBlock,o)
  evaluate(k,dofs,o)
end

function Arrays.return_cache(k::Broadcasting{typeof(_sum_if_first_positive)},dofs::VectorBlock,o)
  i = first(findall(dofs.touched))
  dofsi = dofs.array[i]
  ci = return_cache(k,dofsi,o)
  odofsi = evaluate!(ci,k,dofsi,o)
  array = Vector{typeof(odofsi)}(undef,length(dofs.array))
  c = Vector{typeof(ci)}(undef,length(dofs.array))
  for i in 1:length(dofs.array)
    if dofs.touched[i]
      c[i] = return_cache(k,dofs.array[i],o)
    end
  end
  ArrayBlock(array,dofs.touched), c
end

function Arrays.evaluate!(cache,k::Broadcasting{typeof(_sum_if_first_positive)},dofs::VectorBlock,o)
  r,c = cache
  for i in 1:length(dofs.array)
    if dofs.touched[i]
      r[i] = evaluate!(c[i],k,dofs.array[i],o)
    end
  end
  r
end

# API for multi field case

"""
    num_fields(f::MultiFieldFESpace)
"""
function num_fields(f::MultiFieldFESpace)
  length(f.spaces)
end

Base.iterate(m::MultiFieldFESpace) = iterate(m.spaces)

Base.iterate(m::MultiFieldFESpace,state) = iterate(m.spaces,state)

Base.getindex(m::MultiFieldFESpace,field_id::Integer) = m.spaces[field_id]

Base.length(m::MultiFieldFESpace) = length(m.spaces)

# API for the ConsecutiveMultiFieldStyle
import Gridap.FESpaces: interpolate
import Gridap.FESpaces: interpolate_everywhere
import Gridap.FESpaces: interpolate_dirichlet

@deprecate(
  interpolate(fs::MultiFieldFESpace, object),
  interpolate(object, fs::MultiFieldFESpace)
)

@deprecate(
  interpolate_everywhere(fs::MultiFieldFESpace, object),
  interpolate_everywhere(object, fs::MultiFieldFESpace)
)

@deprecate(
  interpolate_dirichlet(fs::MultiFieldFESpace, object),
  interpolate_dirichlet(object, fs::MultiFieldFESpace)
)

"""
The resulting MultiFieldFEFunction is in the space (in particular it fulfills Dirichlet BCs
even in the case that the given cell field does not fulfill them)
"""
function FESpaces.interpolate(objects, fe::MultiFieldFESpace)
  free_values = zero_free_values(fe)
  interpolate!(objects,free_values,fe)
end

function FESpaces.interpolate!(objects,free_values::AbstractVector,fe::MultiFieldFESpace)
  blocks = SingleFieldFEFunction[]
  for (field, (U,object)) in enumerate(zip(fe.spaces,objects))
    free_values_i = restrict_to_field(fe,free_values,field)
    uhi = interpolate!(object, free_values_i,U)
    push!(blocks,uhi)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end

"""
like interpolate, but also compute new degrees of freedom for the dirichlet component.
The resulting MultiFieldFEFunction does not necessary belongs to the underlying space
"""
function FESpaces.interpolate_everywhere(objects, fe::MultiFieldFESpace)
  free_values = zero_free_values(fe)
  blocks = SingleFieldFEFunction[]
  for (field, (U,object)) in enumerate(zip(fe.spaces,objects))
    free_values_i = restrict_to_field(fe,free_values,field)
    dirichlet_values_i = zero_dirichlet_values(U)
    uhi = interpolate_everywhere!(object, free_values_i,dirichlet_values_i,U)
    push!(blocks,uhi)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end

"""
"""
function FESpaces.interpolate_dirichlet(objects, fe::MultiFieldFESpace)
  free_values = zero_free_values(fe)
  blocks = SingleFieldFEFunction[]
  for (field, (U,object)) in enumerate(zip(fe.spaces,objects))
    free_values_i = restrict_to_field(fe,free_values,field)
    dirichlet_values_i = zero_dirichlet_values(U)
    uhi = interpolate_dirichlet!(object, free_values_i,dirichlet_values_i,U)
    push!(blocks,uhi)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end


text:
using IRTools, Test
using IRTools: @meta, IR, CFG, Simple, Loop, Multiple, reloop, explicitbranch!, stackify

relu(x) = (y = x > 0 ? x : 0)

relu_cfg = CFG(@code_ir relu(1))

@test relu_cfg == CFG([[3,2],[4],[4],[]])

@test reloop(relu_cfg) == Simple(1, Multiple([Simple(3), Simple(2)], Simple(4)))

@test sort.(stackify(relu_cfg)) == ([1=>2,1=>3,1=>4],[])

function pow(x, n)
  r = 1
  while n > 0
    n -= 1
    r *= x
  end
  return r
end

pow_cfg = CFG(@code_ir pow(2, 3))

@test pow_cfg == CFG([[2],[4,3],[2],[]])

@test reloop(pow_cfg) == Simple(1, Loop(Simple(2, Simple(3)), Simple(4)))

@test sort.(stackify(pow_cfg)) == ([1=>2,1=>4,2=>3],[3=>2])

# AST recovery

relu(x) = (y = x > 0 ? x : 0)

ir = explicitbranch!(IR(@meta(relu(1)), slots = true))
ex = reloop(ir)

@test eval(:(let arg2 = 5; $ex; end)) == 5
@test eval(:(let arg2 = -5; $ex; end)) == 0

relu(x) = x > 0 ? x : 0

ir = explicitbranch!(IR(@meta(relu(1)), slots = true))
ex = reloop(ir)

@test eval(:(let arg2 = 5; $ex; end)) == 5
@test eval(:(let arg2 = -5; $ex; end)) == 0

ir = explicitbranch!(IR(@meta(pow(1,1)), slots = true))
ex = reloop(ir)

@test eval(:(let arg2 = 5, arg3 = 3; $ex; end)) == 125

ir = explicitbranch!(IR(@meta(gcd(1,1)), slots = true))
ex = reloop(ir)

@test eval(:(let arg2 = 85, arg3 = 391; $ex; end)) == gcd(85, 391)

ext: jl


text:
function construct( e::FundamentalMatrix,
                   𝐊₁::AbstractArray,
                   𝐑₁::AbstractArray,
                   𝐭₁::AbstractArray,
                   𝐊₂::AbstractArray,
                   𝐑₂::AbstractArray,
                   𝐭₂::AbstractArray)

    if size(𝐊₁) != (3,3) || size(𝐊₂) != (3,3) ||
       size(𝐑₁) != (3,3) || size(𝐑₂) != (3,3)
        throw(ArgumentError("Expect 3 x 3 calibration and rotation matrices."))
    end
    if length(𝐭₁) != 3 || length(𝐭₂) != 3
        throw(ArgumentError("Expect length-3 translation vectors."))
    end
    𝐅 = vec2antisym(𝐊₂*𝐑₂*(𝐭₁ .- 𝐭₂))*𝐊₂*𝐑₂/𝐑₁/𝐊₁
    SMatrix{3,3,Float64,3*3}(𝐅)
end

function construct( e::FundamentalMatrix, 𝐏₁::AbstractArray, 𝐏₂::AbstractArray)
    if (size(𝐏₁) != (3,4)) || (size(𝐏₂) != (3,4))
        throw(ArgumentError("Expect 3 x 4 projection matrices."))
    end
    𝐜₁ = SVector{4,Float64}(nullspace(Array(𝐏₁)))
    𝐞₂ = 𝐏₂*𝐜₁
    𝐅 = vec2antisym(𝐞₂)*𝐏₂*pinv(Array(𝐏₁))
    SMatrix{3,3,Float64,3*3}(𝐅)
end


text:
module IDEDelay

    function ntrp3h(tint,t,y,yp,tnew,ynew,ypnew)
    # Hermite extrapolation
    h     = tnew - t;
    s     = (tint - t)/h;
    s2    = s * s;
    s3    = s * s2;
    slope = (ynew - y)/h;
    c     = 3*slope - 2*yp - ypnew;
    d     = yp + ypnew - 2*slope;

    return y + (h*d*s3 + h*c*s2 + h*yp*s);
	end

	function int_simpson(h, y_begin, y_end, y_half)
    # Simpson's method
    return h/6*(y_begin + 4*y_half + y_end);
	end

	function ide_delay_rk(idefun,Core,delays,history,tspan,step)
    #  idefun - right-hand side function
    #    Core - Kernel (integrated function)
    #  delays - delays function (lower integration limit)
    # History - history function
    #   tspan - solution interval
    #    step - step-size

    t0   = tspan[1];    # t begin
    tf   = tspan[2];    # t end
    y0   = history(t0); # initial solution
    neq  = length(y0);  # number of equations
    htry = step;        # constant mesh step-size
    h    = htry;        # current step-size (for the last step)
    #==============================================================#
    # VIDE Runge-Kutta Bel'tyukov-type of Order 4
    A = [0 1/2   0  0
         0   0 1/2  0
         0   0   0  1
         0   0   0  0 ];

    b = [1/6; 1/3; 1/3; 1/6];
    s = length(b);
    c = [0 1/2 1/2 1];
    d = [1/2; 1/2; 1];
    #==============================================================#
    # Calculate integral (F) in history
    F = zeros(neq);

    # Piece from delays(t0) to the next mesh point
    # tj, tj_half, tj_1, Core_tj, Core_tj_1, Core_tj_h are always needed
    # to use Simpson's Method for numerical approximations of integrals
    tj      = delays(t0);             # Begin
    step    = trunc(Int,(t0 - tj)/h); # The number of memorized
                                      # intervals of history
    tj_1    = t0 - step*h;            # End
    tj_half = (tj+tj_1)/2;            # Half segment

    # Calculate Kernel values at the nodes
    Core_tj   = Core(t0,tj, history(tj));
    Core_tj_1 = Core(t0,tj_1, history(tj_1));
    Core_tj_h = Core(t0,tj_half, history(tj_half));

    # Simpson's method
    F   = F + ones(neq).*int_simpson(tj_1-tj,Core_tj,Core_tj_1,Core_tj_h);

    # Main integral over the mesh points
    Core_tj = Core_tj_1;
    for j = step-1:-1:0
        tj_1      = t0 - j*h;
        tj_half   = tj_1 - h/2;

        # Calculate Kernel values at the nodes
        Core_tj_h = Core(t0,tj_half, history(tj_half));
        Core_tj_1 = Core(t0,tj_1, history(tj_1));

        # Simpson's method
        F         = F + ones(neq).*int_simpson(h,Core_tj,Core_tj_1,Core_tj);

        # Kernel of tj_1 is equal to tj of the next step
        Core_tj   = Core_tj_1;
    end
    #==============================================================#
    # Initialization | First Step | Y | K
    t      = t0;
    y      = zeros(neq,1);
    y[:,1] = ones(neq).*y0;
    k      = 1;  # step

    Y        = zeros(neq,s);
    K        = zeros(neq,s,1);
    K[:,1,k] = idefun(t[k],y[:,1],F);
    Core_di  = zeros(neq,s);
    #==============================================================#

    while t[k] < tf
        #==========================================================#
        # Last step
        if t[k] + h > tf
            h = tf - t[k];
        end
        #==========================================================#
        Z      = zeros(neq,s);
        Y[:,1] = y[:,k];

        # Runge-Kutta steps
        for i = 2 : s
            ti = t[k] + c[i]*h;
            #======================================================#
            # Calculate integral (F)
            if i != 3 # c[3] = c[2] so the same F value is used
                F = zeros(neq);

                dtk_begin = delays(ti); # lower integration limit
                if dtk_begin < t0
                    #==============================================#
                    # Integral begins in the history

                    # Step of delays(ti) in the history
                    step = trunc(Int,(t0 - dtk_begin)/h);

                    # Add piece from dtk_begin to the next mesh point in the history
                    tj        = dtk_begin;
                    tj_1      = t0-step*htry;
                    tj_half   = (tj+tj_1)/2;

                    # Calculate Kernel values at the nodes
                    Core_tj   = Core(ti,tj, history(tj));
                    Core_tj_1 = Core(ti,tj_1, history(tj_1));
                    Core_tj_h = Core(ti,tj_half, history(tj_half));

                    # Simpson's method
                    F = F + ones(neq).*int_simpson(tj_1-tj,Core_tj,Core_tj_1,Core_tj_h);

                    # Main integral in the history
                    Core_tj = Core_tj_1;
                    for j = step-1:-1:0
                        tj_1      = t0 - j*htry;
                        tj_half   = tj_1 - htry/2;

                        # Calculate Kernel values at the nodes
                        Core_tj_1 = Core(ti,tj_1, history(tj_1));
                        Core_tj_h = Core(ti,tj_half, history(tj_half));

                        # Simpson's method
                        F = F + ones(neq).*int_simpson(htry,Core_tj,Core_tj_1,Core_tj_h);

                        # Kernel of tj_1 is equal to tj of the next step
                        Core_tj   = Core_tj_1;
                    end

                    # Add integral in the solution to t(k)
                    for j = 2:k
                        tj_half   = t[j] - h/2;
                        y_half    = ntrp3h(tj_half,t[j-1],y[:,j-1],K[:,1,j-1],
                                                   t[j],y[:,j],K[:,1,j]);

                        # Calculate Kernel values at the nodes
                        if neq == 1
                            Core_tj_h = Core(ti,tj_half, y_half[1]);
                            Core_tj_1 = Core(ti,t[j], y[1,j]);
                        else
                            Core_tj_h = Core(ti,tj_half, y_half);
                            Core_tj_1 = Core(ti,t[j], y[:,j]);
                        end

                        # Simpson's method
                        F = F + ones(neq).*int_simpson(htry,Core_tj,Core_tj_1,Core_tj_h);

                        # Kernel of tj_1 is equal to tj of the next step
                        Core_tj   = Core_tj_1;
                    end
                    #==============================================#
                else
                    #==============================================#
                    # Integral only over the solution

                    # Step of delays(ti) in the solution
                    step      = trunc(Int,(dtk_begin-t0)/htry + 1);

                    # Add piece from dtk_begin to the mesh point in the solution
                    tj_half   = (t[step+1] + dtk_begin)/2;

                    y_begin   = ntrp3h(dtk_begin,t[step],y[:,step],K[:,1,step],
                                           t[step+1],y[:,step+1],K[:,1,step+1]);
                    y_begin_h = ntrp3h(tj_half,t[step],y[:,step],K[:,1,step],
                                           t[step+1],y[:,step+1],K[:,1,step+1]);

                    # Calculate Kernel values at the nodes
                    if neq == 1
                        Core_tj   = Core(ti,dtk_begin, y_begin[1]);
                        Core_tj_1 = Core(ti,t[step+1], y[1,step+1]);
                        Core_tj_h = Core(ti,tj_half, y_begin_h[1]);
                    else
                        Core_tj   = Core(ti,dtk_begin, y_begin);
                        Core_tj_1 = Core(ti,t[step+1], y[:,step+1]);
                        Core_tj_h = Core(ti,tj_half, y_begin_h);
                    end

                    # Simpson's method
                    F = F + ones(neq).*int_simpson(t[step+1]-dtk_begin,Core_tj,
                                                  Core_tj_1,Core_tj_h);

                    # Main integral to t(k)
                    Core_tj = Core_tj_1;
                    for j = step+2:k
                        tj_half   = t[j] - htry/2;
                        y_half    = ntrp3h(tj_half,t[j-1],y[:,j-1],K[:,1,j-1],
                                                   t[j],y[:,j],K[:,1,j]);

                        # Calculate Kernel values at the nodes
                        if neq == 1
                            Core_tj_h = Core(ti,tj_half, y_half[1]);
                            Core_tj_1 = Core(ti,t[j], y[1,j]);
                        else
                            Core_tj_h = Core(ti,tj_half, y_half);
                            Core_tj_1 = Core(ti,t[j], y[:,j]);
                        end

                        # Simpson's method
                        F = F + ones(neq).*int_simpson(htry,Core_tj,Core_tj_1,Core_tj_h);

                        # Kernel of tj_1 is equal to tj of the next step
                        Core_tj   = Core_tj_1;
                    end
                    #==============================================#
                end
            end
            #======================================================#
            #println("i = ",i)
            # Y2-4
            Y[:,i] = y[:,k] + h * (K[:,1:i-1,k] * A[1:i-1,i]);

            # Z2-4
            if neq == 1
                Core_di[1,i-1] = Core(t[k]+d[i-1]*h,t[k]+c[i-1]*h,Y[1,i-1]);
            else
                Core_di[:,i-1] = Core(t[k]+d[i-1]*h,t[k]+c[i-1]*h,Y[:,i-1]);
            end
            Z[:,i]         = h * (Core_di[:,1:i-1] * A[1:i-1,i]);

            # K2-4
            K[:, i, k] = idefun(ti, Y[:,i], F+Z[:,i]);
        end
        #==========================================================#
        # Final approximation of RK Method
        t = [t; t[k] + h];
        y = [y y[:,k] + h * ( K[:,:,k] * b )];
        #==========================================================#
        # Calculate K(1) for next step
        # Hermite extrapolation for K(1,k+1)
        y_k_half     = 3/4*y[:,k] + 1/4*y[:,k+1] + h/4*K[:,1,k];
        if neq == 1
            Core_tk_half = Core(t[k+1], t[k]+h/2, y_k_half[1]);
            Core_tk      = Core(t[k+1], t[k+1], y[1,k+1]);
        else
            Core_tk_half = Core(t[k+1], t[k]+h/2, y_k_half);
            Core_tk      = Core(t[k+1], t[k+1], y[:,k+1]);
        end

        F  = F + ones(neq).*int_simpson(h,Core_tj,Core_tk,Core_tk_half);
        K1 = zeros(neq,s,1);
        K  = cat(K, K1, dims=3);
        K[:,1,k+1] = idefun(t[k+1], y[:,k+1], F);

        k = k + 1;
    end
    return t, y';
end

end # module


text:
module CUJulia

using Base.LinAlg.BLAS

if is_windows()
    const libcuda = Libdl.find_library(["nvcuda"])
else
    const libcuda = Libdl.find_library(["libcuda"])
end
isempty(libcuda) && throw("CUDA driver library cannot be found.")

function check_curesult(status)
    status == CUDA_SUCCESS && return
    p = Ptr{UInt8}[0]
    cuGetErrorString(status, p)
    throw(unsafe_string(p[1]))
end

p = Cint[0]
ccall((:cuDriverGetVersion,libcuda), UInt32, (Ptr{Cint},), p)
const driver_version = Int(p[1])
info("CUDA driver version: $driver_version")

const major = div(driver_version, 1000)
const minor = div(driver_version - major*1000, 10)
include("../lib/$(major)$(minor)/libcuda.jl")
include("../lib/$(major)$(minor)/libcuda_types.jl")

include("device.jl")
include("function.jl")
initctx()
infodevices()

ctype(::Type{Int64}) = :int
ctype(::Type{Float32}) = :float
ctype(::Type{Float64}) = :double

include("pointer.jl")
include("array.jl")
include("subarray.jl")
#include("math.jl")
#include("broadcast.jl")
#include("reducedim.jl")

include("Interop.jl")
include("NVRTC.jl")
include("CUBLAS.jl")
include("CUDNN.jl")

end


text:
"""
# Summary
  A Structure for the irfs/simulations of DRIPs

# Fields
    T     : length of IRFs/simulation
    x     : IRFs/simulation of the fundamental shocks
    x_hat : IRFs/simulation of beliefs
    a     : IRFs/simulation of actions

In particular, if `n` is the
dimension of `x`, `m` is the dimension of `a` and `k` is the number of
structural shocks, then

* `x` has dimension `n*k*T` where `x(i,j,:)` is the impulse response function of
    the `i`'th dimension of `x` to the `j`'th structural shock.
* `x_hat` has dimension `n*k*T` where `x_hat(i,j,:)` is the impulse response
    function of the agent's average belief about the `i`'th dimension of `x` to
    the `j`'th structural shock.
* `a` has dimension `m*k*T` where `a(i,j,:)` is the impulse response function of
    the `i`'th action to the `j`'th structural shock.
"""
struct Path
    T     :: Integer        # length of IRFs
    x     :: Array{Float64} # IRFs of the fundamental shocks
    x_hat :: Array{Float64} # IRFs of beliefs
    a     :: Array{Float64} # IRFs of actions
end

"""
              irfs(p :: Drip;        # Steady state of the DRIP
                   T :: Integer = 40   # Optional: length of impulse response functions
                   ) -> Path
Returns a [`Path`](@ref) structure with the impulse response functions of the fundamental (`x`), beliefs (`x_hat`)
    and actions (`a`) to all the structural shocks
    **under the steady state information strucutre**. In particular, if `n` is the
    dimension of `x`, `m` is the dimension of `a` and `k` is the number of
    structural shocks, then

* `x` has dimension `n*k*T` where `x(i,j,:)` is the impulse response function
    of the `i`'th dimension of `x` to the `j`'th structural shock.
* `x_hat` has dimension `n*k*T` where `x_hat(i,j,:)` is the impulse response function
    of the agent's average belief about the `i`'th dimension of `x` to the `j`'th
    structural shock.
* `a` has dimension `m*k*T` where `a(i,j,:)` is the impulse response function
    of the `i`'th action to the `j`'th structural shock.
"""
function irfs(p :: Drip;        # Steady state of the DRIP
              T :: Integer = 40   # Optional: length of impulse response functions
              )
    (n,m) = size(p.H)
    (_,k) = size(p.Q)
    x     = zeros(n,k,T);
    x_hat = zeros(n,k,T);
    a     = zeros(m,k,T);
    for kk in 1:k
        e_k = zeros(k,1); e_k[kk] = 1;
        for ii in 1:T
            if ii==1
                x[:,kk,ii]     = p.Q*e_k;
                x_hat[:,kk,ii] = (p.ss.K*p.ss.Y')*(x[:,kk,ii]);
            else
                x[:,kk,ii]     =p.A*x[:,kk,ii-1];
                x_hat[:,kk,ii] =p.A*x_hat[:,kk,ii-1]+(p.ss.K*p.ss.Y')*(x[:,kk,ii]-p.A*x_hat[:,kk,ii-1]);
            end
            a[:,kk,ii]  .= p.H'*x_hat[:,kk,ii];
        end
    end
    return(Path(T,x,x_hat,a))
end

## get irfs with optimal signals acquired under pt::Trip
"""
          irfs(pt :: Trip;       # Transition dynamics of the DRIP
               T  :: Integer = 40  # Optional: length of impulse response functions
               ) -> Path
Returns a [`Path`](@ref) structure with the impulse response functions of the fundamental (`x`), beliefs (`x_hat`)
    and actions (`a`) to all the structural shocks
    **under the information structure implied by** `P`. In particular, if `n` is the
    dimension of `x`, `m` is the dimension of `a` and `k` is the number of
    structural shocks, then

* `x` has dimension `n*k*T` where `x(i,j,:)` is the impulse response function
    of the `i`'th dimension of `x` to the `j`'th structural shock.
* `x_hat` has dimension `n*k*T` where `x_hat(i,j,:)` is the impulse response function
    of the agent's average belief about the `i`'th dimension of `x` to the `j`'th
    structural shock.
* `a` has dimension `m*k*T` where `a(i,j,:)` is the impulse response function
    of the `i`'th action to the `j`'th structural shock.
"""
function irfs(pt :: Trip;       # Transition dynamics of the DRIP
              T  :: Integer = 40  # Optional: length of impulse response functions
              )
    p     = pt.p; # steady state Drip
    (n,m) = size(p.H)
    (_,k) = size(p.Q)
    L     = size(pt.Ds,2) # length of Trip

    eye   = Matrix(I,n,n);
    x     = zeros(n,k,T); # initialize IRFs of state
    x_hat = zeros(n,k,T); # initialize IRFs of beliefs
    a     = zeros(m,k,T); # initialize IRFs of actions
    for kk in 1:k
        e_k = zeros(k,1); e_k[kk] = 1;
        for ii in 1:T
            if ii==1
                x[:,kk,ii]     = p.Q*e_k;
                x_hat[:,kk,ii] = (eye-pt.Σ_ps[:,:,ii]*pinv(pt.Σ_1s[:,:,ii]))*(x[:,kk,ii]);
            elseif ii <= L
                x[:,kk,ii]     =p.A*x[:,kk,ii-1];
                x_hat[:,kk,ii] =p.A*x_hat[:,kk,ii-1]+(eye-pt.Σ_ps[:,:,ii]*pinv(pt.Σ_1s[:,:,ii]))*(x[:,kk,ii]-p.A*x_hat[:,kk,ii-1]);
            else
                x[:,kk,ii]     =p.A*x[:,kk,ii-1];
                x_hat[:,kk,ii] =p.A*x_hat[:,kk,ii-1]+(eye-pt.Σ_ps[:,:,end]*pinv(pt.Σ_1s[:,:,end]))*(x[:,kk,ii]-p.A*x_hat[:,kk,ii-1]);
            end
            a[:,kk,ii]  .= p.H'*x_hat[:,kk,ii];
        end
    end
    return(Path(T,x,x_hat,a))
end

## get irfs with information treatment S at time zero.
"""
         irfs(p          :: Drip,                          # Steady state of the DRIP (when treatment happens)
              S          :: Signal;                        # Signal for treatment
              T          :: Integer = 40,                    # optional: length of irfs
              reoptimize :: Bool = true,                   # if true gives the irfs with reoptimized signals, if false with steady state signals
              trip       :: Union{Trip, Nothing} = nothing # if false solves the trip, if = P::trip then takes P as the trip
              )
Returns a [`Path`](@ref) structure with the impulse response functions of the fundamental (`x`), beliefs (`x_hat`)
    and actions (`a`) to all the structural shocks
    **under the information structure implied by a one time information treatment
    with** `S` **in the steady state of the DRIP** `P`. In particular, if `n` is the
    dimension of `x`, `m` is the dimension of `a` and `k` is the number of
    structural shocks, then

* `x` has dimension `n*k*T` where `x(i,j,:)` is the impulse response function
    of the `i`'th dimension of `x` to the `j`'th structural shock.
* `x_hat` has dimension `n*k*T` where `x_hat(i,j,:)` is the impulse response function
    of the agent's average belief about the `i`'th dimension of `x` to the `j`'th
    structural shock.
* `a` has dimension `m*k*T` where `a(i,j,:)` is the impulse response function
    of the `i`'th action to the `j`'th structural shock.
"""
function irfs(p          :: Drip,                          # Steady state of the DRIP (when treatment happens)
              S          :: Signal;                        # Signal for treatment
              T          :: Integer = 40,                    # optional: length of irfs
              reoptimize :: Bool = true,                   # if true gives the irfs with reoptimized signals, if false with steady state signals
              trip       :: Union{Trip, Nothing} = nothing # if false solves the trip, if = P::trip then takes P as the trip
              )

    (n,m) = size(p.H)
    (_,k) = size(p.Q)
    eye   = Matrix(I,n,n);

    if reoptimize == true
        if trip == nothing
            trip = Trip(p,S;T = T);
        end
        L     = size(trip.Ds,2) # length of Trips
    end

    x     = zeros(n,k,T); # initialize IRFs of state
    x_hat = zeros(n,k,T); # initialize IRFs of beliefs
    a     = zeros(m,k,T); # initialize IRFs of actions
    for kk in 1:k
        K0   = p.ss.Σ_1*S.L/(S.L'*p.ss.Σ_1*S.L+S.Σ_z); # Get prior covariance matrix from new signal
        Σ0   = p.ss.Σ_1-K0*S.L'*p.ss.Σ_1;
        e_k = zeros(k,1); e_k[kk] = 1;
        for ii in 1:T
            if ii==1
                x[:,kk,ii]     = p.Q*e_k;
                x_hat_p        = K0*S.L'*(x[:,kk,ii]);
                if reoptimize == true
                    x_hat[:,kk,ii] = x_hat_p + (eye-trip.Σ_ps[:,:,ii]*pinv(trip.Σ_1s[:,:,ii]))*(x[:,kk,ii]-x_hat_p);
                else
                    K0 = Σ0*p.ss.Y/(p.ss.Y'*Σ0*p.ss.Y+p.ss.Σ_z);
                    Σ0 = p.Q*p.Q' + p.A*(Σ0-K0*p.ss.Y'*Σ0)*p.A';
                    x_hat[:,kk,ii] = x_hat_p + (K0*p.ss.Y')*(x[:,kk,ii]-x_hat_p);
                end
            else
                x[:,kk,ii]     =p.A*x[:,kk,ii-1];
                if reoptimize  == true
                    if ii <= L
                        x_hat[:,kk,ii] =p.A*x_hat[:,kk,ii-1]+(eye-trip.Σ_ps[:,:,ii]*pinv(trip.Σ_1s[:,:,ii]))*(x[:,kk,ii]-p.A*x_hat[:,kk,ii-1]);
                    else
                        x_hat[:,kk,ii] =p.A*x_hat[:,kk,ii-1]+(eye-trip.Σ_ps[:,:,end]*pinv(trip.Σ_1s[:,:,end]))*(x[:,kk,ii]-p.A*x_hat[:,kk,ii-1]);
                    end
                else
                    K0 = Σ0*p.ss.Y/(p.ss.Y'*Σ0*p.ss.Y+p.ss.Σ_z);
                    Σ0 = p.Q*p.Q' + p.A*(Σ0-K0*p.ss.Y'*Σ0)*p.A';
                    x_hat[:,kk,ii] =p.A*x_hat[:,kk,ii-1]+(K0*p.ss.Y')*(x[:,kk,ii]-p.A*x_hat[:,kk,ii-1]);
                end
            end
            a[:,kk,ii]  .= p.H'*x_hat[:,kk,ii];
        end
    end
    return(Path(T,x,x_hat,a))
end


text:
ENV["TRAVIS_REPO_SLUG"] = "github.com/TAJD/sail_route.jl.git"
ENV["TRAVIS_PULL_REQUEST"] = "false"
ENV["TRAVIS_OS_NAME"] = "osx"
ENV["TRAVIS_JULIA_VERSION"] = "nightly"
ENV["TRAVIS_BRANCH"]="master"


text:
using WriteVTK

#=
This is the 1D WriteMesh routine
=#
function writemesh(base_name, x1; fields=(), realelems=1:size(x1)[end])
  (Nqr, _) = size(x1)
  Nsubcells = (Nqr-1)

  cells = Array{MeshCell{Array{Int,1}}, 1}(undef, Nsubcells * length(realelems))
  for e ∈ realelems
    offset = (e-1) * Nqr
    for i = 1:Nqr-1
        cells[i + (e-1)*Nsubcells] =
          MeshCell(VTKCellTypes.VTK_LINE, offset .+ [i, i+1])
    end
  end

  vtkfile = vtk_grid("$(base_name)", @view(x1[:]), cells; compress=false)
  for (name, v) ∈ fields
    vtk_point_data(vtkfile, v, name)
  end
  outfiles = vtk_save(vtkfile)
end

#=
This is the 2D WriteMesh routine
=#
function writemesh(base_name, x1, x2; x3=nothing, fields=(), realelems=1:size(x1)[end])
  @assert size(x1) == size(x2)
  (Nqr, Nqs, _) = size(x1)
  Nsubcells = (Nqr-1) * (Nqs-1)

  cells = Array{MeshCell{Array{Int,1}}, 1}(undef, Nsubcells * length(realelems))
  ind = LinearIndices((1:Nqr, 1:Nqs))
  for e ∈ realelems
    offset = (e-1) * Nqr * Nqs
    for j = 1:Nqs-1
      for i = 1:Nqr-1
        cells[i + (j-1)*(Nqr-1) + (e-1)*Nsubcells] =
        MeshCell(VTKCellTypes.VTK_PIXEL, offset .+ ind[i:i+1,j:j+1][:])
      end
    end
  end

  if x3 == nothing
    vtkfile = vtk_grid("$(base_name)", @view(x1[:]), @view(x2[:]), cells;
                       compress=false)
  else
    vtkfile = vtk_grid("$(base_name)", @view(x1[:]), @view(x2[:]), @view(x3[:]),
                       cells; compress=false)
  end
  for (name, v) ∈ fields
    vtk_point_data(vtkfile, v, name)
  end
  outfiles = vtk_save(vtkfile)
end

#=
This is the 3D WriteMesh routine
=#
function writemesh(base_name, x1, x2, x3; fields=(), realelems=1:size(x1)[end])
  (Nqr, Nqs, Nqt, _) = size(x1)
  (Nr, Ns, Nt) = (Nqr-1, Nqs-1, Nqt-1)
  Nsubcells = Nr * Ns * Nt
  cells = Array{MeshCell{Array{Int,1}}, 1}(undef, Nsubcells * length(realelems))
  ind = LinearIndices((1:Nqr, 1:Nqs, 1:Nqt))
  for e ∈ realelems
    offset = (e-1) * Nqr * Nqs * Nqt
    for k = 1:Nt
      for j = 1:Ns
        for i = 1:Nr
          cells[i + (j-1) * Nr + (k-1) * Nr * Ns + (e-1) * Nsubcells] =
            MeshCell(VTKCellTypes.VTK_VOXEL,
                     offset .+ ind[i:i+1, j:j+1, k:k+1][:])
        end
      end
    end
  end

  vtkfile = vtk_grid("$(base_name)", @view(x1[:]), @view(x2[:]), @view(x3[:]),
                     cells; compress=false)
  for (name, v) ∈ fields
    vtk_point_data(vtkfile, v, name)
  end
  outfiles = vtk_save(vtkfile)
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
Base.load_juliarc()
include("MxArrays.jl")
using MxArrays
include("Mex.jl")
using Mex


text:
using Pkg
Pkg.activate(@__DIR__)
using Flux, Statistics, Random, Distributions, StatsFuns, Zygote, CUDA, ChainRulesCore

function ChainRulesCore.rrule(::typeof(cpu), x::CuArray{Float32})
    pullback(Δ::Array{Float32}) = NoTangent(), cu(Δ)
    return cpu(x), pullback
end
use_cuda = false #cpu by default

function forward(m, x)
    #this just adds a softmax to the policy output
    if use_cuda
        y = m( Zygote.@ignore cu(x) ) |> cpu #RNN step
    else
        y = m(x)
    end
    logπ = y[1:Naction]
    logπ = logπ .- StatsFuns.logsumexp(logπ) #softmax
    return [logπ; y[Naction+1]]
end

function calc_deltas(rews, Vs)
    #compute RPEs
    R = 0
    δs = []
    for t = 0:(T-1) #compute RPE
        R = rews[T-t] + γ*R #Reward
        δs = [R - Vs[T-t]; δs] #error
    end
    return Float32.(δs)
end
Zygote.@nograd calc_deltas #don't take gradients of this

function run_episode(m; ps = zeros(2), state = zeros(3), hidden = true)
    ### initialize reward probabilities and state ###
    ps, state = initialize(ps, state)

    #arrays for storing output
    #ys = Matrix{Float32}(undef, Nout,0)
    ys = Matrix{Float32}(undef, Nout,0)
    if hidden
        hs = Matrix{Float32}(undef, Nhidden, 0)
        states = Matrix{Float32}(undef, Nstate_rep, 0) #Nstate_rep is the non-one-hot size
    end
    as = Matrix{Int32}(undef, 1,0)
    rews = Matrix{Float32}(undef, 1,0)

    ### reset model and initialize input!! ###
    Flux.reset!(m)
    x = Float32.(zeros(Nin))
    for t = 1:T
        if hidden
            hs = [hs m[1].state[1]] #store hidden state
            states = [states state] #store state
        end

        y = forward(m, x) #RNN step
        a, rew, x, ps, state = env(y, t, state, ps) #output -> action, reward, next input

        ys = [ys y] #store output
        rews = [rews rew] #store reward
        as = [as a] #store action
    end

    δs = calc_deltas(rews, ys[Naction+1, :])
    L = Float32(0.)
    for t = 1:T
        logπ=ys[1:Naction, t]
        H = -sum(exp.(logπ) .* logπ)
        RPE_term = δs[t]*(ys[as[t], t] + βv*ys[Naction+1, t])
        L -= (RPE_term + βe*H) #loss
    end

    hidden && return L, ys, rews[1,:], as[1,:], ps, hs, states#, xs
    return L, ys, rews[1,:], as[1,:], ps
end

function model_loss((), ())
    #wrapper for Flux that takes empty data
    loss, _, _, _, _ = run_episode(m, hidden = false)
    return loss
end

max_stars_repo_name: KrisJensen/metaRL

max_issues_repo_name: KrisJensen/metaRL

max_forks_repo_name: KrisJensen/metaRL


text:
module ProjectFlow

const PF_ROOT = joinpath(homedir(), ".projectflow/")

include("loader.jl")
include("manager.jl")
include("project.jl")
include("logger.jl")
include("template.jl")

export Project, initiate

"""
function initiate(p::Project)

This is the main function compose all other function together and
upon giving value of type project function will validate the
initialisers and create a new project or load a existing project.

# Argument
- p: value of type Project

# Example
```
p = Project(
    id="xyz",
    name="My Fancy? *Project1 2 ",
    template="jl",
    profile="default"
)
initiate(p)
```
"""
function initiate(p::Project)
    properties = select_property(PF_ROOT, p.profile)
    set_workflow(properties)
    meta_info = consolidate(p)
    name = meta_info["full_name"]
    pexist, datalake, iviz, idata = pathfinder(name, properties)
    if !pexist
        log_project(PF_ROOT, meta_info)
        pmap = Dict(
            "i" => p.id,
            "n" => p.name,
            "t" => p.template,
            "p" => p.profile,
        )
        p = project_path(name, properties)
        save_template(p, pmap)
    end
    return datalake, iviz, idata
end

end # module


text:
@testset "Test constrained NLP vector score" begin
###############################################################################
#
# The Stopping structure eases the implementation of algorithms and the
# stopping criterion.
# We illustrate here the basic features of Stopping.
#
# -> the case where the score is a Vector and not a Number.
# a) tol_check is a Number;
# b) tol_check is a Vector of the same size.
# Note that if tol_check is a Vector and the score is a Number, only
# the smallest value of tol_check is considered.
#
#Warning: see https://github.com/JuliaSmoothOptimizers/NLPModels.jl/blob/master/src/autodiff_model.jl
#for the proper way of defining an ADNLPModel
#
###############################################################################

#We consider a regularized MPCC problem:
x0 = ones(2)
t  = 1.0
c(x) = [(x[1] - t) .* (x[2] - t)]
meta = NLPModelMeta(2, x0=x0, lvar = zeros(2), uvar = Inf*ones(2),
                    ncon = 1, y0=zeros(1), lcon = -Inf*ones(1), ucon = zeros(1))
nlp = ADNLPModel(meta, Counters(), ADNLPModels.ForwardDiffAD(1, 1), x -> x[1]+x[2], c)

#We consider a vectorized optimality_check KKT function:
function KKTvect(pb    :: AbstractNLPModel,
                 state :: NLPAtX;
                 kwargs...)

    #Check the gradient of the Lagrangian
    gLagx      = Stopping._grad_lagrangian(pb, state)
    #Check the complementarity condition for the bounds
    dual_res_bounds = Stopping._sign_multipliers_bounds(pb, state)
    #Check the complementarity condition for the constraints
    res_nonlin = Stopping._sign_multipliers_nonlin(pb, state)
    #Check the feasibility
    feas       = Stopping._feasibility(pb, state)

    res = vcat(gLagx, feas, dual_res_bounds, res_nonlin)

    return res
end

#One solution is [1.0, 0.0] while [1. + 1e-7, 0] is an approximate solution.
sol  = [1.0, 0.0]
esol = [1.0 + 1e-7, 0.0]

#As usual, we initialize the State and the Stopping
nlp_at_x_c = NLPAtX(x0, NaN*ones(nlp.meta.ncon), Array{Float64,1}(undef, 14))
stop_nlp = NLPStopping(nlp, nlp_at_x_c, optimality_check = KKTvect)
fill_in!(stop_nlp, x0)
OK = stop!(stop_nlp)
@test !OK #the initial guess (1.0, 1.0) is not a solution
reinit!(stop_nlp, rstate = true) #reinitialize the Stopping and the State
fill_in!(stop_nlp, sol)
OK = stop!(stop_nlp)
@test OK #sol is a solution
reinit!(stop_nlp, rstate = true)
fill_in!(stop_nlp, esol)
OK = stop!(stop_nlp)
@test OK #esol is also an (approximate) solution.
#KKTvect(stop_nlp.pb, stop_nlp.current_state)

#Our aim is to verify that a solution satisfies the strong epsilon-stationarity
#condition, i.e. the feasibility and the complementarity condition have to be
#satisfied exactly.
tol_check1(atol,rtol,opt0) = vcat(max(atol,rtol*opt0) .* ones(6), zeros(8))

stop_nlp = NLPStopping(nlp, nlp_at_x_c, optimality_check = KKTvect, tol_check = tol_check1)
fill_in!(stop_nlp, esol)
OK = stop!(stop_nlp)
@test !OK #esol is not a solution in this case as it satisfies the feasibility approximately
reinit!(stop_nlp, rstate = true) #reinitialize the Stopping and the State
fill_in!(stop_nlp, sol)
OK = stop!(stop_nlp)
@test OK #However, sol is a solution.

#We can also use an asymetric test on the optimality condition using tol_check_neg
#by default tol_check_neg = - tol_check
tol_check1_neg(atol,rtol,opt0) = - 1e-3 * ones(14)
stop_nlp = NLPStopping(nlp, nlp_at_x_c, optimality_check = KKTvect,
                       tol_check = tol_check1, tol_check_neg = tol_check1_neg)
fill_in!(stop_nlp, sol)
OK = stop!(stop_nlp)
@test OK

#An error is returned if size(KKT) is different from size(tol_check).
try
   stop_nlp2 = NLPStopping(nlp, nlp_at_x_c, optimality_check = KKT, tol_check = tol_check1)
   fill_in!(stop_nlp2, esol)
   stop!(stop_nlp2)
   @test false
catch
   @test true
end

end


text:
begin
# Jump Search

# Julia code to jump search x in arr[].
# If x is present then return its location,
# otherwise return -1
import Base:length,show,AbstractArray
arr = [ 0, 1, 1, 2, 3, 5, 8, 13, 21,
    34, 55, 89, 144, 233, 377, 610 ]
x = 610
n = length(arr)
function jumpsearch(arr::AbstractArray{},n,x)

    # Finding block size to be jumped
    step = sqrt(n)

    # Finding the block where element is
    # present (if it is present)
    prev = 1
    while arr[round(Int,floor(min(step, n)))] < x
        prev = step
        step += sqrt(n)
        if prev >= n
            return -1
		end
	end
    # Doing a linear search for x in
    # block beginning with prev.
    while arr[round(Int,floor(prev))] < x
        prev += 1

        # If we reached next block or end
        # of array, element is not present.
        if prev == min(step, n)
            return -1
		end
	end
    # If element is found
    if arr[round(Int,floor(prev))] == x
        return prev
	end
    return -1
end
global result=-1

result=jumpsearch(arr,n, x)

println("Result after function execution:",result)
if result == -1
    println("Element:",x," is not present in array",arr)
else
    println("Element:",x," is present at index:", result," in array",arr)
end
end

max_stars_repo_name: Learning-Julia/Algorithms

max_issues_repo_name: Learning-Julia/Algorithms

max_forks_repo_name: Learning-Julia/Algorithms


text:
module Gate

using PyCall
import ..@pyfunc
import ..@pyclass

# qulacs = PyNULL()
gate = PyNULL()

const gate_classes = [
    :Adaptive, :RX,
    :AmplitudeDampingNoise, :RY,
    :BitFlipNoise, :RZ,
    :CNOT, :RandomUnitary,
    :CP, :ReversibleBoolean,
    :CPTP, :S,
    :CZ, :SWAP,
    :DenseMatrix, :Sdag,
    :DephasingNoise, :SparseMatrix,
    :DepolarizingNoise, :StateReflection,
    :DiagonalMatrix, :T,
    :FREDKIN, :TOFFOLI,
    :H, :Tdag,
    :Identity, :TwoQubitDepolarizingNoise,
    :IndependentXZNoise, :U1,
    :Instrument, :U2,
    :Measurement, :U3,
    :P0, :X,
    :P1, :Y,
    :ParametricPauliRotation, :Z,
    :ParametricRX,               #:add,
    :ParametricRY,               #:merge,
    :ParametricRZ,               #:sqrtX,
    :Pauli,                      #:sqrtXdag,
    :PauliRotation,              #:sqrtY,
    :Probabilistic,              #:sqrtYdag,
    :ProbabilisticInstrument,    #:to_matrix_gate,
]

#=
for class in gate_classes
    @eval begin
        struct $class
            pyobj::PyObject
            $(class)(args...) = new(qulacs.gate.$(class)(args...))
        end

        PyObject(t::$(class)) = t.pyobj

        function Base.propertynames(t::$(class))
            propertynames(getfield(t, :pyobj))
        end

        function Base.getproperty(t::$(class), s::Symbol)
            if s ∈ fieldnames($(class))
                return getfield(t, s)
            else
                return getproperty(getfield(t, :pyobj), s)
            end
        end

        export $(class)
    end
end
=#

for class in gate_classes
    @eval begin
        @pyclass gate $(class)
        export $(class)
    end
end

const gate_functions = [
    :add,
    :merge,
    :sqrtX,
    :sqrtXdag,
    :sqrtY,
    :sqrtYdag,
    :to_matrix_gate,
]

#=
for func in gate_functions
    @eval begin
        function $(func)(args...; kwargs...)
            qulacs.gate.$(func)(args...; kwargs...)
        end
        export $(func)
    end
end
=#

for func in gate_functions
    @eval begin
        @pyfunc gate $(func)
        export $(func)
    end
end

function __init__()
    #copy!(qulacs, pyimport("qulacs"))
    copy!(gate, pyimport("qulacs.gate"))
end

end # module Gate


text:
include("utils.jl")

@testset "Stiefel" begin
    @testset "Real" begin
        M = Stiefel(3, 2)
        @testset "Basics" begin
            @test repr(M) == "Stiefel(3, 2, ℝ)"
            x = [1.0 0.0; 0.0 1.0; 0.0 0.0]
            @test representation_size(M) == (3, 2)
            @test manifold_dimension(M) == 3
            base_manifold(M) === M
            @test_throws DomainError is_manifold_point(M, [1.0, 0.0, 0.0, 0.0], true)
            @test_throws DomainError is_manifold_point(
                M,
                1im * [1.0 0.0; 0.0 1.0; 0.0 0.0],
                true,
            )
            @test !is_tangent_vector(M, x, [0.0, 0.0, 1.0, 0.0])
            @test_throws DomainError is_tangent_vector(
                M,
                x,
                1 * im * zero_tangent_vector(M, x),
                true,
            )
        end
        @testset "Embedding and Projection" begin
            x = [1.0 0.0; 0.0 1.0; 0.0 0.0]
            y = similar(x)
            z = embed(M, x)
            @test z == x
            embed!(M, y, x)
            @test y == z
            a = [1.0 0.0; 0.0 2.0; 0.0 0.0]
            @test !is_manifold_point(M, a)
            b = similar(a)
            c = project(M, a)
            @test c == x
            project!(M, b, a)
            @test b == x
            X = [0.0 0.0; 0.0 0.0; -1.0 1.0]
            Y = similar(X)
            Z = embed(M, x, X)
            embed!(M, Y, x, X)
            @test Y == X
            @test Z == X
        end

        types = [Matrix{Float64}, ]
        TEST_FLOAT32 && push!(types, Matrix{Float32})
        TEST_STATIC_SIZED && push!(types, MMatrix{3, 2, Float64})

        @testset "Type $T" for T in types
            x = [1.0 0.0; 0.0 1.0; 0.0 0.0]
            y = exp(M, x, [0.0 0.0; 0.0 0.0; 1.0 1.0])
            z = exp(M, x, [0.0 0.0; 0.0 0.0; -1.0 1.0])
            pts = convert.(T, [x, y, z])
            v = inverse_retract(M, x, y, PolarInverseRetraction())
            @test !is_manifold_point(M, 2 * x)
            @test_throws DomainError !is_manifold_point(M, 2 * x, true)
            @test !is_tangent_vector(M, 2 * x, v)
            @test_throws DomainError !is_tangent_vector(M, 2 * x, v, true)
            @test !is_tangent_vector(M, x, y)
            @test_throws DomainError is_tangent_vector(M, x, y, true)
            test_manifold(
                M,
                pts,
                test_exp_log = false,
                default_inverse_retraction_method = PolarInverseRetraction(),
                test_injectivity_radius = false,
                test_is_tangent = true,
                test_project_tangent = true,
                test_vector_transport = false,
                test_forward_diff = false,
                test_reverse_diff = false,
                test_vee_hat = false,
                projection_atol_multiplier = 15.0,
                retraction_atol_multiplier = 10.0,
                is_tangent_atol_multiplier = 4 * 10.0^2,
                retraction_methods = [PolarRetraction(), QRRetraction()],
                inverse_retraction_methods = [
                    PolarInverseRetraction(),
                    QRInverseRetraction(),
                ],
            )

            @testset "inner/norm" begin
                v1 = inverse_retract(M, pts[1], pts[2], PolarInverseRetraction())
                v2 = inverse_retract(M, pts[1], pts[3], PolarInverseRetraction())

                @test real(inner(M, pts[1], v1, v2)) ≈ real(inner(M, pts[1], v2, v1))
                @test imag(inner(M, pts[1], v1, v2)) ≈ -imag(inner(M, pts[1], v2, v1))
                @test imag(inner(M, pts[1], v1, v1)) ≈ 0

                @test norm(M, pts[1], v1) isa Real
                @test norm(M, pts[1], v1) ≈ sqrt(inner(M, pts[1], v1, v1))
            end
        end
    end

    @testset "Complex" begin
        M = Stiefel(3, 2, ℂ)
        @testset "Basics" begin
            @test repr(M) == "Stiefel(3, 2, ℂ)"
            @test representation_size(M) == (3, 2)
            @test manifold_dimension(M) == 8
            @test Manifolds.allocation_promotion_function(M,exp!,(1,)) == complex
            @test !is_manifold_point(M, [1.0, 0.0, 0.0, 0.0])
            @test !is_tangent_vector(M, [1.0 0.0; 0.0 1.0; 0.0 0.0], [0.0, 0.0, 1.0, 0.0])
            x = [1.0 0.0; 0.0 1.0; 0.0 0.0]
        end
        types = [Matrix{ComplexF64}]
        @testset "Type $T" for T in types
            x = [0.5 + 0.5im 0.5 + 0.5im; 0.5 + 0.5im -0.5 - 0.5im; 0.0 0.0]
            y = exp(M, x, [0.0 0.0; 0.0 0.0; 1.0 1.0])
            z = exp(M, x, [0.0 0.0; 0.0 0.0; -1.0 1.0])
            pts = convert.(T, [x, y, z])
            v = inverse_retract(M, x, y, PolarInverseRetraction())
            @test !is_manifold_point(M, 2 * x)
            @test_throws DomainError !is_manifold_point(M, 2 * x, true)
            @test !is_tangent_vector(M, 2 * x, v)
            @test_throws DomainError !is_tangent_vector(M, 2 * x, v, true)
            @test !is_tangent_vector(M, x, y)
            @test_throws DomainError is_tangent_vector(M, x, y, true)
            test_manifold(
                M,
                pts,
                test_exp_log = false,
                default_inverse_retraction_method = PolarInverseRetraction(),
                test_injectivity_radius = false,
                test_is_tangent = true,
                test_project_tangent = true,
                test_vector_transport = false,
                test_forward_diff = false,
                test_reverse_diff = false,
                test_vee_hat = false,
                projection_atol_multiplier = 15.0,
                retraction_atol_multiplier = 10.0,
                is_tangent_atol_multiplier = 4 * 10.0^2,
                retraction_methods = [PolarRetraction(), QRRetraction()],
                inverse_retraction_methods = [
                    PolarInverseRetraction(),
                    QRInverseRetraction(),
                ],
            )

            @testset "inner/norm" begin
                v1 = inverse_retract(M, pts[1], pts[2], PolarInverseRetraction())
                v2 = inverse_retract(M, pts[1], pts[3], PolarInverseRetraction())

                @test real(inner(M, pts[1], v1, v2)) ≈ real(inner(M, pts[1], v2, v1))
                @test imag(inner(M, pts[1], v1, v2)) ≈ -imag(inner(M, pts[1], v2, v1))
                @test imag(inner(M, pts[1], v1, v1)) ≈ 0

                @test norm(M, pts[1], v1) isa Real
                @test norm(M, pts[1], v1) ≈ sqrt(inner(M, pts[1], v1, v1))
            end
        end
    end

    @testset "Quaternion" begin
        M = Stiefel(3, 2, ℍ)
        @testset "Basics" begin
            @test representation_size(M) == (3, 2)
            @test manifold_dimension(M) == 18
        end
    end
end


text:
using MultivariateFunctions
using Dates

tol = 10*eps()
const global_base_date = Date(2000,1,1)
StartDate = Date(2018, 7, 21)

x = StartDate .+ Dates.Day.(2 .* (1:1000 .- 1))

function ff(x::Date)
    days_between = years_from_global_base(x)
    return log(days_between) + sqrt(days_between)
end
y = ff.(x)

spline = create_linear_interpolation(x,y)
# Test if interpolating
all(abs.(evaluate.(Ref(spline), x) .- y) .< tol)

x_float = years_from_global_base.(x)
coefficient_in_first_interval = (y[2] - y[1])/(x_float[2] - x_float[1])
(coefficient_in_first_interval - spline.functions_[1].functions_[2].multiplier_) < tol
other_coefficient = (y[11] - y[10])/(x_float[11] - x_float[10])
(other_coefficient - spline.functions_[10].functions_[2].multiplier_) < tol

left_const = create_constant_interpolation_to_left(x,y)
right_const = create_constant_interpolation_to_right(x,y)
all(abs.(evaluate.(Ref(left_const), x_float .- tol) .- evaluate.(Ref(right_const), x_float)) .< tol)


text:
using Oceananigans.Operators

using Oceananigans.Operators: Δx_vᶜᶠᵃ, Δy_uᶠᶜᵃ, Δxᶠᶜᵃ, Δyᶜᶠᵃ, Az_wᶜᶜᵃ
using Oceananigans.Advection: div_Uu, div_Uv

######
###### Horizontally-vector-invariant formulation of momentum advection
######
###### Follows https://mitgcm.readthedocs.io/en/latest/algorithm/algorithm.html#vector-invariant-momentum-equations
######

@inline ϕ²(i, j, k, grid, ϕ) = @inbounds ϕ[i, j, k]^2
@inline Khᶜᶜᶜ(i, j, k, grid, u, v) = (ℑxᶜᵃᵃ(i, j, k, grid, ϕ², u) + ℑyᵃᶜᵃ(i, j, k, grid, ϕ², v)) / 2

@inbounds ζ₂wᶠᶜᶠ(i, j, k, grid, u, w) = ℑxᶠᵃᵃ(i, j, k, grid, Az_wᶜᶜᵃ, w) * δzᵃᵃᶠ(i, j, k, grid, u) / Δzᵃᵃᶠ(i, j, k, grid) / Azᶠᶜᵃ(i, j, k, grid)
@inbounds ζ₁wᶜᶠᶠ(i, j, k, grid, v, w) = ℑyᵃᶠᵃ(i, j, k, grid, Az_wᶜᶜᵃ, w) * δzᵃᵃᶠ(i, j, k, grid, v) / Δzᵃᵃᶠ(i, j, k, grid) / Azᶜᶠᵃ(i, j, k, grid)

@inline U_dot_∇u(i, j, k, grid, advection::VectorInvariant, U) = (
    - ℑyᵃᶜᵃ(i, j, k, grid, ζ₃ᶠᶠᵃ, U.u, U.v) * ℑxᶠᵃᵃ(i, j, k, grid, ℑyᵃᶜᵃ, Δx_vᶜᶠᵃ, U.v) / Δxᶠᶜᵃ(i, j, k, grid) # Vertical relative vorticity term
    + ∂xᶠᶜᵃ(i, j, k, grid, Khᶜᶜᶜ, U.u, U.v)    # Bernoulli head term
    + ℑzᵃᵃᶜ(i, j, k, grid, ζ₂wᶠᶜᶠ, U.u, U.w))  # Horizontal vorticity / vertical advection term

@inline U_dot_∇v(i, j, k, grid, advection::VectorInvariant, U) = (
    + ℑxᶜᵃᵃ(i, j, k, grid, ζ₃ᶠᶠᵃ, U.u, U.v) * ℑyᵃᶠᵃ(i, j, k, grid, ℑxᶜᵃᵃ, Δy_uᶠᶜᵃ, U.u) / Δyᶜᶠᵃ(i, j, k, grid) # Vertical relative vorticity term
    + ∂yᶜᶠᵃ(i, j, k, grid, Khᶜᶜᶜ, U.u, U.v)   # Bernoulli head term
    + ℑzᵃᵃᶜ(i, j, k, grid, ζ₁wᶜᶠᶠ, U.v, U.w)) # Horizontal vorticity / vertical advection term

######
###### Conservative formulation of momentum advection
######

@inline U_dot_∇u(i, j, k, grid, advection::AbstractAdvectionScheme, U) = div_Uu(i, j, k, grid, advection, U, U.u)
@inline U_dot_∇v(i, j, k, grid, advection::AbstractAdvectionScheme, U) = div_Uv(i, j, k, grid, advection, U, U.v)

######
###### No advection
######

@inline U_dot_∇u(i, j, k, grid::AbstractGrid{FT}, advection::Nothing, U) where FT = zero(FT)
@inline U_dot_∇v(i, j, k, grid::AbstractGrid{FT}, advection::Nothing, U) where FT = zero(FT)


text:
# Nonfunctional tests aren't testing whether it's correct but
# properties of the runs. It's a terrible, but amusing, name.
# Here we create combinatorial interaction test data using
# several algorithms and ask how long it takes to run and
# how many test cases it generates under different conditions.
using CSV
using Dates
using UnitTestDesign
using DataFrames
import Random: MersenneTwister
import StatsBase: trim, trimvar, winsor, mean

function ipog_extra(arity, n_way, M = nothing, rng = nothing)
    size(UnitTestDesign.ipog(arity, n_way), 2)
end

function greedy_extra(arity, n_way, M, rng)
    length(UnitTestDesign.n_way_coverage(arity, n_way, M, rng))
end

fut_map = Dict(:ipog => ipog_extra, :greedy => greedy_extra)

function single_test!(test_case, rng)
    arity = zeros(test_case.k, test_case.n)
    fill!(arity, test_case.r)
    callee = fut_map[test_case.fut]
    @timed callee(arity, test_case.wayness, test_case.M, rng)
end

function test_against_cases!(df, sample_cnt, rng)
    for idx in 1:size(df, 1)
        duration = zeros(sample_cnt)
        memory = zeros(sample_cnt)
        case_cnt = zeros(Int, sample_cnt)
        for sample_idx in 1:sample_cnt
            v = single_test!(df[idx, :], rng)
            duration[sample_idx] = v.time
            memory[sample_idx] = v.bytes / 1024.0^2
            case_cnt[sample_idx] = v.value
        end
        df[idx, :duration] = mean(trim(duration, prop = 0.2))
        df[idx, :duration_max] = maximum(duration)
        df[idx, :mbytes] = mean(trim(memory, prop = 0.2))
        df[idx, :mbytes_max] = maximum(memory)
        df[idx, :cases] = minimum(case_cnt)
        df[idx, :cases_max] = maximum(case_cnt)
    end
end

# n is the number of parameters
# r is the number of values for each parameter.
# wayness is 2-way or n-way coverage
# k is the data type to use to store the coverage
# fut is the function-under-test.
# M is a parameter for the greedy algorithm.
cases = DataFrame[]

# These tests come from comparison with
# Calvagna, A., and A. Gargantini. 2009. “IPO-S: Incremental Generation of
# Combinatorial Interaction Test Data Based on Symmetries of Covering Arrays.”
# In 2009 International Conference on Software Testing, Verification, and
# Validation Workshops, 10–18.
N = 10
calvagna_table_3 = DataFrame(
    n = [10i for i in 1:N],
    r = fill(4, N),
    wayness = fill(2, N),
    cases = [31, 34, 41, 42, 48, 48, 51, 51, 51, 53]
)
N = 6
calvagna_table4 = DataFrame(
    n = fill(10, N),
    r = [5i for i in 1:N],
    wayness = fill(2, N),
    cases = [47, 169, 361, 618, 956, 1355]
)
N = 5
# calvagna_table4 = DataFrame(
#     n = fill(3, N),
#     r = [3, 4, 5, 6, 7],
#     wayness = fill(2, N),
#     cases = [9, 9, 9, 9, 10, 9]
# )

N = 10
push!(cases, DataFrame(
    n = [10i for i in 1:N],
    r = fill(4, N),
    wayness = fill(2, N),
    k = fill(Int, N),
    fut = fill(:ipog, N),
    M = fill(0, N)
))
N = 3
push!(cases, DataFrame(
    n = [10i for i in 1:N],
    r = fill(4, N),
    wayness = fill(2, N),
    k = fill(Int, N),
    fut = fill(:greedy, N),
    M = fill(50, N)
))
push!(cases, DataFrame(
    n = [10i for i in 1:N],
    r = fill(4, N),
    wayness = fill(2, N),
    k = fill(Int, N),
    fut = fill(:greedy, N),
    M = fill(10, N)
))
push!(cases, DataFrame(
    n = [10i for i in 1:N],
    r = fill(4, N),
    wayness = fill(2, N),
    k = fill(Int8, N),  # Using Int8 for reduced memory consumption.
    fut = fill(:greedy, N),
    M = fill(10, N)
))
# Add tests from Lei, Yu, Raghu Kacker, D. Richard Kuhn, Vadim Okun,
# and James Lawrence. 2008. “IPOG/IPOG-D: Efficient Test Generation
# for Multi-Way Combinatorial Testing.” Software Testing, Verification
# & Reliability 18 (3): 125–48.
N = 4
lei_table_1 = DataFrame(
    n = fill(15, N),
    r = fill(4, N),
    wayness = [3, 4, 5, 6],
    cases = [181, 924, 4519, 20384],
    time = [0.56, 16.57, 230, 2152]
)
N = 10
lei_table_2 = DataFrame(
    n = 11:20,
    r = fill(4, N),
    wayness = fill(5, N),
    cases = [3287, 3703, 4001, 4260, 4519, 4787, 5018, 5245, 5471, 5685],
    time = [23.3, 44, 80, 139, 230, 368, 565, 839, 1206, 1739]
)
N = 6
lei_table_3 = DataFrame(
    n = fill(15, N),
    r = 2:7,
    wayness = fill(5, N),
    cases = [134, 1123, 4531, 15095, 37748, 81814],
    time = [4.08, 48, 234, 997, 3273, 9040]
)

# Test 1: 15 4-value parameters
N = 2
push!(cases, DataFrame(
    n = fill(15, N),
    r = fill(4, N),
    wayness = 2:(2 + N - 1),
    k = fill(Int, N),
    fut = fill(:ipog, N),
    M = fill(0, N)
))

df = vcat(cases...)
# Add a place in the dataframe to store results.
df[:, :duration] = zeros(N),
df[:, :duration_max] = zeros(N),
df[:, :mbytes] = zeros(Float64, N),
df[:, :mbytes_max] = zeros(Float64, N),
df[:, :cases] = zeros(Int, N),
df[:, :cases_max] = zeros(Int, N),

rng = MersenneTwister(947234)
test_against_cases!(df, 10, rng)
df

CSV.write("nonfunctional_cases_$(today()).csv", df)

match_options = [
    UnitTestDesign.case_compatible_with_tuple,
    UnitTestDesign.case_partial_cover,
    UnitTestDesign.case_covers_tuple
]
opt_cnt = length(match_options)^3
opt_res = zeros(Int, 4, opt_cnt)
for opt_choice in 1:opt_cnt
    strat_idx = digits(opt_choice - 1, base = 3, pad = 3) .+ 1
    if strat_idx[1] != 3
        strategy = Dict((a, match_options[b]) for (a, b) in
            zip([:lastparam, :missingvals, :expand], strat_idx))

        ret = UnitTestDesign.ipog_instrumented(fill(4, 20), 2, strategy)
        opt_res[:, opt_choice] = vcat(size(ret[1], 2), strat_idx)
    end
end
opt_res

match_options = [
    UnitTestDesign.case_compatible_with_tuple,
    UnitTestDesign.case_partial_cover,
    UnitTestDesign.case_covers_tuple
]
strat_cnt = 2
opt_cnt = length(match_options)^strat_cnt
opt_res = zeros(Int, 1 + strat_cnt, opt_cnt)
for opt_choice in 1:opt_cnt
    strat_idx = digits(opt_choice - 1, base = 3, pad = strat_cnt) .+ 1
    if strat_idx[1] != 3
        strategy = Dict((a, match_options[b]) for (a, b) in
            zip([:lastparam, :expand], strat_idx))

        ret = UnitTestDesign.ipog_bytuple_instrumented(fill(4, 20), 2, strategy)
        opt_res[:, opt_choice] = vcat(size(ret[1], 2), strat_idx)
    end
end
opt_res

strat_idx = [2, 1]
strategy = Dict((a, match_options[b]) for (a, b) in
    zip([:lastparam, :expand], strat_idx))
ret = UnitTestDesign.ipog_bytuple_instrumented(fill(4, 30), 2, strategy)


text:
"""
    ClimaCore.Topologies

Objects describing the horizontal connections between elements.

All elements are quadrilaterals, using the face and vertex numbering
convention from [p4est](https://p4est.github.io/papers/BursteddeWilcoxGhattas11.pdf):
```
          4
      3-------4
 ^    |       |
 |  1 |       | 2
x2    |       |
      1-------2
          3
        x1-->
```
"""
module Topologies

import ..Geometry
import ..Domains: Domains, coordinate_type
import ..Meshes:
    Meshes, EquispacedRectangleMesh, Mesh2D, AbstractMesh, TensorProductMesh

# TODO: seperate types for MPI/non-MPI topologies
"""
   AbstractTopology

Subtypes of `AbstractHorizontalTopology` define connectiveness of a
mesh in the horizontal domain.
"""
abstract type AbstractTopology end

"""
    domain(topology)

The `domain` underlying the topology.
"""
function domain end

coordinate_type(topology::AbstractTopology) = coordinate_type(domain(topology))

"""
    nlocalelems(topology)

The number of local elements in `topology`.
"""
function nlocalelems end

"""
    (c1,c2,c3,c4) = vertex_coordinates(topology, elem)

The coordinates of the 4 vertices of element `elem`.
"""
function vertex_coordinates end

"""
    (opelem, opface, reversed) = opposing_face(topology, elem, face)

The opposing face of face number `face` of element `elem` in `topology`.

- `opelem` is the opposing element number, 0 for a boundary, negative for a ghost element
- `opface` is the opposite face number, or boundary face number if a boundary
- `reversed` indicates whether the opposing face has the opposite orientation.
"""
function opposing_face end

"""
    i,j = face_node_index(face, Nq, q, reversed=false)

The node indices of the `q`th node on face `face`, where `Nq` is the number of
face nodes in each direction.
"""
function face_node_index(face, Nq, q, reversed = false)
    if reversed
        q = Nq - q + 1
    end
    if face == 1
        return 1, q
    elseif face == 2
        return Nq, q
    elseif face == 3
        return q, 1
    else
        return q, Nq
    end
end

"""
    i,j = vertex_node_index(vertex_num, Nq)

The node indices of `vertex_num`, where `Nq` is the number of face nodes in
each direction.
"""
function vertex_node_index(vertex_num, Nq)
    if vertex_num == 1
        return 1, 1
    elseif vertex_num == 2
        return Nq, 1
    elseif vertex_num == 3
        return 1, Nq
    else
        return Nq, Nq
    end
end

"""
    interior_faces(topology::AbstractTopology)

An iterator over the interior faces of `topology`. Each element of the iterator
is a 5-tuple the form

    (elem1, face1, elem2, face2, reversed)

where `elemX, faceX` are the element and face numbers, and `reversed` indicates
whether they have opposing orientations.
"""
function interior_faces(topology)
    InteriorFaceIterator(topology)
end

struct InteriorFaceIterator{T <: AbstractTopology}
    topology::T
end

"""
    boundaries(topology)

A `Tuple` or `NamedTuple` of the boundary tags of the topology. A boundary tag
is an integer that uniquely identifies a boundary.
"""
function boundaries end

"""
    boundary_faces(topology, boundarytag)

An iterator over the faces of `topology` which face the boundary with tag
`boundarytag`. Each element of the iterator is an `(elem, face)` pair.
"""
function boundary_faces(topology, boundarytag::Integer)
    BoundaryFaceIterator(topology, boundarytag)
end

struct BoundaryFaceIterator{T}
    topology::T
    boundary::Int
end

"""
    vertices(topology)

An iterator over the unique (shared) vertices of the topology `topology`.
Each vertex is an iterator over `(element, vertex_number)` pairs.
"""
function vertices(topology)
    VertexIterator(topology)
end
struct VertexIterator{T <: AbstractTopology}
    topology::T
end
struct Vertex{T <: AbstractTopology, V}
    topology::T
    num::V
end

# implementations
include("grids.jl")
include("grid2d.jl")

end # module


text:
pad(io, n)::Void = print(io, " "^n)

# function Base.show(io::IO, q::QueryNode, leftmargin=0)::Void
#     println(io, name(typeof(q)))
#     pad(io, leftmargin+2); println(io, "arguments:")
#     for (i, arg) in enumerate(q.args)
#         pad(io, leftmargin+6); println(io, "$i)  ", arg)
#     end
#     pad(io, leftmargin+2); println(io, "inputs:")
#     pad(io, leftmargin+6); print(io, "1)  ")
#     show(io, q.input, leftmargin+10)
#     return
# end

function Base.show(io::IO, q::QueryNode, leftmargin=0)::Void
    println(io, name(typeof(q)))
    pad(io, leftmargin+2); println(io, "arguments:")
    for (i, arg) in enumerate(q.args)
        pad(io, leftmargin+6); println(io, "$i)  ", arg)
    end
    pad(io, leftmargin+2); println(io, "inputs:")
    for (i, input) in enumerate(q.inputs)
        pad(io, leftmargin+6); print(io, "$i)  ")
        show(io, input, leftmargin+10)
    end
    return
end

function Base.show(io::IO, d::DataNode, leftmargin)::Void
    # println("Data source: ", d.input)
    println(io, name(typeof(d)))
    pad(io, leftmargin+2); print(io, "source:  ")
    if isdefined(d, :input)
        println(io, "source of type ", name(typeof(d.input)))
    else
        println(io, "unset source")
    end
    return
end


text:
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Date: December 2020

mutable struct SPG_params
    verbose::Integer
    optTol::Real
    progTol::Real
    maxIter::Integer
    suffDec::Real
    memory::Integer
    useSpectral::Bool
    curvilinear::Bool
    feasibleInit::Bool
    testOpt::Bool
    bbType::Bool
    testInit::Bool
    store_trace::Bool
    optNorm::Union{Real, Integer}
    iniStep::Real
    maxLinesearchIter::Integer
end
"""
    spg_options(;verbose=3,optTol=1f-5,progTol=1f-7,
                maxIter=20,suffDec=1f-4,memory=2,
                useSpectral=true,curvilinear=false,
                feasibleInit=false,testOpt=true,
                bbType=true,testInit=false, store_trace=false,
                optNorm=Inf,iniStep=1f0, maxLinesearchIter=10)

Options structure for Spectral Project Gradient algorithm.

# Arguments

- `verbose`: level of verbosity (0: no output, 1: iter (default))
- `optTol`: tolerance used to check for optimality (default: 1e-5)
- `progTol`: tolerance used to check for lack of progress (default: 1e-9)
- `maxIter`: maximum number of iterations (default: 20)
- `suffDec`: sufficient decrease parameter in Armijo condition (default: 1e-4)
- `memory`: number of steps to look back in non-monotone Armijo condition
- `useSpectral`: use spectral scaling of gradient direction (default: 1)
- `curvilinear`: backtrack along projection Arc (default: 0)
- `testOpt`: test optimality condition (default: 1)
- `feasibleInit`: if 1, then the initial point is assumed to be feasible
- `bbType`: type of Barzilai Borwein step (default: 1)
- `testInit`: Whether to test the initial estimate for optimality (default: false)
- `store_trace`: Whether to store the trace/history of x (default: false)
- `optNorm`: First-Order Optimality Conditions norm (default: Inf)
- `iniStep`: Initial step length estimate (default: 1)
- `maxLinesearchIter`: Maximum number of line search iteration (default: 20)
"""
function spg_options(;verbose=1,optTol=1f-10,progTol=1f-10,
                     maxIter=20,suffDec=1f-8,memory=2,
                     useSpectral=true,curvilinear=false,
                     feasibleInit=false,testOpt=true,
                     bbType=true,testInit=false, store_trace=false,
                     optNorm=Inf,iniStep=1, maxLinesearchIter=20)
    return SPG_params(verbose,optTol,progTol,
                      Int64(maxIter),suffDec,memory,
                      useSpectral,curvilinear,
                      feasibleInit,testOpt, bbType,testInit, store_trace,
                      optNorm,iniStep, Int64(maxLinesearchIter))
end

"""
    spg(funObj, x, funProj, options)

Function for using Spectral Projected Gradient to solve problems of the form
  min funObj(x) s.t. x in C

# Arguments

- `funObj(x)`:function to minimize (returns gradient as second argument)
- `funProj(x)`: function that returns projection of x onto C
- `x`: Initial guess
- `options`: spg_options structure

# Notes:

- if the projection is expensive to compute, you can reduce the
  number of projections by setting testOpt to 0 in the options

- Adapted fromt he matlab implementation of minConf_SPG
"""
function spg(funObj::Function, x::AbstractArray{T}, funProj::Function,
             options::SPG_params, ls=nothing) where {T}
    # Result structure
    sol = result(x)
    # Initialize array for gradient
    G = similar(x)

    # Setup Function to track number of evaluations
    projection(x) = (sol.n_project +=1; return funProj(x))
    grad!(g, x) = (sol.n_ϕeval +=1; sol.n_geval +=1 ; g .= funObj(x)[2])
    objgrad!(g, x) = (sol.n_ϕeval +=1; sol.n_geval +=1 ;(obj, g0) = funObj(x); g .= g0; return obj)
    obj(x) = objgrad!(G, x)

    # Solve optimization
    return _spg(obj, grad!, objgrad!, projection, x, G, sol, ls, options)
end

"""
    spg(f, g!, fg!, x, funProj, options)

Function for using Spectral Projected Gradient to solve problems of the form
min funObj(x) s.t. x in C

# Arguments
- `f(x)`: function to minimize (returns objective only)
- `g!(g, x)`: gradient of function (in place)
- `fg!(g, x)`: objective and gradient (in place)
- `funProj(x)`: function that returns projection of x onto C
- `x`: Initial guess
- `options`: spg_options structure

# Notes:
- if the projection is expensive to compute, you can reduce the
      number of projections by setting testOpt to 0 in the options

- Adapted fromt he matlab implementation of minConf_SPG
"""
function spg(f::Function, g!::Function, fg!::Function, x::AbstractArray{T},
             funProj::Function, options::SPG_params, ls=nothing) where {T}
    # Result structure
    sol = result(x)
    # Initialize array for gradient
    G = similar(x)

    # Setup Function to track number of evaluations
    projection(x) = (sol.n_project +=1; return funProj(x))
    obj(x) = (sol.n_ϕeval +=1 ; return f(x))
    grad!(g, x) = (sol.n_geval +=1 ; return g!(g, x))
    objgrad!(g, x) = (sol.n_ϕeval +=1;sol.n_geval +=1 ; return fg!(g, x))

    # Solve optimization
    return _spg(obj, grad!, objgrad!, projection, x, G, sol, ls, options)
end

"""
Low level SPG solver
"""
function _spg(obj::Function, grad!::Function, objgrad!::Function, projection::Function,
              x::AbstractArray{T}, g::AbstractArray{T}, sol::result, ls, options::SPG_params) where {T}
    # Initialize local variables
    nVars = length(x)
    old_ϕvals = -T(Inf)*ones(T, options.memory)
    d = similar(x)
    optCond = 0
    # Best solution
    x_best = x

    # Line search function
    isnothing(ls) && (ls = BackTracking{T}(order=3, iterations=options.maxLinesearchIter))
    checkls(ls)

    # Evaluate Initial Point and objective function
    ~options.feasibleInit && (x = projection(x))
    ϕ = objgrad!(g, x)
    ϕ_best = ϕ
    update!(sol; iter=1, ϕ=ϕ, g=g, x=x, store_trace=options.store_trace)

    # Output Log
    options.testOpt ? optCond = norm(projection(x-g)-x, options.optNorm) : optCond = 0
    init_log(ϕ, norm(projection(x-g)-x, options.optNorm), options, ls)

    # Optionally check optimality
    if options.testOpt && options.testInit
        if optCond < optTol
            verbose > 0 &&  @printf("First-Order Optimality Conditions Below optTol at Initial Point, norm g is %5.4f \n", norm(g))
            return
        end
    end

    # Start iterations
    for i = 1:options.maxIter
        # Compute Step Directional
        if i == 1 || ~options.useSpectral
            alpha = T(.1*norm(x, Inf)/norm(g, Inf))
        else
            y = g - sol.g
            s = x - sol.x
            options.bbType ? alpha = dot(s,s)/dot(s,y) : alpha = dot(s,y)/dot(y,y)
        end
        # Make sure alpha value is valid
        (alpha <= 1e-10 || alpha > 1e10 || ~isLegal(alpha)) && (alpha = T(1))
        # Update log
        i>1 && update!(sol; iter=i, ϕ=ϕ, x=x, g=g, store_trace=options.store_trace)
        # Compute Step
        @. d = -T(alpha).*g

        # Compute Projected Step
        ~options.curvilinear && (d .= projection(x + d) - x)

        # Check that Progress can be made along the direction
        gtd = dot(g,d)

        # Select Initial Guess to step length
        t = T(options.iniStep)

        # Compute reference function for non-monotone condition
        old_ϕvals[i%options.memory + 1] = T(ϕ)

        # Line search
        t, ϕ = linesearch(ls, sol, d, obj, grad!, objgrad!, t, maximum(old_ϕvals), gtd, g)
        x .= projection(sol.x + t*d)
        g == sol.g && grad!(g, x)

        # Check conditioning
        options.testOpt ? optCond = norm(projection(x-g)-x, options.optNorm) : optCond = Inf

        # Check if terminate
        i>1 && (terminate(options, optCond, t, d, ϕ, sol.ϕ) && break)

        # Check if better than best solution
        ϕ < ϕ_best && (x_best = x; ϕ_best = ϕ)

        # Output Log
        iter_log(i, sol, t, alpha, ϕ, optCond, options)
    end

    # Restore best iteration
    update!(sol; iter=options.maxIter+1, ϕ=ϕ_best, x=x_best, g=g, store_trace=options.store_trace)
    return sol
end

"""
Loging utilities
"""
function init_log(ϕ, optCond, options, ls)
    if options.verbose > 0
        @printf("Running SPG...\n");
        @printf("Number of objective function to store: %d\n",options.memory);
        @printf("Using  spectral projection : %s\n",options.useSpectral);
        @printf("Maximum number of iterations: %d\n",options.maxIter);
        @printf("SPG optimality tolerance: %.2e\n",options.optTol);
        @printf("SPG progress tolerance: %.2e\n",options.progTol);
        @printf("Line search: %s\n", typeof(ls))
     end
    if options.verbose > 0
        if options.testOpt
            @printf("%10s %10s %10s %10s %15s %15s %15s %15s\n","Iteration","FunEvals", "GradEvals","Projections","Step Length","alpha","Function Val","Opt Cond")
            @printf("%10d %10d %10d %10d %15.5e %15.5e %15.5e %15.5e\n",0,0,0,0,0,0,ϕ,optCond)
        else
            @printf("%10s %10s %10s %10s %15s %15s %15s\n","Iteration","FunEvals","GradEvals","Projections","Step Length","alpha","Function Val")
            @printf("%10d %10d %10d %10d %15.5e %15.5e %15.5e\n",0,0,0,0,0,0,ϕ)
        end
    end
end

function iter_log(i, sol, t, alpha, ϕ, optCond, options)
    if options.verbose > 0
        if options.testOpt
            @printf("%10d %10d %10d %10d %15.5e %15.5e %15.5e %15.5e\n",i,sol.n_ϕeval,sol.n_geval,sol.n_project,t,alpha,ϕ,optCond)
        else
            @printf("%10d %10d %10d %10d %15.5e %15.5e %15.5e\n",i,sol.n_ϕeval,sol.n_geval,sol.n_project,t,alpha,ϕ)
        end
    end
end


text:
module Viewshed_Test

profile1 = [
    (0,12),
    (25,12),
    (50,8),
    (75,11),
    (100,13),
    (125,11),
    (150,15),
    (175,14),
    (200,23)
]

profile2 = [
    (0,500),
    (25,12),
    (50,8),
    (75,10),
    (100,300),
    (125,11),
    (150,15),
    (175,14),
    (200,150)
]

profile3 = [
    (0,12),
    (25,12),
    (50,12),
    (75,12),
    (100,12),
    (125,12),
    (150,8),
    (175,10),
    (200,23)
]

veiwshed.([ profile1, profile2, profile3 ])

end # module

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
function build_sysimg(
    m::Module,
    options::Options.Comonicon;
    incremental = options.sysimg.incremental,
    cpu_target = options.sysimg.cpu_target,
    filter_stdlibs = options.sysimg.filter_stdlibs,
)

    precompile_execution_file = String[pkgdir(m, x) for x in options.sysimg.precompile.execution_file]
    precompile_statements_file =
        String[pkgdir(m, x) for x in options.sysimg.precompile.statements_file]

    create_sysimage(
        [nameof(m)];
        sysimage_path = sysimg_dylib(m, options),
        project = create_command_env(m),
        incremental,
        cpu_target,
        filter_stdlibs,
        precompile_execution_file,
        precompile_statements_file,
    )
    return
end

"""
    create_command_env(m::Module, envpath::String=mktempdir(); test_deps=true)

Create an environment to execute the CLI command.

# Arguments

- `m`: the CLI module.
- `envpath`: the generated environment path.

# Keyword Arguments

- `test_deps`: include test deps or not.
"""
function create_command_env(m::Module, envpath::String = mktempdir(); test_deps::Bool = true)
    project = Pkg.Types.projectfile_path(pkgdir(m))
    ctx = Pkg.Types.Context(env = Pkg.Types.EnvCache(project))
    cmd_project = Pkg.Types.Project()
    merge!(cmd_project.deps, ctx.env.project.deps)
    merge!(cmd_project.compat, ctx.env.project.compat)
    # add the package to dependencies
    cmd_project.deps[ctx.env.project.name] = ctx.env.project.uuid

    if test_deps
        merge!(cmd_project.deps, ctx.env.project.extras)
        test_project = Pkg.Types.projectfile_path(pkgdir(m, "test"))
        if !isnothing(test_project) # add test deps if has it
            test_env = Pkg.Types.EnvCache(test_project)
            merge!(cmd_project.deps, test_env.project.deps)
            merge!(cmd_project.compat, test_env.project.compat)
        end
    end

    cmd_project_path = Pkg.Types.projectfile_path(envpath)
    cmd_manifest_path = Pkg.Types.manifestfile_path(envpath)
    Pkg.Types.write_project(cmd_project, cmd_project_path)

    # update Manifest from project
    if v"1.6" ≤ VERSION < v"1.7-"
        pkg_manifest = Pkg.Operations.abspath!(ctx, ctx.env.manifest)
    elseif v"1.7" ≤ VERSION < v"1.8-"
        pkg_manifest = Pkg.Operations.abspath!(ctx.env, ctx.env.manifest)
    else
        error("unsupported Julia version: $VERSION")
    end
    # TODO: merge test into the package manifest
    # how does TestENV do it? it's a bit unclear
    # cmd_env = Pkg.Types.EnvCache(cmd_project_path)
    # cmd_manifest = cmd_env.manifest
    # for (name, uuid) in ctx.env.project.deps
    #     entry = get(pkg_manifest, uuid, nothing)
    #     if !isnothing(entry) && Pkg.Operations.isfixed(entry)
    #         subgraph = Pkg.Operations.prune_manifest(pkg_manifest, [uuid])
    #         for (uuid, entry) in subgraph
    #             if haskey(cmd_manifest, uuid)
    #                 @show cmd_manifest[uuid]
    #                 Pkg.Types.pkgerror("can not merge projects")
    #             end
    #             cmd_manifest[uuid] = entry
    #         end
    #     end
    # end

    Pkg.Types.write_manifest(pkg_manifest, cmd_manifest_path)
    current_project = Base.current_project()
    Pkg.activate(envpath)
    Pkg.develop(Pkg.PackageSpec(path = pkgdir(m)))
    Pkg.update()
    Pkg.activate(current_project)
    return envpath
end

function build_sysimg_tarball(m::Module, options::Options.Comonicon)
    dylib = sysimg_dylib(m, options)
    tarball = tarball_name(m, options.name, "sysimg")

    tmp_dir = mktempdir()
    cd(tmp_dir) do
        mkdir("sysimg")
        cp(dylib, joinpath("sysimg", basename(dylib)))
        run(`tar -czvf $tarball sysimg`)
    end
    mv(joinpath(tmp_dir, tarball), joinpath(pwd(), tarball))
    return
end

function download_sysimg(m::Module, options::Options.Comonicon)
    url = sysimg_url(m, options)
    isnothing(url) && error("$m does not have a host")
    tarball = download(url)
    dylib = sysimg_dylib(m, options)
    Pkg.unpack(tarball, dylib)
    # NOTE: sysimg won't be shared, so we can just remove it
    isfile(tarball) && rm(tarball)
    return
end

function sysimg_url(m::Module, options::Options.Comonicon)
    isnothing(options.download) && return
    host = options.download.host
    if host == "github.com"
        url =
            "https://github.com/" *
            options.download.user *
            "/" *
            options.download.repo *
            "/releases/download/"
    else
        error("host $host is not supported, please open an issue at $COMONICON_URL")
    end

    tarball = tarball_name(mod, name, "sysimg")
    url *= "v$(get_version(mod))/$tarball"
    return url
end

function sysimg_dylib(m::Module, options::Options.Comonicon)
    dylib_name = "lib$(options.name).$(Libdl.dlext)"
    sysimg_dir = get_scratch!(m, "sysimg")
    return joinpath(sysimg_dir, dylib_name)
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, FerriteViz, WGLMakie

makedocs(sitename="FerriteViz",
         modules=[FerriteViz],
         authors="Maximilian Köhler",
         pages=["Home"=> "index.md",
                "Tutorial" => "tutorial.md",
                "Advanced Topics" => "atopics.md",
                "API Reference" => "api.md",]
)

deploydocs(repo = "github.com/Ferrite-FEM/FerriteViz.jl.git",
           push_preview=true,
           forcepush=true,)


text:
"""
	brlength(p, q)

Get the length of the branch between the nodes `p` and `q`.
"""
brlength(p::AbstractNode, q::AbstractNode) = getbranch(p, q).length

"""
    brlength!(p, q, length)

Set the `length` (`Float64` or `Nothing`) of the branch between the nodes `p` and `q`.
"""
function brlength!(p::AbstractNode, q::AbstractNode, v::Union{Nothing, Float64})
    getbranch(p, q).length = v

    return nothing
end

"""
brlabel(p, q)

Get the label of the branch between nodes `p` and `q`.
"""
brlabel(p::AbstractNode, q::AbstractNode) = getbranch(p, q).label

"""
brlabel!(p, q, v)

Set the string `v` as label of branch between the nodes `p` and `q`.
"""
function brlabel!(p::AbstractNode, q::AbstractNode, v::String)
	getbranch(p, q).label = v

	return nothing
end

"""
	brannotation(p, q, k)

Get the value of the annotation with key `k` in the branch between nodes `p` and `q`.
"""
brannotation(p::AbstractNode, q::AbstractNode, k::AbstractString) =
	getbranch(p, q).annotations[k]

"""
	brannotate!(node, key, value)

Set the `value` for an annotation `key` in the branch between the nodes `p` and `q`.

Annotation keys can only be `AbstractString`s.
"""
function brannotate!(p::AbstractNode, q::AbstractNode, k::AbstractString, v::Any)
	getbranch(p, q).annotations[k] = v

	return nothing
end

"""
	brdeannotate!(p, q, key)

Remove an annotation `key` from the branch between the nodes `p` and `q`.
"""
function brdeannotate!(p::AbstractNode, q::AbstractNode, k::AbstractString)
	delete!(getbranch(p, q).annotations, k)

	return nothing
end


text:
using DataDeps

register(DataDep(
  "planning-relax",
  "http://archive.ics.uci.edu/ml/datasets/Planning+Relax",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/00230/plrx.txt",
  "88d6db30747ae5090525ecbb512eb32f552eafb9e78b180ca0dfa37667e93c74",
  post_fetch_method=(path -> begin
    # First char of each row is whitespace, remove
    write(path, join(strip.(readlines(path)), "\n"))

    UCIData.process_dataset(path,
        target_index=13,
        feature_indices=1:12,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
function tarball_name(m::Module, name::String, postfix = "sysimg")
    return "$name-$postfix-$(get_version(m))-julia-$VERSION-$(osname())-$(Sys.ARCH).tar.gz"
end

"""
    osname()

Return the name of OS, will be used in building tarball.
"""
function osname()
    return Sys.isapple() ? "darwin" :
           Sys.islinux() ? "linux" :
           error("unsupported OS, please open an issue to request support at $COMONICON_URL")
end

function pack_tarball(tarball::String, root::String, path::String)
    cd(root) do
        run(`tar -czvf $tarball $path`)
    end
end


text:
export Fractal, Triangle, Carpet
"""
A `Fractal` contains an `IFS` (iterated function system) and a
`PixelPicture` (bitmap image) that the `IFS` generates.

Standard constructors:
* `Fractal(F::IFS,P::PixelPicture)`
* `Fractal(F::IFS,n::Int)`
"""
type Fractal
  F::IFS
  P::PixelPicture
  function Fractal(FF::IFS,PP::PixelPicture)
    new(FF,PP)
  end
end

function Fractal(F::IFS, n::Int)
  @assert n>0
  P = PixelPicture(n)
  return Fractal(F,n)
end

draw(F::Fractal) = draw(F.P)
size(F::Fractal) = size(F.P)

function show(io::IO,FF::Fractal)
  n = size(FF)
  println(io,"$n-by-$n fractal image with this iterated function system:")
  show(io,FF.F)
end

randomize!(FF::Fractal) = randomize!(FF.P)

export step!

"""
`step!(FR)` applies a fractal's IFS to its image returning `true`
if there has been a change and `false` if no change.
"""
function step!(FR::Fractal)
  n = size(FR.P)
  pts = get_points(FR.P)
  newP = PixelPicture(n)
  for xy in pts
    for vv in FR.F(xy)
      set_point!(newP,vv)
    end
  end
  if FR.P.data != newP.data
    FR.P.data = newP.data
    return true
  end
  return false
end

"""
`Triangle(n)` set up `n`-by-`n` Triangle triangle.
"""
function Triangle(n::Int)
  F = Triangle_IFS()
  P = PixelPicture(n)
  FF = Fractal(F,P)

  return FF
end

"""
`Carpet(n)` set up an `n`-by-`n` Carpet gasket.
"""
function Carpet(n::Int)
  F = Carpet_IFS()
  P = PixelPicture(n)
  FF = Fractal(F,P)

  return FF
end

set_point!(F::Fractal,xy::Array{Float64,1})=set_point!(F.P,xy)

get_points(F::Fractal)=get_points(F.P)


text:
@testset  "Test wrapper for greedysolvers" begin

    ψ = randn(100, 20)
    u = randn(100)
    m, n = size(ψ)
    ctrue = ψ\u

    idxls, cls, ϵls = lsomp(ψ, u; invert = true, verbose = true, ϵrel = eps())
    idxqr, cqr, ϵqr = qromp(ψ, u; invert = true, verbose = true, ϵrel = eps())
    idxpivot, cpivot, ϵpivot = pivotedqr(ψ, u; invert = true, verbose = true, ϵrel = eps())

    # Test LS OMP
    idxgreedy, cgreedy, ϵgreedy = greedysolver("lsomp", ψ, u; invert = true, verbose = true, ϵrel = eps())

    @test norm(idxls - idxgreedy)<1e-14
    @test norm(cls - cgreedy)<1e-14
    @test norm(ϵls - ϵgreedy)<1e-14

    # Test QR OMP
    idxgreedy, cgreedy, ϵgreedy = greedysolver("qromp", ψ, u; invert = true, verbose = true, ϵrel = eps())

    @test norm(idxqr - idxgreedy)<1e-14
    @test norm(cqr - cgreedy)<1e-14
    @test norm(ϵqr - ϵgreedy)<1e-14

    # Test Pivoted QR
    idxgreedy, cgreedy, ϵgreedy = greedysolver("pivotedqr", ψ, u; invert = true, verbose = true, ϵrel = eps())

    @test norm(idxpivot - idxgreedy)<1e-14
    @test norm(cpivot - cgreedy)<1e-14
    @test norm(ϵpivot - ϵgreedy)<1e-14
end

max_forks_repo_name: mleprovost/QROMP