text:
using CSV, DataFrames, DataFramesMeta

## Loading data
datadir="/Volumes/SSD Hans/Github/Research_ideas/Testing Market Power/data/Carvajal et al 2013/CDFQ_ECMA_Files/oildata.csv"

data=CSV.read(datadir, DataFrame)
data[:, Between(:algeria, :uk)]
data.id=1:size(data, 1)
ld=stack(data, Between(:algeria, :uk))
country=names(data, Between(:algeria, :uk))
@transform!(ld,q=:value ./ :totalworld)

df=unstack(ld, [:year,:variable], :month, :q, allowduplicates=true)

p=unstack(data[data.year.<2009,:], :year, :month, :oilpricereal2009)

# p array: 36 years x 12 months
p=Array(p[:,Not(:year)])

df[df.year.<2009, Not([:year, :variable])]

# q array, 12 months x 19 countries x 36 years
q=reshape(Array(ld[ld.year.<2009,:q]),12,19,36)

q[:,1,1].==@linq ld |>
    where(:year.==1973, :variable.=="algeria") |>
    select(:q)

@linq ld |>
        where(:year.==1973, :variable.=="algeria") |>
        select(:q)

q[1,1,:]

p=nothing
q=nothing

##
p=@linq data |>
    where(:year.<2009) |>
    select(:oilpricereal2009)

q = @linq data |>
        where(:year.<2009)

q = q[:,Between(:algeria, :uk)]

max_stars_repo_name: hans-mtz/SCR

max_issues_repo_name: hans-mtz/SCR

max_forks_repo_name: hans-mtz/SCR


text:
module IrisFloquet

end # module

ext: jl


text:
"""
    NewtonParameters

A composite type for the parameters of the parameters of simplified Newton in [`ImplicitRungeKuttaSolver`](@ref).

# Constructors
```julia
NewtonParameters(; ϵ = 1e-3, K = 10)
```

# Arguments
- `ϵ :: Real`    : relative tolerance
- `K :: Integer` : maximum number of iterations.

# Functions
- [`show`   ](@ref) : shows name and contents.
- [`summary`](@ref) : shows name.
"""
mutable struct NewtonParameters{ϵ_T, K_T}
    ϵ::ϵ_T
    K::K_T
end

NewtonParameters(; ϵ = 1e-3, K = 10) = NewtonParameters(ϵ, K)

############################################################################################
#                                         FUNCTIONS                                        #
############################################################################################

"""
    show(io::IO, newton::NewtonParameters)

prints a full description of `newton` and its contents to a stream `io`.
"""
Base.show(io::IO, newton::NewtonParameters) = NSDEBase._show(io, newton)

"""
    summary(io::IO, newton::NewtonParameters)

prints a brief description of `newton` to a stream `io`.
"""
Base.summary(io::IO, newton::NewtonParameters) = NSDEBase._summary(io, newton)


text:
export AbstractAgent, @agent, GraphAgent, GridAgent, ContinuousAgent, OSMAgent

"""
    YourAgentType <: AbstractAgent
Agents participating in Agents.jl simulations are instances of user-defined Types that
are subtypes of `AbstractAgent`. It is almost always the case that mutable Types make
for a simpler modelling experience.

Your agent type(s) **must have** the `id` field as first field.
Depending on the space structure there might be a `pos` field of appropriate type
and a `vel` field of appropriate type.
Each space structure quantifies precicely what extra fields (if any) are necessary,
however we recommend to use the [`@agent`](@ref) macro to help you create the agent type.

Your agent type may have other additional fields relevant to your system,
for example variable quantities like "status" or other "counters".

As an example, a [`GraphSpace`](@ref) requires an `id::Int` field and a `pos::Int` field.
To make an agent with two additional properties, `weight, happy`, we'd write
```julia
mutable struct ExampleAgent <: AbstractAgent
    id::Int
    pos::Int
    weight::Float64
    happy::Bool
end
```
"""
abstract type AbstractAgent end

"""
    @agent YourAgentType{X, Y} AgentSupertype begin
        some_property::X
        other_extra_property::Y
        # etc...
    end

Create a struct for your agents which includes the mandatory fields required to operate
in a particular space. Depending on the space of your model, the `AgentSupertype` is
chosen appropriately from [`GraphAgent`](@ref), [`GridAgent`](@ref),
[`ContinuousAgent`](@ref), [`OSMAgent`](@ref).

## Example
Using
```julia
@agent Person{T} GridAgent{2} begin
    age::Int
    moneyz::T
end
```
will in fact create an agent appropriate for using with 2-dimensional [`GridSpace`](@ref)
```julia
mutable struct Person{T} <: AbstractAgent
    id::Int
    pos::NTuple{2, Int}
    age::Int
    moneyz::T
end
```
"""
macro agent(name, base, fields)
    base_type = Core.eval(@__MODULE__, base)
    base_fieldnames = fieldnames(base_type)
    base_types = [t for t in base_type.types]
    base_fields = [:($f::$T) for (f, T) in zip(base_fieldnames, base_types)]
    res = :(mutable struct $(esc(name)) <: AbstractAgent end)
    push!(res.args[end].args, base_fields...)
    push!(res.args[end].args, map(esc,fields.args)...)
    return res
end

"""
    GraphAgent
Combine with [`@agent`](@ref) to create an agent type for [`GraphSpace`](@ref).
It attributes the fields `id::Int, pos::Int` to the start of the agent type.
"""
mutable struct GraphAgent <: AbstractAgent
    id::Int
    pos::Int
end

"""
    GridAgent{D}
Combine with [`@agent`](@ref) to create an agent type for `D`-dimensional
[`GridSpace`](@ref). It attributes the fields `id::Int, pos::NTuple{D,Int}`
to the start of the agent type.
"""
mutable struct GridAgent{D} <: AbstractAgent
    id::Int
    pos::Dims{D}
end

"""
    ContinuousAgent{D}
Combine with [`@agent`](@ref) to create an agent type for `D`-dimensional
[`ContinuousSpace`](@ref). It attributes the fields
`id::Int, pos::NTuple{D,Float64}, vel::NTuple{D,Float64}`
to the start of the agent type.
"""
mutable struct ContinuousAgent{D} <: AbstractAgent
    id::Int
    pos::NTuple{D,Float64}
    vel::NTuple{D,Float64}
end

"""
    OSMAgent
Combine with [`@agent`](@ref) to create an agent type for [`OpenStreetMapSpace`](@ref).
It attributes the fields
`id::Int, pos::Tuple{Int,Int,Float64}, route::Vector{Int}, destination::Tuple{Int,Int,Float64}`
to the start of the agent type.
"""
mutable struct OSMAgent <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int,Float64}
    route::Vector{Int}
    destination::Tuple{Int,Int,Float64}
end


text:
function readinput(file)
    # Length of the input file
    num = countlines(file)
    # Actually reading of the file and saving to an array
    input = readdlm(file)
    # Extracting of numerical values
    values = []
    arguments = ["length", "time", "delta_t", "gravity", "surface_tension"]
    for arg in arguments
        push!(values, findfirst(x -> x == arg, vec(input)))
    end
    length = input[values[1], 2]
    maxtime = input[values[2], 2]
    deltatime =  input[values[3], 2]
    γ = input[values[5], 2]
    g = input[values[4], 2]

    return length, maxtime, deltatime, g, γ
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

"""
Formulation described in:

    Carrión, M., & Arroyo, J. M. (2006). A computationally efficient
    mixed-integer linear formulation for the thermal unit commitment problem.
    IEEE Transactions on power systems, 21(3), 1371-1378.
    DOI: https://doi.org/10.1109/TPWRS.2006.876672
"""
module CarArr2006

import ..PiecewiseLinearCostsFormulation

struct PwlCosts <: PiecewiseLinearCostsFormulation end

end


text:
# SDE.jl

"""
    abstract type AbstactSDE <: Any

subtype: SDE.

A supertype for implementing SDEs.
"""
abstract type AbstractSDE end

"""
    struct SDE{S, T} <: AbstractSDE where {S<:Function, T<:SDEModel}
      b::S
      model::T
    end

# Constructors
```julia
SDE(b::S, model::T) where {S<:Function, T<:SDEModel}
SDE(model::T, b::S) where {S<:Function, T<:SDEModel}
```

Implements a stochastic differential equation dX_t=b(X_t)dt+σ(X_t)dW_t on time
interval [0,model.endtime], with W_t a Brownian motion. The drift function b is
a function on the real line. model.σ (either a function or a number) gives the diffusion.
The begin value X_0=model.beginvalue and the sample path is discretised via the
Euler-Maruyama scheme with time steps model.Δ.

# Example

```julia
model = SDEModel(x->2+sin(x), 0.0,10.0,0.001)
sde = SDE(sin, model)
```
"""
struct SDE{S, T} <: AbstractSDE where {S<:Function, T<:SDEModel}
  b::S
  model::T
end

# Constructor
SDE(model::SDEModel, b::Function) = SDE(b, model)

"""
    calculatenextsamplevalue(prevXval, σ::Float64, b, Δ, BMincrement)
    calculatenextsamplevalue(prevXval, σ, b, Δ, BMincrement)

Internal function, not exported.

Given the previous samplepath value prevXval, it calculates the next sample value with the Euler-Maruyama scheme.
"""
calculatenextsamplevalue(prevXval, σ::Float64, b, Δ, BMincrement) = prevXval + b(prevXval)*Δ + σ*BMincrement
calculatenextsamplevalue(prevXval, σ, b, Δ, BMincrement) = prevXval + b(prevXval)*Δ + σ(prevXval)*BMincrement

"""
  rand(sde::SDE)

Returns a sample path from the SDE sde, from time 0.0 to time sde.model.endtime,
discretised with precision sde.model.Δ with the help of the Euler-Maruyama scheme.

# Examples

```julia
model = SDEModel(1.0,0.0,10.0,0.01)
sde = SDE(sin, model)
X = rand(sde)
```
"""
function rand(sde::SDE)
    timeinterval = 0.0:sde.model.Δ:sde.model.endtime
    lengthoftimeinterval = length(timeinterval)
    BMincrements = sqrt(sde.model.Δ) .* randn(lengthoftimeinterval)
    samplevalues = Array{Float64}(undef, lengthoftimeinterval)
    prevXval = samplevalues[1] = sde.model.beginvalue
    for k in 2:lengthoftimeinterval
        prevXval = samplevalues[k] = calculatenextsamplevalue(prevXval, sde.model.σ, sde.b, sde.model.Δ, BMincrements[k])
    end
    return SamplePath(timeinterval, samplevalues)
end


text:
using TikzGraphs
using TikzPictures
using GraphRecipes
using LightGraphs
using Plots

default(size=(1000, 1000))

f = :(x * y + sin(x))

vertices = Dict([
1 => (name = "x", 			x = 1, y = 0, color = :lightblue, ),
2 => (name = "y", 			x = 0, y = 1, color = :lightblue, ),
3 => (name = "h₁ = sin", 	x = 1, y = 2, color = :lightblue, ),
4 => (name = "h₂ = *", 		x = 0, y = 3, color = :lightblue, ),
5 => (name = "h₃ = +", 		x = 1, y = 4, color = :lightblue, ),
6 => (name = "z", 			x = 1, y = 5, color = :lightblue, ),
7 => (name = "∂z/∂h₃", 		x = 2.5, y = 4, color = :orange, ),
8 => (name = "∂h₃/∂h₂", 	x = 2.5, y = 3, color = :orange, ),
9 => (name = "∂z/∂h₂", 		x = 3.5, y = 3, color = :orange, ),
10 => (name = "∂h₃/∂h₁", 	x = 2.5, y = 2, color = :orange, ),
11 => (name = "∂z/∂h₁", 	x = 3.5, y = 2, color = :orange, ),
12 => (name = "∂h₂/∂y", 	x = 4.5, y = 1, color = :orange, ),
13 => (name = "∂h₂/∂x + ∂h₁/∂x", x = 4.5, y = 0, color = :orange, ),
14 => (name = "∂z/∂x", x = 6, y = 0, color = :orange, ),
])

n = length(vertices)
g = LightGraphs.DiGraph(n)
add_edge!(g, 1, 3)  # x -> sin
add_edge!(g, 1, 4)	# x ->  *
add_edge!(g, 2, 4)	# y ->  *
add_edge!(g, 3, 5)	# sin(x) -> sin(x) + x*y
add_edge!(g, 4, 5)  # x*y -> sin(x) + x*y
add_edge!(g, 5, 6)  # sin(x) + x*y -> z

#add_edge!(g, 5, 7,)	# ∂z/∂h₃ -> z
#add_edge!(g, 4, 8,)	# ∂h₃/∂h₂ -> h₂ = *
#add_edge!(g, 8, 9,)	# ∂z/∂h₂ -> ∂h₃/∂h₂
#add_edge!(g, 7, 9,)	# ∂z/∂h₂ -> ∂h₃/∂h₂
#add_edge!(g, 3, 10,)	# ∂h₃/∂h₁ -> h₁
#add_edge!(g, 10, 11,)# ∂z/∂h₁  -> ∂h₃/∂h₁
#add_edge!(g, 7, 11,) # ∂z/∂h₁  -> ∂z/∂h₃
#add_edge!(g, 2, 12,) # ∂h₂/∂y  -> y
#add_edge!(g, 9, 12,) # ∂h₂/∂y  -> y
#add_edge!(g, 1, 13,) # ∂h₂/∂x + ∂h₁/∂x  -> y
#add_edge!(g, 9, 14,) # "∂z/∂x"  -> ∂z/∂h₂
#add_edge!(g, 11, 14,) # "∂z/∂x"  -> ∂z/∂h₁
#add_edge!(g, 13, 14,) # "∂z/∂x"  -> ∂h₂/∂x + ∂h₁/∂x
# graphplot(adjacency_matrix(g),
# 	names = [vertices[i].name for i in 1:n],
# 	x = 0.25 .* [vertices[i].x for i in 1:n],
# 	y = 0.25 .* [vertices[i].y falseor i in 1:n],
# 	curves=false,
# 	markercolor = [vertices[i].color for i in 1:n],
# )

for n in [6, 7, 9, 11, 13, 14]
    #graphplot(adjacency_matrix(g)[1:n, 1:n],
	#	names = [vertices[i].name for i in 1:n],
	#	x = 0.25 .* [vertices[i].x for i in 1:n],
	#	y = 0.25 .* [vertices[i].y for i in 1:n],
	#	curves=false,
	#	markercolor = [vertices[i].color for i in 1:n],
	#)
	names = [vertices[i].name for i in 1:length(vertices)]
    t = TikzGraphs.plot(g, names)
    TikzPictures.save(SVG("graphdiff_$(n).svg"),t)
end


text:
struct DivConformity <: Conformity end

"""
    RaviartThomasRefFE(::Type{et},p::Polytope,order::Integer) where et

The `order` argument has the following meaning: the divergence of the  functions in this basis
is in the Q space of degree `order`.
"""
function RaviartThomasRefFE(::Type{et},p::Polytope,order::Integer) where et

  D = num_dims(p)

  prebasis = QCurlGradMonomialBasis{D}(et,order)

  nf_nodes, nf_moments = _RT_nodes_and_moments(et,p,order)

  face_own_dofs = _face_own_dofs_from_moments(nf_moments)

  face_dofs = face_own_dofs

  dof_basis = MomentBasedDofBasis(nf_nodes, nf_moments)

  ndofs = num_dofs(dof_basis)

  metadata = nothing

  reffe = GenericRefFE(
    ndofs,
    p,
    prebasis,
    dof_basis,
    DivConformity(),
    metadata,
    face_dofs)

  reffe
end

function get_face_own_dofs(reffe::GenericRefFE{DivConformity}, conf::DivConformity)
  get_face_dofs(reffe)
end

function _RT_nodes_and_moments(::Type{et}, p::Polytope, order::Integer) where et

  @notimplementedif ! is_n_cube(p)

  D = num_dims(p)
  ft = VectorValue{D,et}
  pt = Point{D,et}

  nf_nodes = [ zeros(pt,0) for face in 1:num_faces(p)]
  nf_moments = [ zeros(ft,0,0) for face in 1:num_faces(p)]

  fcips, fmoments = _RT_face_values(p,et,order)
  frange = get_dimrange(p,D-1)
  nf_nodes[frange] = fcips
  nf_moments[frange] = fmoments

  if (order > 0)
    ccips, cmoments = _RT_cell_values(p,et,order)
    crange = get_dimrange(p,D)
    nf_nodes[crange] = ccips
    nf_moments[crange] = cmoments
  end

  nf_nodes, nf_moments
end

# Ref FE to faces geomaps
function _ref_face_to_faces_geomap(p,fp)
  cfvs = get_face_coordinates(p,num_dims(fp))
  nc = length(cfvs)
  freffe = LagrangianRefFE(Float64,fp,1)
  fshfs = get_shapefuns(freffe)
  cfshfs = fill(fshfs, nc)
  fgeomap = lincomb(cfshfs,cfvs)
end

function _nfaces_evaluation_points_weights(p, fgeomap, fips, wips)
  nc = length(fgeomap)
  c_fips = fill(fips,nc)
  c_wips = fill(wips,nc)
  pquad = evaluate(fgeomap,c_fips)
  c_fips, pquad, c_wips
end

function _broadcast(::Type{T},n,b) where T
  c = Array{T}(undef,size(b))
  for (ii, i) in enumerate(b)
    c[ii] = i*n
  end
  return c
end

function _RT_face_moments(p, fshfs, c_fips, fcips, fwips)
  nc = length(c_fips)
  cfshfs = fill(fshfs, nc)
  cvals = evaluate(cfshfs,c_fips)
  cvals = [fwips[i].*cvals[i] for i in 1:nc]
  # fns, os = get_facet_normals(p)
  fns = get_facet_normals(p)
  os = get_facet_orientations(p)
  # @santiagobadia : Temporary hack for making it work for structured hex meshes
  cvals = [ _broadcast(typeof(n),n*o,b) for (n,o,b) in zip(fns,os,cvals)]
  #cvals = [ _broadcast(typeof(n),n,b) for (n,b) in zip(fns,cvals)]
  return cvals
end

# It provides for every face the nodes and the moments arrays
function _RT_face_values(p,et,order)

  # Reference facet
  @assert is_simplex(p) || is_n_cube(p) "We are assuming that all n-faces of the same n-dim are the same."
  fp = Polytope{num_dims(p)-1}(p,1)

  # geomap from ref face to polytope faces
  fgeomap = _ref_face_to_faces_geomap(p,fp)

  # Nodes are integration points (for exact integration)
  # Thus, we define the integration points in the reference
  # face polytope (fips and wips). Next, we consider the
  # n-face-wise arrays of nodes in fp (constant cell array c_fips)
  # the one of the points in the polytope after applying the geopmap
  # (fcips), and the weights for these nodes (fwips, a constant cell array)
  # Nodes (fcips)
  degree = (order+1)*2
  fquad = Quadrature(fp,degree)
  fips = get_coordinates(fquad)
  wips = get_weights(fquad)

  c_fips, fcips, fwips = _nfaces_evaluation_points_weights(p, fgeomap, fips, wips)

  # Moments (fmoments)
  # The RT prebasis is expressed in terms of shape function
  fshfs = MonomialBasis(et,fp,order)

  # Face moments, i.e., M(Fi)_{ab} = q_RF^a(xgp_RFi^b) w_Fi^b n_Fi ⋅ ()
  fmoments = _RT_face_moments(p, fshfs, c_fips, fcips, fwips)

  return fcips, fmoments

end

function _RT_cell_moments(p, cbasis, ccips, cwips)
  # Interior DOFs-related basis evaluated at interior integration points
  ishfs_iips = evaluate(cbasis,ccips)
  return cwips.*ishfs_iips
end

# It provides for every cell the nodes and the moments arrays
function _RT_cell_values(p,et,order)
  # Compute integration points at interior
  degree = 2*(order+1)
  iquad = Quadrature(p,degree)
  ccips = get_coordinates(iquad)
  cwips = get_weights(iquad)

  # Cell moments, i.e., M(C)_{ab} = q_C^a(xgp_C^b) w_C^b ⋅ ()
  cbasis = QGradMonomialBasis{num_dims(p)}(et,order-1)
  cmoments = _RT_cell_moments(p, cbasis, ccips, cwips )

  return [ccips], [cmoments]

end

function _face_own_dofs_from_moments(f_moments)
  face_dofs = Vector{Int}[]
  o = 1
  for moments in f_moments
    ndofs = size(moments,2)
    dofs = collect(o:(o+ndofs-1))
    push!(face_dofs,dofs)
    o += ndofs
  end
  face_dofs
end

struct MomentBasedDofBasis{P,V} <: Dof
  nodes::Vector{P}
  face_moments::Vector{Array{V}}
  face_nodes::Vector{UnitRange{Int}}

  function MomentBasedDofBasis(nodes,f_moments,f_nodes)
    P = eltype(nodes)
    V = eltype(eltype(f_moments))
    new{P,V}(nodes,f_moments,f_nodes)
  end

  function MomentBasedDofBasis(f_nodes,f_moments)
    P = eltype(eltype(f_nodes))
    V = eltype(eltype(f_moments))
    nodes = P[]
    face_nodes = UnitRange{Int}[]
    nfaces = length(f_nodes)
    n = 1
    for fi in 1:nfaces
      nodes_fi = f_nodes[fi]
      nini = n
      for node_fi in nodes_fi
        push!(nodes,node_fi)
        n += 1
      end
      nend = n-1
      push!(face_nodes,nini:nend)
    end
    new{P,V}(nodes,f_moments,face_nodes)
  end
end

get_nodes(b::MomentBasedDofBasis) = b.nodes
get_face_moments(b::MomentBasedDofBasis) = b.face_moments
get_face_nodes_dofs(b::MomentBasedDofBasis) = b.face_nodes

function num_dofs(b::MomentBasedDofBasis)
  n = 0
  for m in b.face_moments
    n += size(m,2)
  end
  n
end

function dof_cache(b::MomentBasedDofBasis,field)
  cf = field_cache(field,b.nodes)
  vals = evaluate_field!(cf,field,b.nodes)
  ndofs = num_dofs(b)
  r = _moment_dof_basis_cache(vals,ndofs)
  c = CachedArray(r)
  (c, cf)
end

function _moment_dof_basis_cache(vals::AbstractVector,ndofs)
  T = eltype(vals)
  r = zeros(eltype(T),ndofs)
end

function _moment_dof_basis_cache(vals::AbstractMatrix,ndofs)
  _, npdofs = size(vals)
  T = eltype(vals)
  r = zeros(eltype(T),ndofs,npdofs)
end

function evaluate_dof!(cache,b::MomentBasedDofBasis,field)
  c, cf = cache
  vals = evaluate_field!(cf,field,b.nodes)
  dofs = c.array
  _eval_moment_dof_basis!(dofs,vals,b)
  dofs
end

function _eval_moment_dof_basis!(dofs,vals::AbstractVector,b)
  o = 1
  z = zero(eltype(dofs))
  face_nodes = b.face_nodes
  face_moments = b.face_moments
  for face in 1:length(face_moments)
    moments = face_moments[face]
    if length(moments) != 0
      nodes = face_nodes[face]
      ni,nj = size(moments)
      for j in 1:nj
        dofs[o] = z
        for i in 1:ni
          dofs[o] += moments[i,j]⋅vals[nodes[i]]
        end
        o += 1
      end
    end
  end
end

function _eval_moment_dof_basis!(dofs,vals::AbstractMatrix,b)
  o = 1
  na = size(vals,2)
  z = zero(eltype(dofs))
  face_nodes = b.face_nodes
  face_moments = b.face_moments
  for face in 1:length(face_moments)
    moments = face_moments[face]
    if length(moments) != 0
      nodes = face_nodes[face]
      ni,nj = size(moments)
      for j in 1:nj
        for a in 1:na
          dofs[o,a] = z
          for i in 1:ni
            dofs[o,a] += moments[i,j]⋅vals[nodes[i],a]
          end
        end
        o += 1
      end
    end
  end
end


text:
mutable struct IndLSTM <: Block
    # input control gate params
    wi::VarOrNil
    ui::VarOrNil
    bi::VarOrNil
    # forget control gate params
    wf::VarOrNil
    uf::VarOrNil
    bf::VarOrNil
    # out control gate params
    wo::VarOrNil
    uo::VarOrNil
    bo::VarOrNil
    # new cell info params
    wc::VarOrNil
    uc::VarOrNil
    bc::VarOrNil
    h  # hidden variable
    c  #   cell variable
    function IndLSTM(inputSize::Int, hiddenSize::Int; type::Type=Array{Float32})
        T  = eltype(type)
        A  = T(1E-1)

        wi = randn(T, hiddenSize, inputSize) .* sqrt( T(2/inputSize) )
        ui = randn(T, hiddenSize, 1) .* A
        bi = 3ones(T, hiddenSize, 1)

        wf = randn(T, hiddenSize, inputSize) .* sqrt( T(2/inputSize) )
        uf = randn(T, hiddenSize, 1) .* A
        bf =-3ones(T, hiddenSize, 1)

        wo = randn(T, hiddenSize, inputSize) .* sqrt( T(2/inputSize) )
        uo = randn(T, hiddenSize, 1) .* A
        bo = 3ones(T, hiddenSize, 1)

        wc = randn(T, hiddenSize, inputSize) .* sqrt( T(2/inputSize) )
        uc = randn(T, hiddenSize, 1) .* A
        bc = zeros(T, hiddenSize, 1)

        new(Variable{type}(wi,true,true,true), Variable{type}(ui,true,true,true), Variable{type}(bi,true,true,true),
            Variable{type}(wf,true,true,true), Variable{type}(uf,true,true,true), Variable{type}(bf,true,true,true),
            Variable{type}(wo,true,true,true), Variable{type}(uo,true,true,true), Variable{type}(bo,true,true,true),
            Variable{type}(wc,true,true,true), Variable{type}(uc,true,true,true), Variable{type}(bc,true,true,true), nothing, nothing)
    end
    function IndLSTM()
        new(nothing, nothing, nothing,
            nothing, nothing, nothing,
            nothing, nothing, nothing,
            nothing, nothing, nothing, nothing, nothing)
    end
end

function clone(this::IndLSTM; type::Type=Array{Float32})
    cloned = IndLSTM()

    cloned.wi = clone(this.wi, type=type)
    cloned.bi = clone(this.bi, type=type)
    cloned.ui = clone(this.ui, type=type)

    cloned.wf = clone(this.wf, type=type)
    cloned.bf = clone(this.bf, type=type)
    cloned.uf = clone(this.uf, type=type)

    cloned.wo = clone(this.wo, type=type)
    cloned.bo = clone(this.bo, type=type)
    cloned.uo = clone(this.uo, type=type)

    cloned.wc = clone(this.wc, type=type)
    cloned.bc = clone(this.bc, type=type)
    cloned.uc = clone(this.uc, type=type)

    return cloned
end

mutable struct IndLSTMs <: Block
    layers::Vector{IndLSTM}
    function IndLSTMs(topology::Vector{Int}; type::Type=Array{Float32})
        n = length(topology) - 1
        layers = Vector{IndLSTM}(undef, n)
        for i = 1:n
            layers[i] = IndLSTM(topology[i], topology[i+1]; type=type)
        end
        new(layers)
    end
end

Base.getindex(m::IndLSTMs,     k...) =  m.layers[k...]
Base.setindex!(m::IndLSTMs, v, k...) = (m.layers[k...] = v)
Base.length(m::IndLSTMs)       = length(m.layers)
Base.lastindex(m::IndLSTMs)    = length(m.layers)
Base.firstindex(m::IndLSTMs)   = 1
Base.iterate(m::IndLSTMs, i=firstindex(m)) = i>length(m) ? nothing : (m[i], i+1)

function Base.show(io::IO, m::IndLSTM)
    SIZE = size(m.wi)
    TYPE = typeof(m.wi.value)
    print(io, "IndLSTM($(SIZE[2]), $(SIZE[1]); type=$TYPE)")
end

function Base.show(io::IO, model::IndLSTMs)
    for m in model
        show(io, m)
    end
end

function resethidden(model::IndLSTM)
    model.h = nothing
    model.c = nothing
end

function resethidden(model::IndLSTMs)
    for m in model
        resethidden(m)
    end
end

function forward(model::IndLSTM, x::Variable{T}) where T
    wi = model.wi
    ui = model.ui
    bi = model.bi

    wf = model.wf
    uf = model.uf
    bf = model.bf

    wo = model.wo
    uo = model.uo
    bo = model.bo

    wc = model.wc
    uc = model.uc
    bc = model.bc

    h = model.h ≠ nothing ? model.h : Variable{T}(Zeros(T, size(wi,1), size(x,2)))
    c = model.c ≠ nothing ? model.c : Variable{T}(Zeros(T, size(wc,1), size(x,2)))

    z = tanh(    matAddVec(wc * x + matMulVec(h,uc), bc) )
    i = sigmoid( matAddVec(wi * x + matMulVec(h,ui), bi) )
    f = sigmoid( matAddVec(wf * x + matMulVec(h,uf), bf) )
    o = sigmoid( matAddVec(wo * x + matMulVec(h,uo), bo) )
    c = dotMul(f, c) + dotMul(i, z)
    h = dotMul(o, tanh(c))

    model.c = c
    model.h = h

    return h
end

function forward(model::IndLSTMs, x::Variable)
    for m in model
        x = forward(m, x)
    end
    return x
end

function predict(model::IndLSTM, x::T) where T
    wi = model.wi.value
    ui = model.ui.value
    bi = model.bi.value

    wf = model.wf.value
    uf = model.uf.value
    bf = model.bf.value

    wo = model.wo.value
    uo = model.uo.value
    bo = model.bo.value

    wc = model.wc.value
    uc = model.uc.value
    bc = model.bc.value

    h = model.h ≠ nothing ? model.h : Zeros(T, size(wi,1), size(x,2))
    c = model.c ≠ nothing ? model.c : Zeros(T, size(wc,1), size(x,2))

    z = tanh(    wc * x + h .* uc .+ bc )
    i = sigmoid( wi * x + h .* ui .+ bi )
    f = sigmoid( wf * x + h .* uf .+ bf )
    o = sigmoid( wo * x + h .* uo .+ bo )
    c = f .* c + i .* z
    h = o .* tanh(c)

    model.c = c
    model.h = h

    return h
end

function predict(model::IndLSTMs, x)
    for m in model
        x = predict(m, x)
    end
    return x
end

"""
    unbiasedof(m::IndLSTM)

unbiased weights of IndLSTM block
"""
function unbiasedof(m::IndLSTM)
    weights = Vector(undef, 4)
    weights[1] = m.wi.value
    weights[2] = m.wf.value
    weights[3] = m.wo.value
    weights[4] = m.wc.value
    return weights
end

function weightsof(m::IndLSTM)
    weights = Vector{Variable}(undef,12)
    weights[1] = m.wi.value
    weights[2] = m.ui.value
    weights[3] = m.bi.value

    weights[4] = m.wf.value
    weights[5] = m.uf.value
    weights[6] = m.bf.value

    weights[7] = m.wo.value
    weights[8] = m.uo.value
    weights[9] = m.bo.value

    weights[10] = m.wc.value
    weights[11] = m.uc.value
    weights[12] = m.bc.value
    return weights
end

"""
    unbiasedof(model::IndLSTMs)

unbiased weights of IndLSTMs block
"""
function unbiasedof(model::IndLSTMs)
    weights = Vector(undef, 0)
    for m in model
        append!(weights, unbiasedof(m))
    end
    return weights
end

function weightsof(model::IndLSTMs)
    weights = Vector(undef,0)
    for m in model
        append!(weights, weightsof(m))
    end
    return weights
end

function gradsof(m::IndLSTM)
    grads = Vector{Variable}(undef,12)
    grads[1] = m.wi.delta
    grads[2] = m.ui.delta
    grads[3] = m.bi.delta

    grads[4] = m.wf.delta
    grads[5] = m.uf.delta
    grads[6] = m.bf.delta

    grads[7] = m.wo.delta
    grads[8] = m.uo.delta
    grads[9] = m.bo.delta

    grads[10] = m.wc.delta
    grads[11] = m.uc.delta
    grads[12] = m.bc.delta
    return grads
end

function gradsof(model::IndLSTMs)
    grads = Vector(undef,0)
    for m in model
        append!(grads, gradsof(m))
    end
    return grads
end

function zerograds!(m::IndLSTM)
    for v in gradsof(m)
        v .= zero(v)
    end
end

function zerograds!(m::IndLSTMs)
    for v in gradsof(m)
        v .= zero(v)
    end
end

function paramsof(m::IndLSTM)
    params = Vector{Variable}(undef,12)
    params[1] = m.wi
    params[2] = m.ui
    params[3] = m.bi

    params[4] = m.wf
    params[5] = m.uf
    params[6] = m.bf

    params[7] = m.wo
    params[8] = m.uo
    params[9] = m.bo

    params[10] = m.wc
    params[11] = m.uc
    params[12] = m.bc
    return params
end

function xparamsof(m::IndLSTM)
    xparams = Vector{XVariable}(undef,12)
    xparams[1] = ('w', m.wi)
    xparams[2] = ('u', m.ui)
    xparams[3] = ('b', m.bi)

    xparams[4] = ('w', m.wf)
    xparams[5] = ('u', m.uf)
    xparams[6] = ('b', m.bf)

    xparams[7] = ('w', m.wo)
    xparams[8] = ('u', m.uo)
    xparams[9] = ('b', m.bo)

    xparams[10] = ('w', m.wc)
    xparams[11] = ('u', m.uc)
    xparams[12] = ('b', m.bc)
    return xparams
end

function paramsof(model::IndLSTMs)
    params = Vector{Variable}(undef,0)
    for m in model
        append!(params, paramsof(m))
    end
    return params
end

function xparamsof(model::IndLSTMs)
    xparams = Vector{XVariable}(undef,0)
    for m in model
        append!(xparams, xparamsof(m))
    end
    return xparams
end

function nparamsof(m::IndLSTM)
    lw = length(m.wi)
    lu = length(m.ui)
    lb = length(m.bi)
    return (lw+lu+lb)*4
end

function bytesof(model::IndLSTM, unit::String="MB")
    n = nparamsof(model) * elsizeof(model.wi)
    return blocksize(n, uppercase(unit))
end

function nparamsof(model::IndLSTMs)
    num = 0
    for m in model
        num += nparamsof(m)
    end
    return num
end

function bytesof(model::IndLSTMs, unit::String="MB")
    n = nparamsof(model) * elsizeof(model[1].wi)
    return blocksize(n, uppercase(unit))
end

function to(type::Type, m::IndLSTM)
    m.wi = to(type, m.wi)
    m.ui = to(type, m.ui)
    m.bi = to(type, m.bi)

    m.wf = to(type, m.wf)
    m.uf = to(type, m.uf)
    m.bf = to(type, m.bf)

    m.wo = to(type, m.wo)
    m.uo = to(type, m.uo)
    m.bo = to(type, m.bo)

    m.wc = to(type, m.wc)
    m.uc = to(type, m.uc)
    m.bc = to(type, m.bc)
    return m
end

function to!(type::Type, m::IndLSTM)
    m = to(type, m)
    return nothing
end

function to(type::Type, m::IndLSTMs)
    for layer in m
        layer = to(type, layer)
    end
    return m
end

function to!(type::Type, m::IndLSTMs)
    for layer in m
        to!(type, layer)
    end
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
import PackageCompiler, Pkg
Pkg.activate(".")
Pkg.instantiate()
Pkg.precompile()

PackageCompiler.create_sysimage(
    ["Plots", "CSV", "StatsPlots", "Markdown", "DataFrames", "Query", "ColorSchemes", "StatsBase", "Random", "Distributions"];
    sysimage_path = "DSSysimage.so",
    precompile_execution_file = "src/EnvBaker.jl"
)

max_stars_repo_name: Jblew/towards-julia-notebook

max_issues_repo_name: Jblew/towards-julia-notebook

max_forks_repo_name: Jblew/towards-julia-notebook


text:
using MathOptInterfaceMosek

using Base.Test

const MOI = MathOptInterface

include(joinpath(Pkg.dir("MathOptInterface"),"test","contlinear.jl"))
@testset "Continuous linear problems" begin
    contlineartest(MosekSolver(QUIET = true))
end

# include("contquadratic.jl")
# @testset "Continuous quadratic problems" begin
#     # contquadratictest(GurobiSolver())
# end

include(joinpath(Pkg.dir("MathOptInterface"),"test","contconic.jl"))
@testset "Continuous conic problems" begin
    contconictest(MosekSolver(QUIET = true))
end

include(joinpath(Pkg.dir("MathOptInterface"),"test","intlinear.jl"))
@testset "Mixed-integer linear problems" begin
    intlineartest(MosekSolver(QUIET = true))
end


text:
@testset "pointer" begin

# inner constructors

voidptr_a = Ptr{Cvoid}(Int(0xDEADBEEF))
generic_voidptr_a = AMDGPUnative.DevicePtr{Cvoid,AS.Generic}(voidptr_a)
global_voidptr_a = AMDGPUnative.DevicePtr{Cvoid,AS.Global}(voidptr_a)
local_voidptr_a = AMDGPUnative.DevicePtr{Cvoid,AS.Local}(voidptr_a)

voidptr_b = Ptr{Cvoid}(Int(0xCAFEBABE))
generic_voidptr_b = AMDGPUnative.DevicePtr{Cvoid,AS.Generic}(voidptr_b)
global_voidptr_b = AMDGPUnative.DevicePtr{Cvoid,AS.Global}(voidptr_b)
local_voidptr_b = AMDGPUnative.DevicePtr{Cvoid,AS.Local}(voidptr_b)

intptr_b = convert(Ptr{Int}, voidptr_b)
generic_intptr_b = AMDGPUnative.DevicePtr{Int,AS.Generic}(intptr_b)
global_intptr_b = AMDGPUnative.DevicePtr{Int,AS.Global}(intptr_b)
local_intptr_b = AMDGPUnative.DevicePtr{Int,AS.Local}(intptr_b)

# outer constructors
@test AMDGPUnative.DevicePtr{Cvoid}(voidptr_a) == generic_voidptr_a
@test AMDGPUnative.DevicePtr(voidptr_a) == generic_voidptr_a

# getters
@test eltype(generic_voidptr_a) == Cvoid
@test eltype(global_intptr_b) == Int
@test AMDGPUnative.addrspace(generic_voidptr_a) == AS.Generic
@test AMDGPUnative.addrspace(global_voidptr_a) == AS.Global
@test AMDGPUnative.addrspace(local_voidptr_a) == AS.Local

# comparisons
@test generic_voidptr_a != global_voidptr_a
@test generic_voidptr_a != generic_intptr_b

@testset "conversions" begin

# between regular and device pointers

@test convert(Ptr{Cvoid}, generic_voidptr_a) == voidptr_a
@test convert(AMDGPUnative.DevicePtr{Cvoid}, voidptr_a) == generic_voidptr_a
@test convert(AMDGPUnative.DevicePtr{Cvoid,AS.Global}, voidptr_a) == global_voidptr_a

# between device pointers

@test_throws ArgumentError convert(typeof(local_voidptr_a), global_voidptr_a)
@test convert(typeof(generic_voidptr_a), generic_voidptr_a) == generic_voidptr_a
@test convert(typeof(global_voidptr_a), global_voidptr_a) == global_voidptr_a
@test Base.unsafe_convert(typeof(local_voidptr_a), global_voidptr_a) == local_voidptr_a

@test convert(typeof(global_voidptr_a), global_intptr_b) == global_voidptr_b
@test convert(typeof(generic_voidptr_a), global_intptr_b) == generic_voidptr_b
@test convert(typeof(global_voidptr_a), generic_intptr_b) == global_voidptr_b

@test convert(AMDGPUnative.DevicePtr{Cvoid}, global_intptr_b) == global_voidptr_b

end

end


text:
energy(z::AbstractMatrix, θ) = energy(z, θ, EnergyBuffer(θ, 1))
energy(z::AbstractArray{<:Any,3}, θ) = energy(z, θ, EnergyBuffer(θ, size(z, 3)))
function energy(z, θ::SpinModel, buffer)
    (; E, A, A₁₂, B, B₁) = buffer
    (; q, N) = θ
    qN = q*N

    # E[m] = h[μ,i] * Z[μ,i,m]
    if hasproperty(θ, :h)
        mul!(vec(E), transpose(reshape(z, qN, :)), vec(θ.h))
    else
        fill!(E, zero(eltype(E)))
    end

    # E[m]+= Z[μ,i,m] J[μ,i,ν,j] Z[ν,j,m]
    if hasproperty(θ, :J)
        mul!(reshape(A, qN,:), reshape(θ.J, qN,qN), reshape(z, qN,:))
        A .*= z
        sum!(A₁₂, A)
        E .+= vec(A₁₂)
    end

    # E[m]+= logcosh(B[a,m])  &  B[a,m] = W[a,μ,i] Z[μ,i,m] + b[a]
    if hasproperty(θ, :W)
        mul!(B, reshape(θ.W, :, qN), reshape(z, qN, :))
        @. B = logcosh(θ.b + B)
        sum!(B₁, B)
        rmul!(B₁, θ.rbm_wt)
        E .+= vec(B₁)
    end

    return rmul!(E, -1)
end

∂x_energy(z::AbstractMatrix, θ) = ∂x_energy(z, θ, EnergyBuffer(θ, 1))
∂x_energy(z::AbstractArray{<:Any,3}, θ) = ∂x_energy(z, θ, EnergyBuffer(θ, size(z, 3)))
function ∂x_energy(z, θ::SpinModel, buffer)
    (; A, B, K) = buffer
    (; q, N) = θ
    qN = q*N
    ∇E = A         # non-allocating alias

    if hasproperty(θ, :h)
        @. ∇E = θ.h    # ∇E[μim] = h[μi]
    else
        fill!(∇E, zero(eltype(∇E)))
    end

    # ∇E[μim]+= (J[μiνj] + J[νjμi]) * Z[νjm]
    if hasproperty(θ, :J)
        K .= reshape(θ.J, qN,qN) .+ transpose(reshape(θ.J, qN,qN))
        mul!(reshape(∇E, qN, :), K, reshape(z, qN, :), 1, 1)
    end

    # ∇E[μim]+= W[aμi] * tanh(W[aνj]*Z[νjm] + b[a])
    if hasproperty(θ, :W)
        mul!(B, reshape(θ.W, :, qN), reshape(z, qN, :))
        @. B = tanh_fast(B + θ.b)
        mul!(reshape(∇E, qN, :), transpose(reshape(θ.W, :, qN)), B, θ.rbm_wt, 1)
    end

    return rmul!(∇E, -1)
end

struct EnergyBuffer{V,T3,M}
    E::V
    A::T3
    A₁₂::T3
    B::M
    B₁::M
    K::M
    function EnergyBuffer(θ::SpinModel{T}, M::Int) where T
        (; q, N) = θ
        has_h, has_J, has_W = (hasproperty(θ, s) for s ∈ (:h, :J, :W))
        X = has_h ? θ.h : has_J ? θ.J : θ.W
        E = similar(X, M)
        A = similar(X, q, N, M)
        A₁₂ = similar(X, 1, 1, M)
        B = has_W ? similar(X, θ.P, M) : similar(X, 0, 0)
        B₁ = has_W ? similar(X, 1, M) : similar(X, 0, 0)
        K = has_J ? similar(X, q*N, q*N) : similar(X, 0, 0)
        return new{typeof(E),typeof(A),typeof(B)}(E, A, A₁₂, B, B₁, K)
    end
    EnergyBuffer(z::AbstractArray{<:Any,3}, θ::SpinModel) = EnergyBuffer(θ, size(z,3))
end

struct ΔEnergyBuffer{T3,M,T4}
    ΔE::T3      # q × N × M
    A::T3       # q × N × M
    A₁::T3      # 1 × N × M     (for J & h)
    K::M        # qN × qN       (for J)
    B::M        # P × M         (for W)
    B₁::M       # 1 × M         (for W)
    B̃::T4       # P × q × N × M (for W)
    B̃₁::T4      # 1 × q × N × M (for W)
    B̃₂::T4      # P × 1 × N × M (for W)
    function ΔEnergyBuffer(θ::SpinModel, M::Int)
        (;q, N) = θ
        has_h, has_J, has_W = (hasproperty(θ, s) for s ∈ (:h, :J, :W))
        X = has_h ? θ.h : has_J ? θ.J : θ.W
        ΔE, A = ntuple(i -> similar(X, q,N,M), 2)
        A₁ = (has_h | has_J) ? similar(X, 1,N,M) : similar(X, 0,0,0)
        K  = has_J ? similar(X, q*N,q*N) : similar(X, 0,0)
        B  = has_W ? similar(X, θ.P,M) : similar(X, 0,0)
        B₁ = has_W ? similar(X, 1,M) : similar(X, 0,0)
        B̃  = has_W ? similar(X, θ.P,q,N,M) : similar(X, 0,0,0,0)
        B̃₁ = has_W ? similar(X, 1,q,N,M) : similar(X, 0,0,0,0)
        B̃₂ = has_W ? similar(X, θ.P,1,N,M) : similar(X, 0,0,0,0)
        return new{typeof(A), typeof(K), typeof(B̃)}(
            ΔE, A, A₁, K, B, B₁, B̃, B̃₁, B̃₂
        )
    end
    ΔEnergyBuffer(z::AbstractArray{<:Any,3}, θ::SpinModel) = ΔEnergyBuffer(θ, size(z,3))
end

# Compute energy changes from one spin flip:    ΔE[ϕ,k,m] = E(zᵠᵏ) - E(z)
ΔE₁(z, θ::SpinModel) = ΔE₁(z, θ, ΔEnergyBuffer(z, θ))
function ΔE₁(z, θ::SpinModel, buffer)
    (; ΔE, A, A₁, K, B, B₁, B̃, B̃₁, B̃₂) = buffer
    fill!(ΔE, zero(eltype(ΔE))) # initialize ΔE
    q, N, M = size(z)

    if hasproperty(θ, :J)
        J = θ.J
        K .= reshape(J, q*N, q*N)
        K.+= transpose(reshape(J, q*N, q*N))
        mul!(reshape(ΔE, :,M), K, reshape(z, :,M), -1, 1)
        @. A = ΔE * z
        ΔE .-= sum!(A₁, A)
    end

    if hasproperty(θ, :h)
        h = θ.h
        @. A = h * z
        sum!(A₁, A)
        @. ΔE += -(h - A₁)
    end

    if hasproperty(θ, :W) && hasproperty(θ, :b)
        (;W, b) = θ
        mul!(B, reshape(W, :, q*N), reshape(z, :, M))   # B[a,m] = ∑ᵨₖ W[a,ρ,k] * z[ρ,k,m]
        B.+= b                                          # B[a,m] = b[a] + ∑ᵨₖ W[a,ρ,k] * z[ρ,k,m]
        B̃ .= W .* reshape(z, 1,q,N,M)   # B̃[a,ρ,k,m] = W[a,ρ,k] * z[ρ,k,m]
        sum!(B̃₂, B̃)                     # B̃₂[a,.,k,m] = ∑ᵨ W[a,ρ,k] * z[ρ,k,m]
        @. B̃ = W - B̃₂                   # B̃[a,σ,k,m] = W[a,σ,k] - ∑ᵨ W[a,ρ,k] * z[ρ,k,m]
        B̃.+= reshape(B, :,1,1,M)        # B̃[a,σ,k,m] = B[a,m] + W[a,σ,k] - ∑ᵨ W[a,ρ,k] * z[ρ,k,m]
        ΔE.+= θ.rbm_wt .* reshape(sum!(logcosh, B₁, B), 1,1,M)
        ΔE.-= θ.rbm_wt .* reshape(sum!(logcosh, B̃₁, B̃), q,N,M)
    end

    return ΔE
end


text:
# Evaluation
function f(xp::AbstractVector, α::T, x::AbstractVector, v::AbstractVector, β::T) where {T}
	Δ = (xp - (x + α*v))
	return 0.5*β*Δ'*Δ
end

function c(x::AbstractVector, ascore::ActiveSetCore, prob::Algames.GameProblem)
	probsize = prob.probsize
	N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	S = probsize.S
	Sh = S + p*(p-1)*(N-1)
	Sv = S + Int(p*(p-1)*(N-1)/2)

	set_primal_dual!(ascore, prob, x)
	residual!(ascore, prob, prob.pdtraj)
	out = deepcopy(ascore.res)
	return out
end

function L(xp::AbstractVector, α::T, λ::AbstractVector, x::AbstractVector,
	v::AbstractVector, β::T, ascore::ActiveSetCore, prob::Algames.GameProblem, vmask::Vector{Int}) where {T}

	out = f(xp, α, x, v, β) + λ[vmask]'*c(xp, ascore, prob)[vmask]
	return out
end

# Gradients
function Algames.residual_jacobian!(ascore::ActiveSetCore, prob::Algames.GameProblem{KN,n,m,T,SVd,SVx},
	pdtraj::PrimalDualTraj{KN,n,m,T,SVd}) where {KN,n,m,T,SVd,SVx}

	probsize = ascore.probsize
	N = probsize.N
	p = probsize.p
	S = probsize.S

	# Reset!
	sparse_zero!(ascore.jac)

	residual_jacobian!(prob, pdtraj)
	# regularize_residual_jacobian!(prob) #############################
	ascore.jac[1:S,1:S] .= prob.core.jac

    # Allocations
	vs = VStamp()
	hs = HStamp()
    cs = CStamp()

	jacobian!(prob.game_con, pdtraj.pr)
	for i = 1:p
		for j ∈ setdiff(1:p,i)
			conval, v = get_collision_conval(prob.game_con, i, j)
			if v
				for (l,k) in enumerate(conval.inds)
					stampify!(vs, :opt, i, :x, 1, k)
					stampify!(cs, :h, :col, i, j, k)
					if valid(cs,N,p) && valid(vs,N,p)
						add2sub(ascore.jac_sub[(vs,cs)], conval.jac[l]')
					end
					if valid(hs,N,p) && valid(cs,N,p)
						stampify!(hs, :x, 1, k)
						stampify!(cs, :v, :col, i, j, k)
						add2sub(ascore.jac_sub[(hs,cs)], conval.jac[l])
					end
				end
			end
		end
	end
    return nothing
end

function ∇c(x::AbstractVector,ascore::ActiveSetCore, prob::Algames.GameProblem) where {T}
	N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	S = probsize.S
	Sh = S + p*(p-1)*(N-1)
	Sv = S + Int(p*(p-1)*(N-1)/2)

	∇ = zeros(Sv,Sh)
	set_primal_dual!(ascore, prob, x)
	residual_jacobian!(ascore, prob, prob.pdtraj)
	@assert all(size(ascore.jac) .== size(∇))
	∇ = ascore.jac
	return ∇
end

function grad(xp::AbstractVector, α::T, λ::AbstractVector, x::AbstractVector,
	v::AbstractVector, β::T, ascore::ActiveSetCore, prob::Algames.GameProblem, vmask::Vector{Int}) where {T}
	probsize = prob.probsize
	N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	S = probsize.S
	Sh = S + p*(p-1)*(N-1)
	Sv = S + Int(p*(p-1)*(N-1)/2)
	nx = Sh
	# nα = 1
	nλ = Sv
	ix = Vector(1:nx)
	# iα = Vector(nx .+ (1:nα))
	iλ = Vector(nx .+ (1:nλ))
	# iλ = Vector(nx + nα .+ (1:nλ))
	jmask = vcat(ix[vmask], iλ[vmask])

	∇ = ∇c(xp, ascore, prob)
	gx = β * (xp - (x+α*v)) + ∇'*λ
	# gx = (-β/length(x)^2 * (xp - (x+α*v)) + ∇'*λ)[vmask]
	# gα = [-β/length(x)^2 * (xp - (x+α*v))'*v - 1.0]
	gλ = c(xp, ascore, prob)
	g = vcat(gx, gλ)[jmask]
	# g = vcat(gx, gα, gλ)
	return g
end

# Jacobians

function jac(xp::AbstractVector, α::T, λ::AbstractVector, x::AbstractVector,
	v::AbstractVector, β::T, ascore::ActiveSetCore, prob::Algames.GameProblem, vmask::Vector{Int}) where {T}
	probsize = prob.probsize
	N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	S = probsize.S
	Sh = S + p*(p-1)*(N-1)
	Sv = S + Int(p*(p-1)*(N-1)/2)
	nx = Sh
	# nα = 1
	nλ = Sv
	ix = Vector(1:nx)
	# iα = Vector(nx .+ (1:nα))
	iλ = Vector(nx .+ (1:nλ))
	# iλ = Vector(nx + nα .+ (1:nλ))
	jmask = vcat(ix[vmask], iλ[vmask])
	# jmask = vcat(ix[vmask], iα, iλ[vmask])

	∇ = ∇c(xp, ascore, prob)
	j = zeros(nx+nλ, nx+nλ)
	# j = zeros(nx+nα+nλ, nx+nα+nλ)
	j[ix,ix] += β*I
	# j[ix,ix] += β/length(x)^2*I
	# j[iα,ix] = -β/length(x)^2*v'
	# j[ix,iα] = -β/length(x)^2*v
	j[iλ,ix] = ∇
	j[ix,iλ] = ∇'
	j[iλ,iλ] -= 1e-3*I

	return j[jmask, jmask]
end

# Newton Solver
function newton_solver(α::T, λ::AbstractVector, x::AbstractVector, v::AbstractVector,
	β::T, ascore::ActiveSetCore, prob::Algames.GameProblem, vmask::Vector{Int}) where {T}
	probsize = prob.probsize
	N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	S = probsize.S
	Sh = S + p*(p-1)*(N-1)
	Sv = S + Int(p*(p-1)*(N-1)/2)
	nx = Sh
	# nα = 1
	nλ = Sv
	ix = Vector(1:nx)
	# iα = Vector(nx .+ (1:nα))
	iλ = Vector(nx .+ (1:nλ))
	# iλ = Vector(nx + nα .+ (1:nλ))
	jmask = vcat(ix[vmask], iλ[vmask])
	# jmask = vcat(ix[vmask], iα, iλ[vmask])

	xp = deepcopy(x + α*v)

	for l = 1:1
		# plot!(prob.model, prob.pdtraj.pr)
		@show scn(mean(abs.(c(xp, ascore, prob_copy2)[vmask])))

		g = grad(xp, α, λ, x, v, β, ascore, prob, vmask)
		println("***********L = ", scn(L(xp, α, λ, x, v, β, ascore, prob, vmask)))
		println("***********res = ",scn(mean(abs.(g))))
		j = jac(xp, α, λ, x, v, β, ascore, prob, vmask) + 1e-18*I
		Δm = - \(lu(j), g)
		println("Δm = ", scn(mean(abs.(Δm))))

		Δ = zeros(nx+nλ)
		# Δ = zeros(nx+nα+nλ)
		# Δ[jmask] = 0.01*Δm
		Δ[jmask] = 0.1*α*Δm/mean(abs.(Δm[ix[vmask]]))
		# Δ[jmask] = 0.05*α*Δm/mean(abs.(Δm[ix[vmask]]))
		println("Δ = ", scn(mean(abs.(Δ[jmask]))))
		println("Δx = ", scn(mean(abs.(Δ[ix[vmask]]))))
		println("α = ", scn(α))
		xp[vmask] += Δ[ix[vmask]]
		# α += Δ[iα][1]
		λ[vmask] += Δ[iλ[vmask]]
		set_primal_dual!(ascore, prob, xp)
		#################
		prob.opts.shift = 0
		prob.opts.α_dual = 0.0
		prob.opts.αx_dual = 0.0*1e-3*ones(p)
		prob.opts.dual_reset = false
		prob.opts.ρ_0 = prob.pen.ρ[1]
		prob.opts.ρ_max = prob.pen.ρ[1]
		prob.opts.inner_print = false
		prob.opts.outer_print = false
		prob.opts.inner_iter = 20
		prob.opts.outer_iter = 60
		prob.opts.reg_0 = 1e-6
		ϵ = 1e-6
		prob.opts.ϵ_opt = ϵ
		prob.opts.ϵ_dyn = ϵ
		prob.opts.ϵ_con = ϵ
		prob.opts.ϵ_sta = ϵ
		newton_solve!(prob)
		xp = get_primal_dual(ascore, prob)
		println("f_no_change", scn(f(x, α, x, v, β)))
		println("f_solution ", scn(f(xp, α, x, v, β)))
		#################
		record!(prob.stats, prob.core, prob.model,
			prob.game_con, prob.pdtraj, prob.opts.outer_iter)
		opt = mean(abs.(prob.stats.opt_vio[end].vio))
		sta = mean(abs.(prob.stats.sta_vio[end].vio))
		dyn = mean(abs.(prob.stats.dyn_vio[end].vio))
		con = mean(abs.(prob.stats.con_vio[end].vio))
		println("opt = "*scn(opt))
		println("sta = "*scn(sta))
		println("dyn = "*scn(dyn))
		println("con = "*scn(con))
	end
	return xp, λ
end

function finite_difference_grad(xp, α, λ, x, v, β, ascore, prob, vmask)
	probsize = prob.probsize
	N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	S = probsize.S
	Sh = S + p*(p-1)*(N-1)
	Sv = S + Int(p*(p-1)*(N-1)/2)
	nx = Sh
	nλ = Sv
	ix = Vector(1:nx)
	iλ = Vector(nx .+ (1:nλ))
	jmask = vcat(ix[vmask], iλ[vmask])

	function Lλ(λ)
		return L(xp, α, λ, x, v, β, ascore, prob, vmask)
	end
	function Lx(xp)
		return L(xp, α, λ, x, v, β, ascore, prob, vmask)
	end

	gx = finite_difference_gradient(Lx, xp)
	gλ = finite_difference_gradient(Lλ, λ)
	g = vcat(gx, gλ)
	return g[jmask]
end

function finite_difference_gradient(f, x; eps::T=1e-5)
	# the gradient dimension will be m = length(x)
	m = length(x)
	x = SVector{m,T}(x)
	grad = zeros(m)
	for k = 1:m
		epsk = zeros(m)
		epsk[k] = eps
		grad[k] = (f(x .+ epsk) - f(x .- epsk))/(2eps)
	end
	return grad
end

function finite_difference_hessian(f, x; eps::T=1e-5)
	# the Hessian dimension will be m.m = length(x).length(x)
	m = length(x)
	x = SVector{m,T}(x)
	hess = zeros(m,m)
	for k = 1:m
		for l = 1:m
			epsk = zeros(m)
			epsl = zeros(m)
			epsk[k] = eps
			epsl[l] = eps
			hess[k,l] = (f(x .+ epsk .+ epsl) - f(x .+ epsk .- epsl) - f(x .- epsk .+ epsl) + f(x .- epsk .- epsl))/(4eps^2)
		end
	end
	return hess
end

function finite_difference_jacobian(f, x; n::Int=length(f(x)), eps::T=1e-5)
	# n is the dimension of the output of f
	# the Jacobian dimension will be n.m = n.length(x)
	m = length(x)
	x = SVector{m,T}(x)
	jac = zeros(n,m)
	for k = 1:m
		epsk = zeros(m)
		epsk[k] = eps
		jac[:,k] = (f(x .+ epsk) - f(x .- epsk))/(2eps)
	end
	return jac
end


text:
using FusionSystems
using Test

@testset "FusionSystems.jl" begin
    # Write your tests here.
end


text:
Base.convert(::Type{Cint}, x::GLFW.Key) = Cint(x)

@enum GlfwClientApi GlfwClientApi_Unknown GlfwClientApi_OpenGL GlfwClientApi_Vulkan

function ImGui_ImplGlfw_Init(window::GLFW.Window, install_callbacks::Bool, client_api::GlfwClientApi)
    g_Window[] = window
    g_Time[] = 0.0

    # setup back-end capabilities flags
    io = GetIO()
    io.BackendFlags |= ImGuiBackendFlags_HasMouseCursors
    io.BackendFlags |= ImGuiBackendFlags_HasSetMousePos
    io.BackendPlatformName = "imgui_impl_glfw"

    # keyboard mapping
    Set_KeyMap(io, ImGuiKey_Tab, GLFW.KEY_TAB)
    Set_KeyMap(io, ImGuiKey_LeftArrow, GLFW.KEY_LEFT)
    Set_KeyMap(io, ImGuiKey_RightArrow, GLFW.KEY_RIGHT)
    Set_KeyMap(io, ImGuiKey_UpArrow, GLFW.KEY_UP)
    Set_KeyMap(io, ImGuiKey_DownArrow, GLFW.KEY_DOWN)
    Set_KeyMap(io, ImGuiKey_PageUp, GLFW.KEY_PAGE_UP)
    Set_KeyMap(io, ImGuiKey_PageDown, GLFW.KEY_PAGE_DOWN)
    Set_KeyMap(io, ImGuiKey_Home, GLFW.KEY_HOME)
    Set_KeyMap(io, ImGuiKey_End, GLFW.KEY_END)
    Set_KeyMap(io, ImGuiKey_Insert, GLFW.KEY_INSERT)
    Set_KeyMap(io, ImGuiKey_Delete, GLFW.KEY_DELETE)
    Set_KeyMap(io, ImGuiKey_Backspace, GLFW.KEY_BACKSPACE)
    Set_KeyMap(io, ImGuiKey_Space, GLFW.KEY_SPACE)
    Set_KeyMap(io, ImGuiKey_Enter, GLFW.KEY_ENTER)
    Set_KeyMap(io, ImGuiKey_Escape, GLFW.KEY_ESCAPE)
    Set_KeyMap(io, ImGuiKey_A, GLFW.KEY_A)
    Set_KeyMap(io, ImGuiKey_C, GLFW.KEY_C)
    Set_KeyMap(io, ImGuiKey_V, GLFW.KEY_V)
    Set_KeyMap(io, ImGuiKey_X, GLFW.KEY_X)
    Set_KeyMap(io, ImGuiKey_Y, GLFW.KEY_Y)
    Set_KeyMap(io, ImGuiKey_Z, GLFW.KEY_Z)

    io.SetClipboardTextFn = g_ImplGlfw_SetClipboardText[]
    io.GetClipboardTextFn = g_ImplGlfw_GetClipboardText[]
    io.ClipboardUserData = Ptr{Cvoid}(g_Window[].handle)

    g_MouseCursors[ImGuiMouseCursor_Arrow+1] = GLFW.CreateStandardCursor(GLFW.ARROW_CURSOR)
    g_MouseCursors[ImGuiMouseCursor_TextInput+1] = GLFW.CreateStandardCursor(GLFW.IBEAM_CURSOR)
    g_MouseCursors[ImGuiMouseCursor_ResizeAll+1] = GLFW.CreateStandardCursor(GLFW.ARROW_CURSOR)
    g_MouseCursors[ImGuiMouseCursor_ResizeNS+1] = GLFW.CreateStandardCursor(GLFW.VRESIZE_CURSOR)
    g_MouseCursors[ImGuiMouseCursor_ResizeEW+1] = GLFW.CreateStandardCursor(GLFW.HRESIZE_CURSOR)
    g_MouseCursors[ImGuiMouseCursor_ResizeNESW+1] = GLFW.CreateStandardCursor(GLFW.ARROW_CURSOR)
    g_MouseCursors[ImGuiMouseCursor_ResizeNWSE+1] = GLFW.CreateStandardCursor(GLFW.ARROW_CURSOR)
    g_MouseCursors[ImGuiMouseCursor_Hand+1] = GLFW.CreateStandardCursor(GLFW.HAND_CURSOR)

    if install_callbacks
        GLFW.SetMouseButtonCallback(window, ImGui_ImplGlfw_MouseButtonCallback)
        GLFW.SetScrollCallback(window, ImGui_ImplGlfw_ScrollCallback)
        GLFW.SetKeyCallback(window, ImGui_ImplGlfw_KeyCallback)
        GLFW.SetCharCallback(window, ImGui_ImplGlfw_CharCallback)
    end

    g_ClientApi[] = client_api
    return true
end

function ImGui_ImplGlfw_InitForOpenGL(window::GLFW.Window, install_callbacks::Bool)
    ImGui_ImplGlfw_Init(window, install_callbacks, GlfwClientApi_OpenGL)
end

function ImGui_ImplGlfw_InitForVulkan(window::GLFW.Window, install_callbacks::Bool)
    ImGui_ImplGlfw_Init(window, install_callbacks, GlfwClientApi_Vulkan)
end

function ImGui_ImplGlfw_Shutdown()
    for cursor_n = 1:ImGuiMouseCursor_COUNT
        GLFW.DestroyCursor(g_MouseCursors[cursor_n])
        g_MouseCursors[cursor_n] = GLFW.Cursor(C_NULL)
    end
    g_ClientApi[] = GlfwClientApi_Unknown
end

function ImGui_ImplGlfw_UpdateMousePosAndButtons()
    # update buttons
    io = GetIO()
    for i = 1:length(g_MouseJustPressed)
        # if a mouse press event came, always pass it as "mouse held this frame",
        # so we don't miss click-release events that are shorter than 1 frame.
        mousedown = g_MouseJustPressed[i] || GLFW.GetMouseButton(g_Window[], GLFW.MouseButton(i-1))
        Set_MouseDown(io, i-1, mousedown)
        g_MouseJustPressed[i] = false
    end

    # update mouse position
    mouse_pos_backup = io.MousePos
    io.MousePos = ImVec2(-FLT_MAX, -FLT_MAX)
    if GLFW.GetWindowAttrib(g_Window[], GLFW.FOCUSED) != 0
        if io.WantSetMousePos
            GLFW.SetCursorPos(g_Window[], Cdouble(mouse_pos_backup.x), Cdouble(mouse_pos_backup.y))
        else
            mouse_x, mouse_y = GLFW.GetCursorPos(g_Window[])
            io.MousePos = ImVec2(Cfloat(mouse_x), Cfloat(mouse_y))
        end
    end
end

function ImGui_ImplGlfw_UpdateMouseCursor()
    io = GetIO()
    if (io.ConfigFlags & ImGuiConfigFlags_NoMouseCursorChange) != 0 ||
        GLFW.GetInputMode(g_Window[], GLFW.CURSOR) == GLFW.CURSOR_DISABLED
        return nothing
    end

    imgui_cursor = GetMouseCursor()
    if imgui_cursor == ImGuiMouseCursor_None || io.MouseDrawCursor
        # hide OS mouse cursor if imgui is drawing it or if it wants no cursor
        GLFW.SetInputMode(g_Window[], GLFW.CURSOR, GLFW.CURSOR_HIDDEN)
    else
        # show OS mouse cursor
        cursor = g_MouseCursors[imgui_cursor+1]
        GLFW.SetCursor(g_Window[], cursor.handle != C_NULL ? g_MouseCursors[imgui_cursor+1] : g_MouseCursors[ImGuiMouseCursor_Arrow+1])
        GLFW.SetInputMode(g_Window[], GLFW.CURSOR, GLFW.CURSOR_NORMAL)
    end

    return nothing
end

function ImGui_ImplGlfw_NewFrame(window::Union{GLFW.Window,Nothing}=nothing)
    window !== nothing && (g_Window[] = window;)
    io = GetIO()
    @assert ImFontAtlas_IsBuilt(io.Fonts) "Font atlas not built! It is generally built by the renderer back-end. Missing call to renderer _NewFrame() function? e.g. ImGui_ImplOpenGL3_NewFrame()."

    # setup display size (every frame to accommodate for window resizing)
    w, h = GLFW.GetWindowSize(g_Window[])
    display_w, display_h = GLFW.GetFramebufferSize(g_Window[])
    io.DisplaySize = ImVec2(Cfloat(w), Cfloat(h))
    w_scale = w > 0 ? Cfloat(display_w / w) : 0f0
    h_scale = h > 0 ? Cfloat(display_h / h) : 0f0
    io.DisplayFramebufferScale = ImVec2(w_scale, h_scale)

    # setup time step
    current_time = ccall((:glfwGetTime, GLFW.libglfw), Cdouble, ())
    io.DeltaTime = g_Time[] > 0.0 ? Cfloat(current_time - g_Time[]) : Cfloat(1.0/60.0)
    g_Time[] = current_time

    ImGui_ImplGlfw_UpdateMousePosAndButtons()
    ImGui_ImplGlfw_UpdateMouseCursor()

    # TODO: Gamepad navigation mapping
    # ImGui_ImplGlfw_UpdateGamepads()
end

# function ImGui_ImplGlfw_UpdateGamepads()
#     io = GetIO()
#     memset(io.NavInputs, 0, sizeof(io.NavInputs));
#     io.ConfigFlags & ImGuiConfigFlags_NavEnableGamepad == 0 && return
#
#     # update gamepad inputs
#     MAP_BUTTON(NAV_NO, BUTTON_NO) = buttons_count > BUTTON_NO && buttons[BUTTON_NO] == GLFW.PRESS && Set_NavInputs(io, NAV_NO, 1.0f0)
#     function MAP_ANALOG(NAV_NO, AXIS_NO, V0, V1)
#         v = (axes_count > AXIS_NO) ? axes[AXIS_NO] : V0
#         v = (v - V0) / (V1 - V0)
#         v > 1.0 && (v = 1.0)
#         Get_NavInputs(io, NAV_NO) < v && Set_NavInputs(io, NAV_NO, v)
#     end
#     axes_count = Cint(0)
#     buttons_count = Cint(0)
#     axes = @c GLFW.GetJoystickAxes(GLFW.JOYSTICK_1, &axes_count)
#     buttons = GLFW.GetJoystickButtons(GLFW.JOYSTICK_1, &buttons_count)
#     MAP_BUTTON(ImGuiNavInput_Activate,   0)     # Cross / A
#     MAP_BUTTON(ImGuiNavInput_Cancel,     1)     # Circle / B
#     MAP_BUTTON(ImGuiNavInput_Menu,       2)     # Square / X
#     MAP_BUTTON(ImGuiNavInput_Input,      3)     # Triangle / Y
#     MAP_BUTTON(ImGuiNavInput_DpadLeft,   13)    # D-Pad Left
#     MAP_BUTTON(ImGuiNavInput_DpadRight,  11)    # D-Pad Right
#     MAP_BUTTON(ImGuiNavInput_DpadUp,     10)    # D-Pad Up
#     MAP_BUTTON(ImGuiNavInput_DpadDown,   12)    # D-Pad Down
#     MAP_BUTTON(ImGuiNavInput_FocusPrev,  4)     # L1 / LB
#     MAP_BUTTON(ImGuiNavInput_FocusNext,  5)     # R1 / RB
#     MAP_BUTTON(ImGuiNavInput_TweakSlow,  4)     # L1 / LB
#     MAP_BUTTON(ImGuiNavInput_TweakFast,  5)     # R1 / RB
#     MAP_ANALOG(ImGuiNavInput_LStickLeft, 0,  -0.3,  -0.9)
#     MAP_ANALOG(ImGuiNavInput_LStickRight,0,  +0.3,  +0.9)
#     MAP_ANALOG(ImGuiNavInput_LStickUp,   1,  +0.3,  +0.9)
#     MAP_ANALOG(ImGuiNavInput_LStickDown, 1,  -0.3,  -0.9)
#     if axes_count > 0 && buttons_count > 0
#         io.BackendFlags |= ImGuiBackendFlags_HasGamepad
#     else
#         io.BackendFlags &= ~ImGuiBackendFlags_HasGamepad
#     end
# end


text:
#=
Please see
https://computationalmindset.com/en/neural-networks/ordinary-differential-equation-solvers.html#sys1
for details
=#

#=Original problem
x′ + x − y = 0
y′ − 4x + y = 0
x(0)=2
y(0)=0
=#

#=Esplicit form of equations
x′ = y - x
y′ = 4x - y
=#

#=Analytical solution
x(t) = e^t + e^−3t
y(t) = 2e^t − 2e^−3t
=#

using DifferentialEquations
using Plots

function ode_fn(du,u,p,t)
    x, y = u
    du[1] = y - x
    du[2] = 4.0 * x - y
end

an_sol_x(t) = exp(t) + exp(-3.0 * t)
an_sol_y(t) = 2.0 * exp(t) - 2.0 * exp(-3.0 * t)

t_begin=0.0
t_end=5.0
tspan = (t_begin,t_end)
x_init=2.0
y_init=0.0

prob = ODEProblem(ode_fn, [x_init, y_init], tspan)
num_sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
x_num_sol = [u[1] for u in num_sol.u]
y_num_sol = [u[2] for u in num_sol.u]

plot(num_sol.t, an_sol_x.(num_sol.t),
    linewidth=2, ls=:dash,
    title="System of 2 ODEs 1st order IVP solved by D.E. package",
    xaxis="t",
    label="analytical x",
    legend=true)
plot!(num_sol.t, an_sol_y.(num_sol.t),
    linewidth=2, ls=:dash,
    label="analytical y",
    legend=true)
plot!(num_sol.t, x_num_sol,
    linewidth=1,
    label="numerical x")
plot!(num_sol.t, y_num_sol,
    linewidth=1,
    label="numerical y")

max_stars_repo_name: ettoremessina/differential-equations

max_issues_repo_name: ettoremessina/differential-equations

max_forks_repo_name: ettoremessina/differential-equations


text:
export i_solve!, solve2

# do not wrap fields
NiLang.AD.GVar(param::AcousticPropagatorParams) = param
NiLang.AD.grad(param::AcousticPropagatorParams) = nothing

@i function i_one_step!(param::AcousticPropagatorParams, u, w, wold, φ, φ0, ψ, ψ0, c::AbstractMatrix{T}) where T
    @routine @invcheckoff begin
        @zeros Float64 Δt hx hy Δt2 Δt_hx Δt_hy Δt_hx2 Δt_hy2 Δt_2
        Δt += param.DELTAT
        hx += param.DELTAX
        hy += param.DELTAY
        Δt2 += Δt ^ 2
        Δt_hx += Δt / hx
        Δt_hy += Δt / hy
        Δt_hx2 += Δt_hx ^ 2
        Δt_hy2 += Δt_hy ^ 2
        Δt_2 += Δt/2
    end

    @invcheckoff @inbounds for j=2:param.NY+1
        for i=2:param.NX+1
            @routine begin
                @zeros T σpτ σpτΔt_2 cΔt_hx2 cΔt_hy2 dwx dwy dφx dψy στ anc1 anc2 anc3 uij
                σpτ += param.Σx[i,j] + param.Σy[i,j]
                σpτΔt_2 += σpτ * Δt_2
                cΔt_hx2 += Δt_hx2 * c[i,j]
                cΔt_hy2 += Δt_hy2 * c[i,j]
                dwx += w[i+1,j] + w[i-1,j]
                dwy += w[i,j+1] + w[i,j-1]
                dφx += φ0[i+1,j] - φ0[i-1,j]
                dψy += ψ0[i,j+1] - ψ0[i,j-1]
                στ += param.Σx[i,j] * param.Σy[i,j]

                anc1 += 2
                anc1 -= στ * Δt2
                anc1 -= 2 * cΔt_hx2
                anc1 -= 2 * cΔt_hy2
                anc2 += Δt_2 * Δt_hx
                anc3 += Δt_2 * Δt_hy
                uij += anc1 * w[i,j]
                uij += cΔt_hx2  *  dwx
                uij += cΔt_hy2  *  dwy
                uij += anc2 * dφx
                uij += anc3 * dψy
                uij -= wold[i,j]
                uij += σpτΔt_2 * wold[i,j]
                σpτΔt_2 += 1
            end
            u[i,j] += uij / σpτΔt_2
            ~@routine
        end
    end
    @invcheckoff @inbounds for j=2:param.NY+1
        for i=2:param.NX+1
            @routine begin
                @zeros T σmτ σmτ_2 dux duy cσmτ_2 σΔt τΔt anc1 anc2
                σmτ += param.Σx[i,j] - param.Σy[i,j]
                σmτ_2 += σmτ / 2
                dux += u[i+1,j] - u[i-1,j]
                duy += u[i,j+1] - u[i,j-1]
                cσmτ_2 += c[i,j] * σmτ_2
                σΔt += Δt * param.Σx[i,j]
                τΔt += Δt * param.Σy[i,j]
                σΔt -= 1
                τΔt -= 1
                anc1 += Δt_hx * cσmτ_2
                anc2 += Δt_hy * cσmτ_2
            end
            φ[i,j] -= σΔt * φ0[i,j]
            φ[i,j] -=  anc1 * dux
            ψ[i,j] -= τΔt * ψ0[i,j]
            ψ[i,j] += anc2 * duy
            ~@routine
        end
    end
    ~@routine
end

@i function i_solve!(param::AcousticPropagatorParams, src,
            srcv::Array{T, 1}, c::Array{T, 2},
            tu::Array{T,3}, tφ::Array{T,3}, tψ::Array{T,3}) where T

    @routine begin
        d2 ← zero(param.DELTAT)
        d2 += param.DELTAT^2
    end
    for i = 3:param.NSTEP+1
        i_one_step!(param, tu |> subarray(:,:,i), tu |> subarray(:,:,i-1), tu |> subarray(:,:,i-2),
            tφ |> subarray(:,:,i), tφ |> subarray(:,:,i-1), tψ |> subarray(:,:,i), tψ |> subarray(:,:,i-1), c)
        tu[src[1], src[2], i] += srcv[i-2]* d2
    end
    ~@routine
end

function solve2(param::AcousticPropagatorParams, srcloc,
            srcv::Array{T, 1}, c::Array{T, 2}) where T
    src = SeismicState(T, param.NX, param.NY)
    for i=1:param.NSTEP-1
        dest = SeismicState(T, param.NX, param.NY)
        src, = bennett_step!(dest, src, param, srcloc, srcv, c)
    end
    return src
end


text:
include(joinpath(pwd(), "src/models/particle.jl"))
include(joinpath(pwd(), "src/constraints/contact.jl"))

# Horizon
T = 51

# Time step
tf = 5.0
h = tf / (T-1)

# Bounds
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= 0.0
_ul = zeros(model.m)
ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
q1 = [0.0, 0.0, 1.0]
v1 = [3.0, 5.0, 0.0]
v2 = v1 - G_func(model,q1) * h
q2 = q1 + 0.5 * h * (v1 + v2)

x1 = [q1; q2]

xl, xu = state_bounds(model, T, x1 = x1)

# Objective
obj = PenaltyObjective(1000.0, model.m)

# Constraints
con_contact = contact_constraints(model, T)

# Problem
prob = problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con_contact
               )

# Trajectory initialization
X0 = [0.01 * rand(model.n) for t = 1:T] #linear_interp(x1, x1, T) # linear interpolation on state
U0 = [0.001 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
Z0 = pack(X0, U0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
@time Z̄ = solve(prob, copy(Z0), tol = 1.0e-5, c_tol = 1.0e-5)

check_slack(Z̄, prob)
X̄, Ū = unpack(Z̄, prob)

include(joinpath(pwd(), "src/models/visualize.jl"))
vis = Visualizer()
open(vis)
visualize!(vis, model, state_to_configuration(X̄), Δt = h)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
### A Pluto.jl notebook ###
# v0.12.19

using Markdown
using InteractiveUtils

# ╔═╡ a9d90294-61a6-11eb-1641-75e225e67315
begin
        using DataFrames
        using Plots
end

# ╔═╡ 61ca784e-5ceb-11eb-1423-5f9d48e72f26
md"# Lab 1 - Electric Field Visualization"

# ╔═╡ f4dbabba-61d4-11eb-1475-716ef6ea6109
md"
##### Info:

- Name: Omar Elrefaei
- ID: 20817426
- Date: 28/1/2021
- Code: [github](https://github.com/Omar-Elrefaei/small-projects/blob/master/ECE106/lab1-Electric_Field.jl)
- Made with: Julia Language & [Pluto Notebooks](https://github.com/fonsp/Pluto.jl/)

"

# ╔═╡ f688a336-61a3-11eb-1f9d-754da19d57ac
md"
## Question 1

![](https://i.imgur.com/nPFMTWX.png)
"

# ╔═╡ 456010ee-619e-11eb-0b5d-c5b88755b425
begin
        dataQ1 = [
                8.81 7.35 4.54 2.76 0.4 ;
                18.3 13.2 6.36 3.04 0.19]
        dataFrameQ1 = DataFrame(
                var"Before -1nC"        = dataQ1[1,:],
                var"After -1nC"         = dataQ1[2,:])
end

# ╔═╡ c5c0813e-61a2-11eb-3a8d-09f01538c980
md"
##### Answer:

When a sensor is in between the two charges: the magnitude of E-field at that sensor increases. But if the sensor is on the left/right of both charges, the value decreases.

When a sensor is (horizontally) between the two opposite charges, they both induce an E-Field in the same direction, and thus increase the magnitude. When they are in the same direction from the sensor, they cancel out instead."

# ╔═╡ 30bcb2f4-61a4-11eb-3765-a365cf2665b5
md"
## Question 2

##### Answer:

The E-field for the first 4 sensors point horizontally to the right. The vertical components of the field have been cancelled out due to both charges being opposite to each other and lying at equal horizontal distance from the sensors. On the other hand, the horizontal components compound and add to each other.

"

# ╔═╡ 0fd26dd4-61a6-11eb-23bb-bf99cdd17325
md"
## Question 3

##### Answer:
The field strength increase linearly proportional to the distance between the
two charges.

"

# ╔═╡ 719bfbe0-61b3-11eb-19ab-d9fab0dff0f0
begin
        dataQ3 = [ 200  400     600     800     1000; 0.19 0.44 0.78 1.21 1.49]
        plot(dataQ3[1,:], dataQ3[2,:], title="Plot 1", label="Real Data",
                lw=3, xlabel="Distance [m]", ylabel="E-field [V/m]")

        m = (dataQ3[2,end] - dataQ3[2,begin]) / (dataQ3[1,end] - dataQ3[1,begin])
        x = [i for i in 200:200:1000]
        y = -0.135 .+ x .* m
        plot!(x,y, label="Best fit", lw=2)
end

# ╔═╡ 06c52ba4-61bb-11eb-27a2-e9219195e2db
        dataFrameQ3 = DataFrame(var"Distance between charges" = dataQ3[1,:],
        var"Feild at remote sensor" = dataQ3[2,:])

# ╔═╡ 613aceaa-61a6-11eb-3ebb-47d2098a108e
md"
## Question 4

##### Answer:
It doubles from `0.19 V/m` to `0.38 V/m`

"

# ╔═╡ 738a33da-61b9-11eb-360f-aff73553fd68
md"
## Question 5

##### Answer:

$
\begin{aligned} E \displaystyle \propto dq
\end{aligned}$

"

# ╔═╡ ab45850e-5cf0-11eb-09de-cdd1a970c6e5
md"## Question 6"

# ╔═╡ 148bb47a-5cea-11eb-2c24-b702525f32e5
md"![](https://i.imgur.com/JRmAd9P.png)"

# ╔═╡ f3deac54-619c-11eb-164d-3b50118d91c5
data = [
        50              100             150     200     250     300     350;
        25.6    5.09    1.71    0.82    0.44    0.27    0.17]

# ╔═╡ 0edcf9f8-5cf3-11eb-3c35-db30a5695de9
dataFrame = DataFrame(Distance = data[1,:], Magnitude= data[2,:])

# ╔═╡ bdd8daf4-5cf0-11eb-1f64-9149e5e74e12
begin
        distances = data[1,:]; magnitudes = data[2,:]
        plot(distances, magnitudes, title="Plot 1", label="Real Data", lw=3)
        line2 = [1/x^2 for x in 1:7]*25
        plot!(distances, line2 , label="1/x^2")

        line2 = [1/x^3 for x in 1:7]*25
        plot!(distances, line2, label="1/x^3")
        plot!(xlabel="Distance [m]", ylabel="E-field [V/m]")

end

# ╔═╡ 9bf1fd0c-5cf6-11eb-10c7-19216015602d
md"### Answer:
the value of `n` is going to be between `2` and `3` as seen in the graph.
"

# ╔═╡ 0369de5e-5cfd-11eb-2421-534f52fe3db1
md"## Question 7
### Answer:

$
\begin{aligned}E_{total}=E_{1}+E_{2}=\dfrac{kq}{\left(r+d\right) ^{2}}-\dfrac{kq}{r^{2}}=\\
=\dfrac{kq}{\left( r+d\right) ^{2}}-\dfrac{kq}{r^{2}}=kq\left( \dfrac{r^{2}-\left( r+d\right) ^{2}}{r^{2}\left( r+d\right) ^{2}}\right) \\ =kq\left( \dfrac{r^{2}-r^{2}-2rd-d^{2}}{r^{2}\left( r+d\right) ^{2}}\right) \\ =kq\left( -d\right) \left( \dfrac{2r+d}{r^{2}\left( r+1\right) ^{2}}\right)\\

\begin{aligned}\because r >>d \ ; \ \therefore E_{tota1}\approx kq\left( \dfrac{-d}{r^{3}}\right)\end{aligned}\approx 0

\end{aligned}$

The equation matches confirms the observation in the graph. As the r increases, the field strength decays exponentially.
"

# ╔═╡ ad446368-5d20-11eb-3cad-c72fd941d953
md"
## Question 8

No, the observation would not change with a change of polarity

## Question 9
![](https://i.imgur.com/DvrqFgI.png)
The charge Q is in an **unstable equilibrium**, because at far away distance, when nudged there is no (neglegible) restoring force.

## Question 10
![](https://i.imgur.com/0CkOshH.png)
The charge Q is in an **unstable equilibrium**, because any small nudge vertically will throw it away. Interestingly, if charges were -ve, Q is stable vertically but unstable horizontally.

## Question 11
![](https://i.imgur.com/vHrfXYp.jpg)
The charge Q is in an **unstable equilibrium**, because any small nudge horizontally will throw it away.

## Question 12
![](https://i.imgur.com/YGXZ1C4.png)
The charge Q is in an **unstable equilibrium**, because any small nudge horizontally will throw it away. Again, one direction is stable with one polarity, and switches with the polarity changes.

"

# ╔═╡ Cell order:
# ╟─61ca784e-5ceb-11eb-1423-5f9d48e72f26
# ╟─a9d90294-61a6-11eb-1641-75e225e67315
# ╟─f4dbabba-61d4-11eb-1475-716ef6ea6109
# ╟─f688a336-61a3-11eb-1f9d-754da19d57ac
# ╠═456010ee-619e-11eb-0b5d-c5b88755b425
# ╟─c5c0813e-61a2-11eb-3a8d-09f01538c980
# ╟─30bcb2f4-61a4-11eb-3765-a365cf2665b5
# ╟─0fd26dd4-61a6-11eb-23bb-bf99cdd17325
# ╠═06c52ba4-61bb-11eb-27a2-e9219195e2db
# ╠═719bfbe0-61b3-11eb-19ab-d9fab0dff0f0
# ╟─613aceaa-61a6-11eb-3ebb-47d2098a108e
# ╟─738a33da-61b9-11eb-360f-aff73553fd68
# ╟─ab45850e-5cf0-11eb-09de-cdd1a970c6e5
# ╟─148bb47a-5cea-11eb-2c24-b702525f32e5
# ╠═f3deac54-619c-11eb-164d-3b50118d91c5
# ╠═0edcf9f8-5cf3-11eb-3c35-db30a5695de9
# ╠═bdd8daf4-5cf0-11eb-1f64-9149e5e74e12
# ╟─9bf1fd0c-5cf6-11eb-10c7-19216015602d
# ╟─0369de5e-5cfd-11eb-2421-534f52fe3db1
# ╟─ad446368-5d20-11eb-3cad-c72fd941d953

max_stars_repo_name: Omar-Elrefaei/small-projects

max_issues_repo_name: Omar-Elrefaei/small-projects

max_forks_repo_name: Omar-Elrefaei/small-projects


text:
#---------------------------------------------------------------------#
#This code computes the Equispaced points & weights
#Written by F.X. Giraldo on 4/2008
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function equispaced_points(P::Integer)

    #Initialize arrays
    xgl=zeros(P)
    wgl=zeros(P)

    #Construct array of points
    xgl=range(-1, length=P, stop=+1)

    return xgl,wgl
end #function