text:
"""
Module Elements --- elements of free modules.

A  `ModuleElt{K,V}`  represents  an  element  of  a free module where basis
elements  are of type  `K` and coefficients  of type `V`.  Usually you want
objects  of type `V` to be elements of  a ring, but it could also be useful
if  they just belong to  an abelian group. This  is similar to the SageMath
CombinatorialFreeModule.  You can  also see  them as  `SparseVector`s where
keys can be type `K` instead of integers.

This  basic  data  structure  is  used  in  my  packages  as  an  efficient
representation   at  many   places.  For   example,  the   `Monomial`  type
representing multivariate monomials is a `ModuleElt{Symbol,Int}`:

`x^2y^-3 ` is represented by `ModuleElt(:x=>2,:y=>-3)`

And  multivariate polynomials are  represented by a `ModuleElt{Monomial,C}`
where `C` is the type of the coefficients:

`x*y-z^2` is represented by ``ModuleElt(x*y=>1,z^2=>-1)

`ModuleElts`  are  also  used  for  cyclotomics, CycPols, elements of Hecke
algebras, etc…

A  `ModuleElt{K,V}` is essentially a  list of `Pairs{K,V}`. The constructor
takes  as argument a list of pairs, or a variable number of pair arguments,
or a generator of pairs.

We provide two implementations:

  - `HModuleElt`, an implementation by `Dict`s

This  requires  that  the  type  `K`  is  hashable.  It  is  a  very simple
implementation  since the interface of the type  is close to that of dicts;
the  only difference is weeding  out keys which have  a zero cofficient ---
which  is necessary since for testing equality of module elements one needs
a canonical form for each element.

-  a faster implementation  `ModuleElt` is obtained  by keeping the list of
pairs  sorted by key. This demands that the type `K` has a `isless` method.
This  implementation is two  to four times  faster than the  `Dict` one and
requires half the memory.

Both  implementations  have  the  same  methods,  which are mostly the same
methods  as a  `Dict` (`haskey`,  `getindex`, `setindex`, `keys`, `values`.
`pairs`,  `first`,  `iterate`,  `length`,  `eltype`), with some exceptions.
Adding  elements is implemented  as `merge(+,...)` which  is a variation on
`merge`  for  `Dict`s  where  keys  with  zero  value are deleted after the
operation (here `+` can be replaced by any operation `op` with the property
that `op(0,x)=op(x,0)=x`).

A module element can also be negated, or multiplied or divided (`/`or `//`)
by  some element (acting on coefficients)  if the method is defined between
type `V` and that element; there are also `zero` and `iszero` methods.

`ModuleElt`s  have  methods  `cmp`  and  `isless` which `HModuleElt`s don't
have. There is also `ModuleElts.merge2` which does the same as merge but is
valid  for more  general operations  (I use  it with  `min` and `max` which
implement `gcd` and `lcm` for `Monomial`s and `CycPol`s).

Here  is an example where basis elements are `Symbol`s and coefficients are
`Int`. As you can see in the examples, at the REPL (or in Jupyter or Pluto,
when  `IO`  has  the  `:limit`  attribute)  the  `show`  method  shows  the
coefficients  (bracketed  if  necessary,  which  is  when  they  have inner
occurences  of `+-*/`), followed by showing  the basis elements. The `repr`
method gives a representation which can be read back in julia:

```julia-repl
julia> a=ModuleElt(:xy=>1,:yx=>-1)
:xy-:yx

julia> repr(a)
"ModuleElt([:xy => 1, :yx => -1])"

julia> ModuleElt([:xy=>1//2,:yx=>-1])
(1//2):xy+(-1//1):yx
```

Setting  the  `IO`  property  `:showbasis`  to  a  custom printing function
changes how the basis elements are printed.

```julia-rep1
julia> show(IOContext(stdout,:showbasis=>(io,s)->string("<",s,">")),a)
3<xy>+2<yx>
```
We illustrate basic operations on `ModuleElt`s:

```julia-repl
julia> a-a
0

julia> a*99
99:xy-99:yx

julia> a//2
(1//2):xy+(-1//2):yx

julia> a/2
0.5:xy-0.5:yx

julia> a+ModuleElt(:yx=>1)
:xy

julia> a[:xy] # indexing by a basis element finds the coefficient
1

julia> a[:xx] # the coefficient of an absent basis element is zero.
0

julia> haskey(a,:xx)
false

julia> first(a)
:xy => 1

julia> collect(a)
2-element Vector{Pair{Symbol, Int64}}:
 :xy => 1
 :yx => -1

julia> collect(keys(a))
2-element Vector{Symbol}:
 :xy
 :yx

julia> collect(values(a))
2-element Vector{Int64}:
  1
 -1

julia> length(a)
2

julia> eltype(a)
Pair{Symbol, Int64}
```

In both implementations the constructor normalizes the constructed element,
removing zero coefficients and merging duplicate basis elements, adding the
corresponding   coefficients  (and   sorting  the   basis  in  the  default
implementation).  If  you  know  this  normalisation is unnecessary, to get
maximum  speed you can disable this  by giving the keyword `check=false` to
the constructor.

```julia-repl
julia> a=ModuleElt(:yy=>1, :yx=>2, :xy=>3, :yy=>-1;check=false)
:yy+2:yx+3:xy-:yy

julia> a=ModuleElt(:yy=>1, :yx=>2, :xy=>3, :yy=>-1)
3:xy+2:yx
```

Adding  or subtracting `ModuleElt`s does promotion  on the type of the keys
and the coefficients if needed:

```julia-repl
julia> a+ModuleElt([:z=>1.0])
3.0:xy+2.0:yx+1.0:z
```
"""
module ModuleElts

export ModuleElt, HModuleElt # data structure

#------------- implementation with Dicts ----------------------
struct HModuleElt{K,V}
  d::Dict{K,V} # the keys K should be hashable
  function HModuleElt(d::Dict{K,V};check::Bool=true)where {K,V}
    if check
      for (k,v) in d if iszero(v) delete!(d,k) end end
    end
    new{K,V}(d)
  end
end

function HModuleElt(x::Vector{Pair{K,V}};check=true) where{K,V}
  if !check || length(x)<=1 return HModuleElt(Dict(x)) end
  res=Dict{K,V}()
  for (k,v) in x
    if haskey(res,k) res[k]+=v
    else res[k]=v
    end
  end
  HModuleElt(res)
end

# forwarded methods
Base.haskey(x::HModuleElt,y...)=haskey(x.d,y...)
Base.keys(x::HModuleElt)=keys(x.d)
Base.values(x::HModuleElt)=values(x.d)

Base.getindex(x::HModuleElt{K,V},i) where{K,V}=haskey(x,i) ?  x.d[i] : zero(V)
Base.setindex!(x::HModuleElt{K,V},v,i) where{K,V}=
  if iszero(v) delete!(x.d,v) else x.d[i]=v end

# Valid for ops such that op(0,x)=x otherwise the result is wrong.
Base.merge(op::Function,a::HModuleElt,b::HModuleElt)=HModuleElt(merge(op,a.d,b.d))

@inline Base.:-(a::HModuleElt,b::HModuleElt)=a+(-b) # since merge2 not available
#-------------- faster implementation -------------------------------------
"""
`ModuleElt{K,V}`  has a  similar interface  to `Dict{K,V}`,  but instead of
assuming that `K` is hashable, it assumes that `K` is sortable. It also has
the advantage that ModuleElts are naturally sortable.

The  only  field,  a  `Vector{Pair{K,V}}`,  is  kept  sorted  by  `K`;  the
constructor  by  default  checks  sorting,  adds  values with same key, and
suppresses  keys  with  zero  value.  This  can  be bypassed by the keyword
`check=false`.
"""
struct ModuleElt{K,V}
  d::Vector{Pair{K,V}}
  function ModuleElt(d::AbstractVector{Pair{K,V}};check::Bool=true)where {K,V}
    if check normalize!(d) end
    new{K,V}(d)
  end
end

function normalize!(d::Vector{<:Pair})
  if isempty(d) return end
  if !issorted(d,by=first) sort!(d,by=first) end
  i=1
  ki,vi=d[i]
@inbounds for j in 2:length(d)
    if first(d[j])==ki
      vi+=last(d[j])
      d[i]=ki=>vi
    else
      if !iszero(vi) i+=1 end
      ki,vi=d[i]=d[j]
    end
  end
  resize!(d,iszero(vi) ? i-1 : i)
end

@inline Base.cmp(x::ModuleElt,y::ModuleElt)=cmp(x.d,y.d)
@inline Base.isless(x::ModuleElt,y::ModuleElt)=cmp(x,y)==-1

"""
`merge(op::Function,a::ModuleElt,b::ModuleElt)`

is  like merge(op,a,b) for Dicts, except that keys with value 0 are deleted
after the operation is done.

The  code is only  valid for `op`s  such that `op(0,x)=op(x,0)=x` otherwise
the result is wrong.
"""
function Base.merge(op::Function,a::ModuleElt,b::ModuleElt)
  (a,b)=promote(a,b)
  la=length(a.d)
  lb=length(b.d)
  res=similar(a.d,la+lb)
  ai=bi=1
  ri=0
@inbounds while ai<=la || bi<=lb
    if     ai>la res[ri+=1]=b.d[bi]; bi+=1
    elseif bi>lb res[ri+=1]=a.d[ai]; ai+=1
    else c=cmp(first(a.d[ai]),first(b.d[bi]))
      if     c>0 res[ri+=1]=b.d[bi]; bi+=1
      elseif c<0 res[ri+=1]=a.d[ai]; ai+=1
      else s=op(last(a.d[ai]),last(b.d[bi]))
        if !iszero(s) res[ri+=1]=first(a.d[ai])=>s end
        ai+=1; bi+=1
      end
    end
  end
  ModuleElt(resize!(res,ri);check=false)
end

"""
`merge2(op::Function,a::ModuleElt,b::ModuleElt)`

does  `op` between coefficients of  the same basis element  in `a` and `b`.
This  version works for  general ops (not  necessarily commutative or which
need not satisfy op(0,x)=op(x,0)=x), but has too much overhead currently to
replace  `merge` for  + or  other ops  such that  op(0,x)==op(x,0)=x. It is
useful for max or min which do lcm and gcd of `Monomial`s or `CycPol`s.
"""
function merge2(op::Function,a::ModuleElt,b::ModuleElt)
  (a,b)=promote(a,b)
  la=length(a.d)
  lb=length(b.d)
  res=similar(a.d,la+lb)
  ai=bi=1
  ri=0
@inbounds while ai<=la || bi<=lb
    if     ai>la
      be=b.d[bi]
      s=op(zero(last(be)),last(be))
      if !iszero(s) res[ri+=1]=first(be)=>s end
      bi+=1
    elseif bi>lb
      ae=a.d[ai]
      s=op(last(ae),zero(last(ae)))
      if !iszero(s) res[ri+=1]=first(ae)=>s end
      ai+=1
    else c=cmp(first(a.d[ai]),first(b.d[bi]))
      if     c>0
        be=b.d[bi]
        s=op(zero(last(be)),last(be))
        if !iszero(s) res[ri+=1]=first(be)=>s end
        bi+=1
      elseif c<0
        ae=a.d[ai]
        s=op(last(ae),zero(last(ae)))
        if !iszero(s) res[ri+=1]=first(ae)=>s end
        ai+=1
      else s=op(last(a.d[ai]),last(b.d[bi]))
        if !iszero(s) res[ri+=1]=first(a.d[ai])=>s end
        ai+=1; bi+=1
      end
    end
  end
  ModuleElt(resize!(res,ri);check=false)
end

Base.:-(a::ModuleElt,b::ModuleElt)=merge2(-,a,b) # faster than a+(-b)

"""
`getindex(x::ModuleElt,key)`  returns  the  value  associated  to the given
`key` in `x`. It returns zero if the key does not occur in `x`.
"""
function Base.getindex(x::ModuleElt{K,V},i) where {K,V}
  r=searchsortedfirst(x.d,i=>zero(V);by=first)
  if r>length(x.d) || first(x.d[r])!=i return zero(V) end
  last(x.d[r])
end

"""
`setindex!(x::ModuleElt,v,key)`  sets `v`  as the  value associated  to the
given `key` in `x`. Setting a value to zero or for a new key is expensive.
"""
function Base.setindex!(x::ModuleElt{K,V},v,key) where {K,V}
  r=searchsortedfirst(x.d,key=>zero(V);by=first)
  if r>length(x.d) || first(x.d[r])!=key
    if !iszero(v) splice!(x.d,r:r-1,[key=>v]) end
  elseif iszero(v) deleteat!(x.d,r)
  else x.d[r]=key=>v
  end
  v
end

function Base.haskey(x::ModuleElt{K,V},i) where {K,V}
  r=searchsortedfirst(x.d,i=>zero(V);by=first)
  r<=length(x.d) && first(x.d[r])==i
end

@inline Base.pairs(x::ModuleElt)=x.d
Base.keys(x::ModuleElt)=(first(p) for p in x.d)
Base.values(x::ModuleElt)=(last(p) for p in x.d)

# copied from Util.jl in order to have a self-contained file.
# Tries to determine which coefficients should be bracketed for unambiguous
# display
function format_coefficient(c::String)
  if c=="1" ""
  elseif c=="-1" "-"
  elseif match(r"^[-+]?([^-+*/]|√-|{-)*(\(.*\))?$",c)!==nothing c
  else "("*c*")"
  end
end

#-------------- methods which have same code in both implementations-------
# we assume that converting the keys from K to K1 does not change hashing
# for the Dict implementation or sorting for the other implementation
function Base.convert(::Type{ModuleElt{K,V}},a::ModuleElt{K1,V1}) where {K,K1,V,V1}
  if K==K1
    if V==V1 a
    elseif iszero(a) zero(ModuleElt{K,V})
    else ModuleElt(k=>convert(V,v) for (k,v) in a.d;check=false)
    end
  else
    if iszero(a) zero(ModuleElt{K,V})
    elseif V==V1  ModuleElt(convert(K,k)=>v for (k,v) in a.d;check=false)
    else ModuleElt(convert(K,k)=>convert(V,v) for (k,v) in a.d;check=false)
    end
  end
end

# I get a bug when putting this routine in the next loop. So here duplicated
function Base.convert(::Type{HModuleElt{K,V}},a::HModuleElt{K1,V1}) where {K,K1,V,V1}
  if K==K1
    if V==V1 a
    elseif iszero(a) zero(HModuleElt{K,V})
    else HModuleElt(k=>convert(V,v) for (k,v) in a.d;check=false)
    end
  else
    if iszero(a) zero(HModuleElt{K,V})
    elseif V==V1  HModuleElt(convert(K,k)=>v for (k,v) in a.d;check=false)
    else HModuleElt(convert(K,k)=>convert(V,v) for (k,v) in a.d;check=false)
    end
  end
end

for M in (:HModuleElt, :ModuleElt)
  @eval begin
@inline $M(x::Pair...;u...)=$M(collect(x);u...)
@inline $M(x::Base.Generator;u...)=$M(collect(x);u...)

@inline Base.:+(a::$M,b::$M)=merge(+,a,b)
Base.:-(a::$M)=iszero(a) ? a : $M(k=>-v for(k,v) in a;check=false)

# multiply module element by scalar
function Base.:*(a::$M{K,V},b)where {K,V}
  if iszero(b) || iszero(a)
    return zero($M{K,promote_type(V,typeof(b))})
  end
  let b=b
    $M(k=>v*b for (k,v) in a;check=false)
  end
end

@inline Base.iszero(x::$M)=isempty(x.d)
Base.zero(x::$M)=$M(empty(x.d))
Base.zero(::Type{$M{K,V}}) where{K,V}=$M(Pair{K,V}[])
# forwarded methods
@inline Base.:(==)(a::$M,b::$M)=a.d==b.d
@inline Base.first(x::$M)=first(x.d)
@inline Base.iterate(x::$M,y...)=iterate(x.d,y...)
@inline Base.length(x::$M)=length(x.d)
@inline Base.eltype(x::$M)=eltype(x.d)
@inline Base.hash(x::$M, h::UInt)=hash(x.d,h)

function Base.promote_rule(a::Type{$M{K1,V1}},
                           b::Type{$M{K2,V2}})where {K1,K2,V1,V2}
  $M{promote_type(K1,K2),promote_type(V1,V2)}
end

function Base.show(io::IO,m::$M)
  if !(get(io,:limit,false) || get(io,:TeX,false))
    print(io,Symbol($M),"(",m.d,")")
    return
  end
  if iszero(m) print(io,"0"); return end
  showbasis=get(io,:showbasis,nothing)
  if isnothing(showbasis)
    showbasis=(io,x)->repr(x;context=io)
  end
  start=true
  res=""
  for (k,v) in m
    v=repr(v;context=io)
    k=showbasis(io,k)
    if !isempty(k) v=format_coefficient(v) end
    if (isempty(v) || v[1]!='-') && !start v="+"*v end
    res*=v*k
    if start start=false end
  end
  if res=="" res="1" end
  print(io,res)
end

function Base.merge(f,m::$M{K,V},b;check=true)where {K,V}
  let b=b, f=f
    $M(k=>f(v,b) for (k,v) in m;check)
  end
end

function Base.merge(f,m::$M;check=true)
  $M(k=>f(v) for (k,v) in m;check)
end

end
end

for M in (:HModuleElt, :ModuleElt), op in (:/,:(//))
  @eval begin
    Base.$op(m::$M,b)=merge($op,m,b)
  end
end
end


text:
using Pkg
using JuMP;
using Plots;
using Distributions
using Random
using StatsBase

function create_serie()
    serie = zeros(120).+7
    d = Normal(0.0, 0.3)
    for t in 1:120
        serie[t] += cos(2*pi*t/12 + 6) + rand(d, 1)[1]
    end
    return serie
end

function create_dummy()
    A = zeros(120, 12)

    A[:,1] .= 1

    for i in 1:120
        col = mod(i,12)
        if col != 0
            A[i,col+1] = 1
        end
    end
    return A
end

function create_trig()
    A = zeros(120, 12)

    A[:,1] .= 1

    for t in 1:120, i = 2:12
        use_cos = (mod(i,2) == 0)

        j = (use_cos ? i/2 : i/2 - 0.5)
        if use_cos
            A[t,i] = cos(pi/6*j*t)
        else
            A[t,i] = sin(pi/6*j*t)
        end
    end
    return A
end

function solve_MQO(A,b)
    x = A'A \ A'b
    return x
end

b = create_serie()
Ad = create_dummy()
At = create_trig()
xd = solve_MQO(Ad,b)
xt = solve_MQO(At,b)

yd = Ad*xd
yt = At*xt

minimum(isapprox.(yd,yt))

p = plot(autocov(yd - b))

savefig(p,"Modelos Estruturais\\Lista 1\\imagens\\FAC.png")

ext: jl


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

pathToFMU = joinpath(dirname(@__FILE__), "..", "model", ENV["EXPORTINGTOOL"], "SpringPendulum1D.fmu")

# load FMU in temporary directory
myFMU = fmiLoad(pathToFMU)
@test isfile(myFMU.zipPath) == true
@test isdir(splitext(myFMU.zipPath)[1]) == true
fmiUnload(myFMU)

# load FMU in source directory
fmuDir = joinpath(splitpath(pathToFMU)[1:end-1]...)
myFMU = fmiLoad(pathToFMU; unpackPath=fmuDir)
@test isfile(splitext(pathToFMU)[1] * ".zip") == true
@test isdir(splitext(pathToFMU)[1]) == true

comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0

# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
@assert fmuStruct != nothing "Unknown fmuStruct, environment variable `FMUSTRUCT` = `$envFMUSTRUCT`"

t_start = 0.0
t_stop = 8.0

# test without recording values (but why?)
success, _ = fmiSimulate(fmuStruct, t_start, t_stop)
@test success

# test with recording values
success, savedValues = fmiSimulate(fmuStruct, t_start, t_stop; recordValues=["mass.s", "mass.v"], setup=true)
@test success
@test length(savedValues.saveval) == 100
@test length(savedValues.saveval[1]) == 2

t = savedValues.t
s = collect(d[1] for d in savedValues.saveval)
v = collect(d[2] for d in savedValues.saveval)
@test t[1] == t_start
@test t[end] == t_stop

# reference values from Simulation in Dymola2020x (Dassl)
@test s[1] == 0.5
@test v[1] == 0.0

if ENV["EXPORTINGTOOL"] == "Dymola/2020x" # ToDo: Linux FMU was corrupted
    @test abs(s[end] - 0.509219) < 0.01
    @test abs(v[end] - 0.314074) < 0.01
end

fmiUnload(myFMU)


text:
module TuringModels

const project_root = dirname(dirname(pathof(TuringModels)))

end # module


text:
function getImportantWordsInTopic(topic,dictionary,numToGet)
# gets the words in the topic with the heighest weights, and outputs the
# actual words in them and their associated weights
    topWords = Array{String}(numToGet, 1);
    weights = zeros(numToGet, 1);
    inds = reverse(sortperm(topic));
    topic = topic[inds];

    for i = 1:numToGet
        if topic[i] == 0
            topWords = topWords[1:i-1];
            weights = weights[1:i-1];
            break;
        end
        topWords[i] = dictionary[inds[i]];
        weights[i] = topic[i];
    end

    return topWords, weights;
end

function getImportantWordsInAllTopics(allTopics, dictionary, numTopics,numToGet)
    topWords = Array{String}(numToGet, numTopics);
    topicWordProbs = zeros(numToGet, numTopics);
    for i = 1:numTopics
        topWords[:,i], topicWordProbs[:,i] = getImportantWordsInTopic(allTopics[:,i], dictionary, numToGet);
    end
    return topWords,topicWordProbs;
end


text:
using TimeSeriesIO

using Base.Test

using TimeSeries

@testset "TimeArray to DataFrame" begin
    filename = joinpath(dirname(@__FILE__), "ford_2012.csv")
    ta2 = readtimearray(filename)
    ta2 = ta2["Open", "High", "Low", "Close"]
    df2 = DataFrame(ta2)
    @test names(df2) == [:Date, :Open, :High, :Low, :Close]

    ## With optional arguments
    df2 = DataFrame(ta2, timestamp=:DateTime, colnames=[:Open, :Close])
    @test names(df2) == [:DateTime, :Open, :Close]
end


text:
"""
	common constants
"""
const dt = 1.0                # [s]
const sensing_range = 1000.0  # [m]
const min_speed = 40.0        # [m/s]
const max_speed = 60.0        # [m/s]
const collision_range = 150.0 # [m]

const rho_dim   = 65
const theta_dim = 25
const phi_dim   = 25
const v_dim     = 3

const rho_min = 0.
const rho_max = 1800.0 # [m]

const theta_min = 0
const theta_max = 2 * pi

const phi_min = 0
const phi_max = 2 * pi

const v_min = min_speed
const v_max = max_speed

const discrete_rho = linspace(rho_min, rho_max, rho_dim)
# It is important to guarantee the appearance of π in angle discretizations!
const discrete_theta = linspace(theta_min, theta_max, theta_dim)
const discrete_phi   = linspace(phi_min, phi_max, phi_dim)
const discrete_v = linspace(v_min, v_max, v_dim)

"""
	sigma-point sampling
	non-cooperative
	increase noise, more robust MDP
"""
const sigma_v = 5.0
const sigma_alert_deg = 5.0
const sigma_coc_deg = 3.0
const sigma_theta_deg = 3.0
const sigma_phi_deg = 7.0

const sigma_alert = deg2rad(sigma_alert_deg)
const sigma_coc = deg2rad(sigma_coc_deg)
const sigma_theta = deg2rad(sigma_theta_deg)
const sigma_phi = deg2rad(sigma_phi_deg)

const sigma_str = "_Sigma_v_" * string(sigma_v) *
				  "_alert_" * string(sigma_alert_deg) *
				  "_theta_" * string(sigma_theta_deg) *
				  "_phi_" * string(sigma_phi_deg) *
				  "_coc_" * string(sigma_coc_deg)

# sigmas are tunable

# const sigmas = [0 sigma_alert -sigma_alert       0        0       0        0;
# 				0           0            0 sigma_v -sigma_v       0        0;
# 				0           0            0       0        0 sigma_v -sigma_v]

# const sigmas = [0 sigma_alert -sigma_alert         0          0       0        0       0        0;
# 				0           0            0 sigma_phi -sigma_phi       0        0       0        0;
# 				0           0            0         0          0 sigma_v -sigma_v       0        0;
# 				0           0            0         0          0       0        0 sigma_v -sigma_v]

const sigmas = [0 sigma_alert -sigma_alert         0          0           0            0       0        0       0        0;
				0           0            0 sigma_phi -sigma_phi           0            0       0        0       0        0;
				0           0            0         0          0 sigma_theta -sigma_theta       0        0       0        0;
				0           0            0         0          0           0            0 sigma_v -sigma_v       0        0;
				0           0            0         0          0           0            0       0        0 sigma_v -sigma_v]

const turn_rate_sigma_ind = 1
const phi_sigma_ind       = 2
const theta_sigma_ind     = 3
const v_own_sigma_ind     = 4
const v_int_sigma_ind     = 5

# const sigmas = [0 sigma_phi -sigma_phi           0            0;
# 				0         0          0 sigma_theta -sigma_theta]

# const phi_sigma_ind       = 1
# const theta_sigma_ind     = 2

const num_sigmas = size(sigmas, 2)

# need to increase the weight of the first sigma that executes the exact action
const actual_weight = 0.5;
const weights = [actual_weight; (1 - actual_weight) / (num_sigmas - 1) * ones(num_sigmas - 1)]

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
@testset "QuadEdgeTests" begin
  @testset "base" begin
    local baseIndex = rand(0:20) * 4 + 1
    @test DelaunayMeshes.base(baseIndex) == baseIndex
    @test DelaunayMeshes.base(baseIndex+1) == baseIndex
    @test DelaunayMeshes.base(baseIndex+2) == baseIndex
    @test DelaunayMeshes.base(baseIndex+3) == baseIndex
  end

  @testset "sym" begin
    local baseIndex = rand(0:20) * 4 + 1
    @test DelaunayMeshes.sym(baseIndex) == baseIndex + 2
    @test DelaunayMeshes.sym(baseIndex+2) == baseIndex
    @test DelaunayMeshes.sym(baseIndex+1) == baseIndex + 3
    @test DelaunayMeshes.sym(baseIndex+3) == baseIndex + 1
  end

  @testset "rot" begin
    local baseIndex = rand(0:20) * 4 + 1
    @test DelaunayMeshes.rot(baseIndex) == baseIndex + 1
    @test DelaunayMeshes.rot(baseIndex+1) == baseIndex + 2
    @test DelaunayMeshes.rot(baseIndex+2) == baseIndex + 3
    @test DelaunayMeshes.rot(baseIndex+3) == baseIndex
  end

  @testset "invrot" begin
    local baseIndex = rand(0:20) * 4 + 1
    @test DelaunayMeshes.invrot(DelaunayMeshes.rot(baseIndex)) == baseIndex
    @test DelaunayMeshes.invrot(DelaunayMeshes.rot(baseIndex+1)) == baseIndex + 1
    @test DelaunayMeshes.invrot(DelaunayMeshes.rot(baseIndex+2)) == baseIndex + 2
    @test DelaunayMeshes.invrot(DelaunayMeshes.rot(baseIndex+3)) == baseIndex + 3
  end

  @testset "makeedge!" begin
    # make empty subdivision and check if make edge adds edges
    local sd = DelaunayMeshes.SubDivision{Int}([],[],[],[])
    @test length(sd.edges)==0

    local ne = DelaunayMeshes.makeedge!(sd)
    @test ne == 1

    @test sd.edges[1].origin == 0
    @test sd.edges[1].next == 1
    @test sd.edges[1].constraintType == DelaunayMeshes.NoConstraint
    @test sd.edges[2].origin == 0
    @test sd.edges[2].next == 4
    @test sd.edges[2].constraintType == DelaunayMeshes.NoConstraint
    @test sd.edges[3].origin == 0
    @test sd.edges[3].next == 3
    @test sd.edges[3].constraintType == DelaunayMeshes.NoConstraint
    @test sd.edges[4].origin == 0
    @test sd.edges[4].next == 2
    @test sd.edges[4].constraintType == DelaunayMeshes.NoConstraint
  end

  @testset "splice!" begin
    local sd = DelaunayMeshes.SubDivision{Int}([],[],[],[])
    local ne1 = DelaunayMeshes.makeedge!(sd)
    local ne2 = DelaunayMeshes.makeedge!(sd)

    # splice edges together to form a single ring
    DelaunayMeshes.splice!(sd, ne1, ne2)
    @test DelaunayMeshes.onext(sd, ne1) == ne2
    @test DelaunayMeshes.onext(sd, ne2) == ne1
    @test DelaunayMeshes.onext(sd, ne1 + 1) == ne2 + 3
    @test DelaunayMeshes.onext(sd, ne1 + 3) == ne1 + 1
    @test DelaunayMeshes.onext(sd, ne2 + 1) == ne1 + 3
    @test DelaunayMeshes.onext(sd, ne2 + 3) == ne2 + 1

    # cut rings again
    DelaunayMeshes.splice!(sd, ne1, ne2)
    @test DelaunayMeshes.onext(sd, ne1) == ne1
    @test DelaunayMeshes.onext(sd, ne2) == ne2
    @test DelaunayMeshes.onext(sd, ne1 + 1) == ne1 + 3
    @test DelaunayMeshes.onext(sd, ne1 + 3) == ne1 + 1
    @test DelaunayMeshes.onext(sd, ne2 + 1) == ne2 + 3
    @test DelaunayMeshes.onext(sd, ne2 + 3) == ne2 + 1
  end

  @testset "endpoints" begin
    local sd = DelaunayMeshes.SubDivision{Float64}([],[],[],[])
    DelaunayMeshes.makeedge!(sd)
    sd.vertexCache = [-1,-1,-1]

    # check setting of end points works correctly
    DelaunayMeshes.endpoints!(sd, 1,1,2,3,4)
    @test DelaunayMeshes.org(sd, 1) == 1
    @test DelaunayMeshes.dest(sd, 1) == 2
    @test DelaunayMeshes.org(sd, DelaunayMeshes.sym(1)) == 2
    @test DelaunayMeshes.dest(sd, DelaunayMeshes.sym(1)) == 1
    @test sd.vertexCache[1] == 1
    @test sd.vertexCache[2] == 3
    @test DelaunayMeshes.lface(sd, 1) == 3
    @test DelaunayMeshes.rface(sd, 1) == 4
    @test DelaunayMeshes.lface(sd, DelaunayMeshes.sym(1)) == 4
    @test DelaunayMeshes.rface(sd, DelaunayMeshes.sym(1)) == 3

    # check end points works correctly for sym edge
    DelaunayMeshes.endpoints!(sd, 3,1,2,3,4)
    @test sd.edges[1].origin == 2
    @test sd.edges[3].origin == 1
    @test sd.vertexCache[1] == 3
    @test sd.vertexCache[2] == 1
    @test DelaunayMeshes.lface(sd, 3) == 3
    @test DelaunayMeshes.rface(sd, 3) == 4
    @test DelaunayMeshes.lface(sd, DelaunayMeshes.sym(3)) == 4
    @test DelaunayMeshes.rface(sd, DelaunayMeshes.sym(3)) == 3
  end

  @testset "subdivision" begin
    local sd = DelaunayMeshes.create_subdivision(7.0,3.0,5.0, 1.0, 2.0)
    @test isa(sd, DelaunayMeshes.SubDivision{Float64})
    testvertexcache(sd)

    # test that created triangle works
    local ea = 1
    @test DelaunayMeshes.org(sd, ea) == 1
    @test DelaunayMeshes.dest(sd, ea) == 2
    @test DelaunayMeshes.org(sd, DelaunayMeshes.rot(ea)) == 2
    @test DelaunayMeshes.dest(sd, DelaunayMeshes.rot(ea)) == 1

    local eb = DelaunayMeshes.lnext(sd, ea)
    @test DelaunayMeshes.org(sd, eb) == 2
    @test DelaunayMeshes.dest(sd, eb) == 3
    @test DelaunayMeshes.org(sd, DelaunayMeshes.rot(eb)) == 2
    @test DelaunayMeshes.dest(sd, DelaunayMeshes.rot(eb)) == 1

    local ec = DelaunayMeshes.lnext(sd, eb)
    @test DelaunayMeshes.org(sd, ec) == 3
    @test DelaunayMeshes.dest(sd, ec) == 1
    @test DelaunayMeshes.org(sd, DelaunayMeshes.rot(ec)) == 2
    @test DelaunayMeshes.dest(sd, DelaunayMeshes.rot(ec)) == 1
  end

  @testset "connect" begin
    local sd = DelaunayMeshes.create_subdivision(7.0,3.0,5.0, 1.0, 2.0)
    local ea = 1
    local eb = 5
    local ec = 9
    local innerFaceIndex = 1
    local outerFaceIndex = 2
    local newFaceIndex = 3

    # add single edge to main triangle starting at vertex 1
    push!(sd.vertices, 8.0)
    push!(sd.vertexCache, -1.)
    local ed = DelaunayMeshes.makeedge!(sd)
    DelaunayMeshes.endpoints!(sd, ed, 1, 4, innerFaceIndex, innerFaceIndex)
    DelaunayMeshes.splice!(sd, ed, ea)

    # and connect it to vertex 2
    local ee = DelaunayMeshes.connect!(sd, ea, DelaunayMeshes.sym(ed), 3.0)

    @test DelaunayMeshes.org(sd, DelaunayMeshes.rot(ea)) == outerFaceIndex
    @test DelaunayMeshes.dest(sd, DelaunayMeshes.rot(ea)) == innerFaceIndex
    @test DelaunayMeshes.org(sd, DelaunayMeshes.rot(eb)) == outerFaceIndex
    @test DelaunayMeshes.dest(sd, DelaunayMeshes.rot(eb)) == newFaceIndex
    @test DelaunayMeshes.org(sd, DelaunayMeshes.rot(ec)) == outerFaceIndex
    @test DelaunayMeshes.dest(sd, DelaunayMeshes.rot(ec)) == newFaceIndex
    @test DelaunayMeshes.org(sd, DelaunayMeshes.rot(ed)) == innerFaceIndex
    @test DelaunayMeshes.dest(sd, DelaunayMeshes.rot(ed)) == newFaceIndex
    @test DelaunayMeshes.org(sd, DelaunayMeshes.rot(ee)) == newFaceIndex
    @test DelaunayMeshes.dest(sd, DelaunayMeshes.rot(ee)) == innerFaceIndex
  end
end


text:
@testset "Reporting" begin
    model = initialise(lake_setup = lake_initial_state(Clear, Scheffer))

    @test nutrients(model) == model.lake.p.nutrients
    @test vegetation(model.lake) == [50.0]
    @test vegetation([model.lake.u[1]], model.lake.p) == [50.0]

    model = initialise(
        experiment = Experiment(nutrient_series = Dynamic()),
        lake_setup = lake_initial_state(Clear, Martin),
        municipalities = Dict(
            "main" => (
                Governance(
                    houseowner_type = Enforced(),
                    interventions = planner(plan(WastewaterTreatment)),
                ),
                100,
            ),
        ),
    )
    @test upgrade_efficiency(model) == 0
end


text:
Sphere(x...) = HyperSphere(x...)

Circle(x...) = HyperSphere(x...)

widths(c::HyperSphere{N, T}) where {N, T} = Vec{N, T}(radius(c)*2)
radius(c::HyperSphere) = c.r
origin(c::HyperSphere) = c.center
minimum(c::HyperSphere{N, T}) where {N, T} = Vec{N, T}(origin(c)) - Vec{N, T}(radius(c))
maximum(c::HyperSphere{N, T}) where {N, T} = Vec{N, T}(origin(c)) + Vec{N, T}(radius(c))


text:
struct DotProduct <: Distances.PreMetric
end

@inline function Distances._evaluate(::DotProduct,a::AbstractVector{T},b::AbstractVector{T}) where {T}
    @boundscheck if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    return dot(a,b)
end

@inline (dist::DotProduct)(a::AbstractArray,b::AbstractArray) = Distances._evaluate(dist,a,b)
@inline (dist::DotProduct)(a::Number,b::Number) = a*b


text:
module SfGeometry

using Unitful
import Elliptic: E, F

import Base: length

export Shape, Ellipsoid, TriaxialEllipsoid, Spheroid, Sphere, SphericalShell, Cylinder, Cuboid, Cube
export volume, radius, area, length
export unit_x, unit_y, unit_z, sphere_volume, sphere_radius, cylinder_volume, cylinder_radius, cylinder_length,
	spherical_cap_solid_angle, is_sphere, is_spheroid, is_triaxial, equatorial_radius, polar_radius,
	cross_sectional_area, spherical_shell_volume, spherical_shell_thickness, shape

const unit_x = [1, 0, 0]
const unit_y = [0, 1, 0]
const unit_z = [0, 0, 1]

sphere_volume(r::Unitful.Length) = (4π/3)r^3
sphere_area(r::Unitful.Length) = 4π * r^2
sphere_radius(v::Unitful.Volume) = cbrt(3v / 4π)

cylinder_volume(r::Unitful.Length, h::Unitful.Length) = π * r^2 * h
cylinder_radius(v::Unitful.Volume, h::Unitful.Length) = sqrt(v / (π * h))
cylinder_length(v::Unitful.Volume, r::Unitful.Length) = v / (π * r^2)

abstract type Shape end

shape(x::Shape) = x

abstract type Ellipsoid <: Shape end

struct TriaxialEllipsoid <: Ellipsoid
	a::Unitful.Length
	b::Unitful.Length
	c::Unitful.Length
end

struct Spheroid <: Ellipsoid
	equatorial_radius::Unitful.Length
	polar_radius::Unitful.Length
end

struct Sphere <: Ellipsoid
	radius::Unitful.Length
end

Sphere(v::Unitful.Volume) = Sphere(sphere_radius(v))

Ellipsoid(r::Unitful.Length) = Sphere(r)
Ellipsoid(equatorial_radius::Unitful.Length, polar_radius::Unitful.Length) =
	equatorial_radius == polar_radius ? Sphere(equatorial_radius) : Spheroid(equatorial_radius, polar_radius)

function Ellipsoid(a::Unitful.Length, b::Unitful.Length, c::Unitful.Length)
	if a == b && b == c
		return Sphere(a)
	elseif a == b || b == c
		return Spheroid(a, c)
	elseif a == c
		return Spheroid(a, b)
	else
		return TriaxialEllipsoid(a, b, c)
	end
end

struct Cylinder <: Shape
	radius::Unitful.Length
	length::Unitful.Length
end

struct Cuboid <: Shape
	length::Unitful.Length
	width::Unitful.Length
	height::Unitful.Length
end

Cube(l::Unitful.Length) = Cuboid(l, l, l)
Cube(v::Unitful.Volume) = Cube(cbrt(v))

struct SphericalShell <: Ellipsoid
	r_inner::Unitful.Length
	r_outer::Unitful.Length
end

SphericalShell(r_inner::Unitful.Length, v::Unitful.Volume) = SphericalShell(r_inner, r_inner + spherical_shell_thickness)
SphericalShell(s1::Sphere, s2::Sphere) = SphericalShell(min(radius(s1), radius(s2)), max(radius(s1), radius(s2)))

volume(x::TriaxialEllipsoid) = (4π/3) * x.a * x.b * x.c
volume(x::Spheroid) = (4π/3) * x.equatorial_radius^2 * x.polar_radius
volume(x::Sphere) = (4π/3) * x.radius^3
volume(x::SphericalShell) = spherical_shell_volume(x.r_inner, x.r_outer)
volume(x::Cylinder) = cylinder_volume(x.length, x.radius)
volume(x::Cuboid) = x.length * x.width * x.height

"""
    radius(x::TriaxialEllipsoid)

Returns the mean radius of `x`.
"""
radius(x::TriaxialEllipsoid) = (x.a + x.b + x.c) / 3

"""
    radius(x::Spheroid)

Returns the mean radius of `x`.
"""
radius(x::Spheroid) = (x.equatorial_radius + x.polar_radius) / 2

"""
    radius(x::Spheroid)

Radius of `x`.
"""
radius(x::Sphere) = x.radius

"""
    radius(x::Cylinder)

Returns the radius of `x`.
"""
radius(x::Cylinder) = x.radius

"""
    radius(x::SphericalShell)

The outer radius of `x`.
"""
radius(x::SphericalShell) = x.r_outer

"""
    equatorial_radius(x::Ellipsoid)

Not defined for triaxial TriaxialEllipsoids. Equivalent to the radius for spheres.
"""
equatorial_radius(x::Ellipsoid) = error("Equatorial radius not well defined for general TriaxialEllipsoids.")
equatorial_radius(x::Sphere) = x.radius
equatorial_radius(x::SphericalShell) = x.r_outer
equatorial_radius(x::Spheroid) = x.equatorial_radius

"""
    polar_radius(x::Ellipsoid)

Not defined for triaxial TriaxialEllipsoids. Equivalent to the radius for spheres.
"""
polar_radius(x::Ellipsoid) = error("Polar radius not well defined for general TriaxialEllipsoids.")
polar_radius(x::Sphere) = x.radius
polar_radius(x::SphericalShell) = x.r_outer
polar_radius(x::Spheroid) = x.polar_radius

area(x::Sphere) = 4π * x.radius^2
area(x::SphericalShell) = 4π * x.r_outer^2

function area(x::Spheroid)
	eq = equatorial_radius(x)
	po = polar_radius(x)

	if eq == po
		return sphere_area(eq)
	elseif eq > po
		# oblate
		e = sqrt(1 - po^2 / eq^2)
		return 2π * eq^2 * (1 + (po^2 / (e * eq^2)) * atanh(e))
	else
		# prolate
		e = sqrt(1 - eq^2 / po^2)
		return 2π * eq^2 * (1 + (po / (e * eq)) * asin(e))
	end
end

function area(x::TriaxialEllipsoid)
	# sort axes ascending
	c, b, a = sort([x.a, x.b, x.c])

	if a == b == c
		return sphere_area(a)
	elseif a == b
		return area(Spheroid(a, c))
	elseif a == c
		return area(Spheroid(a, b))
	elseif b == c
		return area(Spheroid(b, a))
	end

	cos_ϕ = c / a
	ϕ = acos(cos_ϕ)
	sin_ϕ = sin(ϕ)

	k = sqrt((a^2 * (b^2 - c^2))/(b^2 * (a^2 - c^2)))

	e = E(ϕ, k)
	f = F(ϕ, k)

	2π * c^2 + ((2π * a * b) / sin_ϕ) * (e * sin_ϕ^2 + f * cos_ϕ^2)
end

"""
    area(x::Cuboid)

Returns the surface area of `x`.
"""
area(x::Cuboid) = 2 * (x.width * x.length) + 2 * (x.width * x.height) + 2 * (x.height * x.length)

"""
    area(x::Cylinder)

Returns the surface area of `x`, including endcaps.
"""
area(x::Cylinder) = 2π * x.radius * x.length + 2π * x.radius^2

"""
	length(x::Cuboid)

Return the length of the longest side of `x`.

Note that this is not necessarily the same as `x.length`, as it is not enforced that `x.length` be `x`'s longest dimension.
"""
length(x::Cuboid) = max(x.length, x.width, x.height)

"""
    length(x::Cylinder)

Return the length of `x`.
"""
length(x::Cylinder) = x.length

"""
    length(x::Ellipsoid)

The diameter of `x` across its widest part.
"""
length(x::Ellipsoid) = max(equatorial_radius(x), polar_radius(x)) * 2
length(x::TriaxialEllipsoid) = max(x.a, x.b, x.c) * 2

"""
    cross_sectional_area(x::Ellipsoid)

The area of the largest cross-section of `x` in the coronal plane.

Not defined for triaxial TriaxialEllipsoids.
"""
cross_sectional_area(x::Ellipsoid) = equatorial_radius(x) * polar_radius(x) * π
cross_sectional_area(::TriaxialEllipsoid) = error("Cross section not well defined for triaxial TriaxialEllipsoids")

"""
    spherical_cap_solid_angle(θ)

Compute the solid angle of a cone with its apex at the apex of the solid angle and apex angle `2θ`.

See https://en.wikipedia.org/wiki/Solid_angle#Cone,_spherical_cap,_hemisphere
"""
spherical_cap_solid_angle(θ) = 2π * (1 - cos(θ)) * u"sr"

"""
    spherical_cap_solid_angle(h_cone::Unitful.Length, r_cone::Unitful.Length)

Compute the solid angle of a cone with its apex at the apex of the solid angle, height `h_cone` and radius `r_cone`.
"""
spherical_cap_solid_angle(h_cone::Unitful.Length, r_cone::Unitful.Length) = spherical_cap_solid_angle(atan(r_cone, h_cone))

"""
    spherical_shell_volume(r1::Unitful.Length, r2::Unitful.Length)

Return the volume of a spherical shell defined by an inner and outer radius.

Smaller of `r1`, `r2` is used as the inner radius, and the larger is used as the outer radius.
"""
spherical_shell_volume(r1::Unitful.Length, r2::Unitful.Length) = sphere_volume(max(r1, r2)) - sphere_volume(min(r1, r2))

"""
    spherical_shell_thickness(r_inner::Unitful.Length, v::Unitful.Volume)

Returns the thickness of a spherical shell with inner radius `r_inner` and volume `v`.
"""
spherical_shell_thickness(r_inner::Unitful.Length, v::Unitful.Volume) = cbrt((3v/4π) + r_inner^3) - r

end

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics


text:
module SolverTest

using Test, Random, Statistics, SparseArrays
import KLIEPInference
import KLIEPInference: chain, IsingSampler, Ψising, KLIEP, CD_KLIEP
import KLIEPInference: spKLIEP
using Distributions

include("alt_implementation.jl")

import .Alt

@testset "KLIEP low dim" begin
    Random.seed!(3)
    θx = chain(6; lenC=3, lb=0.2, ub=0.4, sgn=0)
    θy = copy(θx)

    nx = 500
    ny = 500
    X = rand(IsingSampler(θx; thin=1000), nx)
    Y = rand(IsingSampler(θy; thin=1000), ny)
    Ψx = Ψising(X)
    Ψy = Ψising(Y)

    scs_out = Alt.KLIEP(Ψx, Ψy, Alt.SCS_KLIEP())
    cd_out = KLIEP(Ψx, Ψy, CD_KLIEP())
    @test scs_out ≈ cd_out atol=0.001

    # mosek_out = Alt.KLIEP(Ψx, Ψy, Alt.Mosek_KLIEP())
    # @test scs_out ≈ mosek_out atol=0.01
end

@testset "spKLIEP delta == 0" begin
    Random.seed!(4)
    θx = chain(30; lenC=10, lb=0.1, ub=0.3, sgn=0)
    θy = copy(θx)

    nx = 300
    ny = 300
    X = rand(IsingSampler(θx; thin=1000), nx)
    Y = rand(IsingSampler(θy; thin=1000), ny)
    Ψx = Ψising(X)
    Ψy = Ψising(Y)

    p = length(θx)
    λ = 1.01 * quantile(Normal(), 1. - 0.05 / p) / sqrt(nx)

    θ = spKLIEP(Ψx, Ψy, λ, CD_KLIEP(); loadings=false)

    scs_out = Alt.spKLIEP(Ψx, Ψy, λ, Alt.SCS_KLIEP())
    @test scs_out ≈ θ atol=0.001

    # mosek_out = Alt.spKLIEP(Ψx, Ψy, λ, Alt.Mosek_KLIEP())
    # @test mosek_out ≈ θ atol=0.001
end

@testset "spKLIEP delta != 0" begin
    Random.seed!(5)
    θx = chain(30; lenC=10, lb=0.1, ub=0.3, sgn=0)
    θy = copy(θx)
    p = length(θx)
    Δ = sprand(p, 0.1)

    du = Uniform(0.1, 0.2)
    I, V = findnz(Δ)
    for i in I
        nv = rand(du)
        nv *= 2. * rand(Bernoulli()) - 1.
        θy[i] += nv
    end
    Δ = sparse(θx - θy)

    nx = 500
    ny = 500
    X = rand(IsingSampler(θx; thin=1000), nx)
    Y = rand(IsingSampler(θy; thin=1000), ny)
    Ψx = Ψising(X)
    Ψy = Ψising(Y)

    λ = 1.01 * quantile(Normal(), 1. - 0.05 / p) / sqrt(nx)

    θ = spKLIEP(Ψx, Ψy, λ, CD_KLIEP(); loadings=false)

    scs_out = Alt.spKLIEP(Ψx, Ψy, λ, Alt.SCS_KLIEP())
    @test scs_out ≈ θ atol=0.001

    # mosek_out = Alt.spKLIEP(Ψx, Ψy, λ, Alt.Mosek_KLIEP())
    # @test mosek_out ≈ θ atol=0.001
end

end


text:
#!/usr/bin/env julia

# using Memoize
# using LRUCache

const DATA_DIR   = "../19_data"
const DQUOTE = "\""

const INPUT_FILE = "input.txt"

# ----------------------------------------

function split_on_empty(xs)
    return split_on_pred(xs, x->isempty(x))
end

# The main logic behind the solution
function is_valid(rules, rule_ids, msg)
    # println("is_valid: rule_ids=$(rule_ids), msg=$(msg)")
    empties = [isempty(rule_ids), isempty(msg)]
    if any(empties)
         return all(empties)  # Better than filtering out on caller side
    end
    rule_first = rules[first(rule_ids)]
    if typeof(rule_first) == String
        # Recurse on the first char of msg
        return msg[1] == rule_first[1] && is_valid(rules, rule_ids[2:end], msg[2:end])
    else
        # Distribute the alternate cases of the first rule into required rule_ids in different calls
        return any(is_valid(rules, vcat(branch, rule_ids[2:end]), msg) for branch in rule_first)
    end
end

function rule_parse(rule_str)
    rule_id, rule_body = split(rule_str, ": ")
    rule_id = parse(Int, rule_id)
    # println("rule_parse: rule_id=$(rule_id), rule_body=$(rule_body)")
    if occursin(DQUOTE, rule_body)
        return (rule_id, string(rule_body[2:2]))
    else
        branches = [[parse(Int, opt) for opt in split_nonempty(branch)] for branch in split(rule_body, "|")]
        # println("branches=$(branches)")
        return (rule_id, branches)
    end
end

split_nonempty = s->filter(z->!isempty(z), split(s, " "))

# TODO: Is there a one-liner for this in Julia?
function split_on_pred(xs, pred)
    result = []
    current = []
    for x in xs
       if pred(x) && !isempty(current)
           push!(result, current)
           current = []
       else
           push!(current, x)
       end
    end
    if !isempty(current)
        push!(result, current)
    end
    return result
end

# ----------------------------------------

function main()
    lines::Array{String,1} = readlines("$DATA_DIR/$INPUT_FILE")
    rule_lines, msgs = split_on_empty(lines)
    rules = Dict([rule_parse(rule_str) for rule_str in rule_lines])
    valid_count = length(filter(b->b, [is_valid(rules, [0], msg) for msg in msgs]))
    println("Part 1: $(valid_count)")

    (_, rules[8])  = rule_parse("8: 42 | 42 8")
    (_, rules[11]) = rule_parse("11: 42 31 | 42 11 31")
    valid_count = length(filter(b->b, [is_valid(rules, [0], msg) for msg in msgs]))
    println("Part 2: $(valid_count)")
end

main()

max_stars_repo_name: jaycoskey/AdventOfCode

max_issues_repo_name: jaycoskey/AdventOfCode

max_forks_repo_name: jaycoskey/AdventOfCode


text:
#
#
#   These methods are for testing purposes only.
#   They are not optimized at all.
#
#

"""
Calculate the condition number from the given singular values
"""
function LinearAlgebra.cond(svs::AbstractVector)
    mi, ma = extrema(svs)
    ma / mi
end

"""
    calc_Bchain(B,N) -> (R, svs)

Calculate B^N as B*B*B*B*...*B from right to left.

Returns the result `R` and the singular values (column of `svs`) of all intermediate powers (rows of `svs`).
"""
function calc_Bchain(B, N)
    bs = size(B,1)
    R = Matrix{eltype(B)}(I, bs, bs)
    svs = zeros(Float64, N,bs)
    svc = 1
    for k in 1:N
        R = B * R
        svs[svc,:] = log.(svdvals(R))
        svc += 1
    end
    return (R, svs)
end

"""
Same as `calc_Bchain` but stabilizes the
matrix products by intermediate matrix decompositions.
"""
function calc_Bchain_svd(B, N; svdalg = svd)
    svs = zeros(Float64, N,size(B,1))
    F = svdalg(B)
    U, S, Vt = F.U, F.S, F.Vt
    svs[1,:] = log.(S)
    svc = 2

    for k in 2:N
        # multiply B from the left to USVt
        # and update USVt
        U = B * U
        U *= Diagonal(S)
        M = svdalg(U)
        U, S, Vtnew = M.U, M.S, M.Vt
        Vt = Vtnew * Vt

        # keep singular values
        svs[svc,:] = log.(S)
        svc += 1
    end

    return SVD(U, S, Vt), svs
end

"""
Same as `calc_Bchain` but stabilizes the
matrix products by intermediate matrix decompositions.
"""
function calc_Bchain_qr(B, N)
    svs = zeros(Float64, N,size(B,1))
    U, D, T = udt(B)
    svs[1,:] = log.(D)
    svc = 2

    for k in 2:N
        # multiply B from the left to UDT
        # and update UDT
        U = B * U
        U *= Diagonal(D)
        U, D, Tnew = udt(U)
        T = Tnew * T

        # keep singular values
        svs[svc,:] = log.(D)
        svc += 1
    end

    return UDT(U, D, T), svs
end

"""
Calculate fake "G(τ, 0)", i.e. [B^-N1 + B^N2]^-1
"""
function calc_tdgf_qr(B, N1, N2; inv_sum_method = inv_sum)
  if N1 != 0
    udtl = calc_Bchain_qr(inv(B), N1)[1]
  else
    udtl = udt!(Matrix(I*one(eltype(B)), size(B)))
  end

  if N2 != 0
    udtr = calc_Bchain_qr(B, N2)[1]
  else
    udtr = udt!(Matrix(I*one(eltype(B)), size(B)))
  end

  inv_sum_method(udtl, udtr)
end

"""
Calculate fake "G(τ, 0)", i.e. [B^-N1 + B^N2]^-1
"""
function calc_tdgf_svd(B, N1, N2; inv_sum_method = inv_sum, svdalg_inv = svd!, svdalg_chain = svd)
  # svdalg_nonmutating = Symbol(replace(string(nameof(svdalg)), "!" => ""))
  if N1 != 0
    svdl = calc_Bchain_svd(inv(B), N1; svdalg = svdalg_chain)[1]
  else
    svdl = svdalg(Matrix(I*one(eltype(B)), size(B)))
  end

  if N2 != 0
    svdr = calc_Bchain_svd(B, N2; svdalg = svdalg_chain)[1]
  else
    svdr = svdalg(Matrix(I*one(eltype(B)), size(B)))
  end

  inv_sum_method(svdl, svdr; svdalg = svdalg_inv)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "CRlibm"
version = v"1.0.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/JuliaIntervals/CRlibm.jl/raw/368ccdebd57e7a7cc96d47ab6b08e7a597495ff2/deps/src/crlibm-1.0beta4.tar.gz", "6836b4299f9421c99da2bdcd5e04a8d35577db4eb61161a401aa93751a96375d"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/crlibm*/

for p in ../patches/*.patch; do
    atomic_patch -p1 "${p}"
done

update_configure_scripts
export CFLAGS="-fPIC"
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -C scs_lib
make -j${nproc} libcrlibm.a
mkdir -p ${libdir}
cc -L. -shared -o ${libdir}/libcrlibm.${dlext} *.o scs_lib/*.o
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

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

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

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Testing

@testset "calculate cross-sectional properties" begin
    mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
    mesh = aster_read_mesh(mesh_file, "CYLINDER_20_TET4")
    # calculate cross-sectional properties A and Iₓ
    fixed1 = Problem(Dirichlet, "left support", 3, "displacement")
    fixed1.elements = create_elements(mesh, "FACE1")
    A = calculate_area(fixed1)
    info("cross-section area: $A")
    # real area is π
    @test isapprox(A, pi; rtol=0.1)
    Xc = calculate_center_of_mass(fixed1)
    info("center of mass: $Xc")
    @test isapprox(Xc, [0.0, 0.0, 0.0]; atol=1.0e-12)
    I = calculate_second_moment_of_mass(fixed1)
    info("moments:")
    info(I)
    I_expected = zeros(3, 3)
    I_expected[2,2] = I_expected[3,3] = pi/4
    rtol = norm(I[2,2]-I_expected[2,2]) / max(I[2,2],I_expected[2,2])
    info("I rtol = $rtol")
    @test isapprox(I, I_expected; rtol = 0.2)
end

#=
test subjects:
- calculate cross-sectional properties
- modal analysis with known solution

Fixed-fixed solution is ωᵢ = λᵢ²√(EI/ρA) , where λᵢ = cosh(λᵢℓ)cos(λᵢℓ)

1:  4.730040744862704
2:  7.853204624095838
3: 10.995607838001671

Youngs modulus is tuned such that lowest eigenfrequency matches 1.0

5 lowest eigenfrequencies using Code Aster and Tet4 elements:
numéro    fréquence (HZ)     norme d'erreur
    1       1.19789E+00        2.20137E-12
    2       1.20179E+00        1.99034E-12
    3       3.07391E+00        3.29226E-13
    4       3.08812E+00        2.91550E-13
    5       4.87370E+00        2.95986E-13

5 lowest eigenfrequencies using Code Aster and Tet10 elements:
numéro    fréquence (HZ)     norme d'erreur
    1       9.65942E-01        1.54950E-11
    2       9.66160E-01        1.62712E-11
    3       2.52127E+00        2.06544E-12
    4       2.52187E+00        1.77970E-12
    5       3.48584E+00        9.96170E-13

[1] De Silva, Clarence W. Vibration: fundamentals and practice. CRC press, 2006, p.355
=#

@testset "long rod natural frequencies" begin
    mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
    mesh = aster_read_mesh(mesh_file, "CYLINDER_20_TET10")
#   for (id, coords) in mesh.nodes
#       mesh.nodes[id][1] *= 5.0
#   end
    body = Problem(Elasticity, "rod", 3)
    body.elements = create_elements(mesh, "CYLINDER")
    #E = 50475.44814745859
    E = 50475.5
    rho = 1.0
    update!(body.elements, "youngs modulus", E)
    update!(body.elements, "poissons ratio", 0.3)
    update!(body.elements, "density", rho)
    # calculate cross-sectional properties A and Iₓ
    fixed1 = Problem(Dirichlet, "left support", 3, "displacement")
    fixed1.elements = create_elements(mesh, "FACE1")
    update!(fixed1.elements, "displacement 1", 0.0)
    update!(fixed1.elements, "displacement 2", 0.0)
    update!(fixed1.elements, "displacement 3", 0.0)
    fixed2 = Problem(Dirichlet, "right support", 3, "displacement")
    fixed2.elements = create_elements(mesh, "FACE2")
    update!(fixed2.elements, "displacement 1", 0.0)
    update!(fixed2.elements, "displacement 2", 0.0)
    update!(fixed2.elements, "displacement 3", 0.0)
    A = calculate_area(fixed1)
    info("cross-section area: $A")
    # using SALOME / SMESH, A = 2.82843
    # real area is π
    @test isapprox(A, pi; rtol=0.1)
    Xc = calculate_center_of_mass(fixed1)
    info("center of mass: $Xc")
    @test isapprox(Xc, [0.0, 0.0, 0.0]; atol=1.0e-5)
    I = calculate_second_moment_of_mass(fixed1)
    info("moments:")
    info(I)
    I_expected = zeros(3, 3)
    r = 1.0
    I_expected[2,2] = I_expected[3,3] = pi/4*r^2
    rtol = norm(I[2,2]-I_expected[2,2]) / max(I[2,2],I_expected[2,2])
    info("I rtol = $rtol")
    @test isapprox(I, I_expected; rtol = 0.2)
#=
    # apply transform Tx + b, in this case move cross-section to
    # xy-plane from yz-plane, i.e.
    # x₁ = y₂
    # y₁ = z₂
    T = [
        0.0 1.0 0.0
        0.0 0.0 1.0]
    b = [0.0, 0.0]
    X 1 = first(cross_section)("geometry", [1/3, 1/3], 0.0)
    apply_affine_transform!(cross_section, T, b)
    X2 = first(cross_section)("geometry", [1/3, 1/3], 0.0)
    info("X1 = $X1, X2 = $X2")
    @test isapprox(T*X1+b, X2)
=#
    c = sqrt(E*I[2,2]/(rho*A))
    info("c = $c")

    # analytical solution is
    l = 20.0
    r = 1.0
    la = 4.730040744862704/l

    # semi-analytical (c numerical)
    freq_sa = (c*la^2)/(2*pi)
    info("freq_sa = $freq_sa")

    A = pi*r^2
    I = pi/4*r^4
    c = sqrt(E*I/(rho*A))
    info("c analytical = $c")
    freq_a = (c*la^2)/(2*pi)

    info("freq_a = $freq_a")

    solver = Solver(Modal, body, fixed1, fixed2)
    solver.properties.nev = 5
    solver()
    freqs_jf = sqrt(solver.properties.eigvals)/(2.0*pi)
    # with Tet4 elements
    #freqs_ca = [1.19789E+00, 1.20179E+00, 3.07391E+00, 3.08813E+00, 4.87370E+00]
    # with Tet10 elements
    freqs_ca = [9.65942E-01, 9.66160E-01, 2.52127E+00, 2.52187E+00, 3.48584E+00]

    # looks that juliafem results are more close to 1.0, maybe different integration order
    rtol1 = norm(freq_sa - freqs_jf[1])/max(freq_sa, freqs_jf[1])
    rtol2 = norm(freq_a - freqs_jf[1])/max(freq_a, freqs_jf[1])
    info("rtol 1 = $rtol1, rtol 2 = $rtol2")
    passed = true
    for (f1, f2) in zip(freqs_jf, freqs_ca)
        rtol = norm(f1-f2) / max(f1,f2)
        @printf "JF: %8.5e | CA: %8.5e | rtol: %8.5e\n" f1 f2 rtol
        passed &= (rtol < 3.0e-2)
    end
    @test rtol2 < 3.0e-2
    @test passed

#=
    result = XDMF()
    for (i, freq) in enumerate(freqs)
        isapprox(freq, 0.0) && continue
        info("$i freq: $freq")
        xdmf_new_result!(result, body, freq)
        xdmf_save_field!(result, body, freq, "displacement"; field_type="Vector")
    end
    xdmf_save!(result, "/tmp/rod_nf.xmf")
=#

end

@testset "eigenvalues of cube (tet4)" begin
    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
    mesh = aster_read_mesh(meshfile, "CUBE_TET4")
    cube = Problem(mesh, Elasticity, "CUBE", 3)
    update!(cube.elements, "youngs modulus", 10000.0)
    update!(cube.elements, "poissons ratio", 0.3)
    update!(cube.elements, "density", 10.0)
    sym23 = create_elements(mesh, "FACE231")
    update!(sym23, "displacement 1", 0.0)
    sym13 = create_elements(mesh, "FACE131")
    update!(sym13, "displacement 2", 0.0)
    sym12 = create_elements(mesh, "FACE121")
    update!(sym12, "displacement 3", 0.0)
    bcs = Problem(Dirichlet, "bcs", 3, "displacement")
    bcs.elements = [sym23; sym13; sym12]
    solver = Solver(Modal)
    solver.properties.nev = 5
    push!(solver, cube, bcs)
    solver()
    freqs_jf = sqrt(solver.properties.eigvals)/(2.0*pi)
    freqs_ca = [3.73724E+00, 3.73724E+00, 4.93519E+00, 6.59406E+00, 7.65105E+00]
    for (f1, f2) in zip(freqs_jf, freqs_ca)
        rtol = norm(f1-f2) / max(f1,f2)
        @printf "JF: %8.5e | CA: %8.5e | rtol: %8.5e\n" f1 f2 rtol
        @test rtol < 1.0e-5
    end
end


text:
using ModelingToolkit, DifferentialEquations
@variables t
D = Differential(t)
sts = @variables x(t) RHS(t)  # RHS is observed
ps = @parameters τ=0.3      # parameters
@named m = ODESystem([ RHS  ~ (1 - x)/τ, D(x) ~ RHS ], t, sts, ps)
ms = structural_simplify(m)
prob = ODEProblem(ms, [x => 1.1], (0.0,1.0), [])
sol = solve(prob, Tsit5())

using BenchmarkTools
const constant = 10
# actually constant can be reassigned by the same type, only issues warnings
add(x) = constant + x
@btime add(3)

var = 10
add2(x) = var::Int + x
@btime add2(3)

add3(x, y=var) = var +x
@btime add3(3, $var)

const semi = Ref(10)
add4(x) = x + semi[]
@btime add4(3)

using ModelingToolkit, DifferentialEquations, LabelledArrays
using MTKHelpers
function samplesystem(;name,τ = 3.0, p1=1.1, p2=1.2)
    @variables t
    D = Differential(t)
    sts = @variables x(t) RHS(t)             # RHS is observed
    ps = @parameters τ=τ p1=p1 p2=p2       # parameters
    ODESystem([ RHS  ~ p1/p2 * (1 - x)/τ, D(x) ~ RHS ], t, sts, ps; name)
end
@named m = samplesystem()
@named sys = embed_system(m)
prob = ODEProblem(sys, [m.x => 0.0], (0.0,10.0))

# setup position matching#
popt = SLVector(m₊x=0.1, m₊p1=2.1)
ps = ProblemParSetter(sys, keys(popt))

# extract optimized
get_paropt_labeled(ps, prob)

# update states and parameters
prob2 = update_statepar(ps, popt, prob)
get_paropt_labeled(ps, prob2) == popt

# Documentation example
using ModelingToolkit, DifferentialEquations, ComponentArrays
using MTKHelpers
function samplesystem(;name,τ = 3.0, p1=1.1, p2=1.2)
    @variables t
    D = Differential(t)
    sts = @variables x(t) RHS(t)             # RHS is observed
    ps = @parameters τ=τ p1=p1 p2=p2       # parameters
    ODESystem([ RHS  ~ p1/p2 * (1 - x)/τ, D(x) ~ RHS ], t, sts, ps; name)
end
@named m = samplesystem()
@named sys = embed_system(m)
prob = ODEProblem(sys, [m.x => 0.0], (0.0,10.0))

# setup position matching
popt = ComponentVector(m₊x=0.1, m₊p1=2.1)
#ps = ProblemParSetter(sys, first(getaxes(popt))) # better type-stability
ps = ProblemParSetter(sys, popt)

# extract optimized
get_paropt_labeled(ps, prob)

# update states and parameters
prob2 = update_statepar(ps, popt, prob)
prob2.p # p is still a plain vector
label_par(ps, prob2.p).m₊p1 == popt.m₊p1 # attach labels and access properties
label_state(ps, prob2.u0).m₊x == popt.m₊x # attach labels and access properties
get_paropt_labeled(ps, prob2) == popt


text:
module DSP

using FFTW
using LinearAlgebra: mul!, rmul!
using IterTools: subsets

export conv, deconv, filt, filt!, xcorr

# This function has methods added in `peridograms` but is not exported,
# so we define it here so one can do `DSP.allocate_output` instead of
# `DSP.Periodograms.allocate_output`.
function allocate_output end

include("dspbase.jl")

include("util.jl")
include("unwrap.jl")
include("windows.jl")
include("periodograms.jl")
include("Filters/Filters.jl")
include("lpc.jl")
include("estimation.jl")

using Reexport
@reexport using .Util, .Windows, .Periodograms, .Filters, .LPC, .Unwrap, .Estimation

include("deprecated.jl")
end


text:
"""
    n_random_modules(n::Int64)

This returns *a function* which, when applied to a network, will randomly assign
every species to one of `n` modules. The correct way to apply this function to a
network `N` is, therefore `n_random_modules(4)(N)` (with four modules).
"""
function n_random_modules(n::Int64)
  n > 0 || throw(ArgumentError("The number of modules must be a positive integer, you used $(n)"))
  n_random = (N::AbstractEcologicalNetwork) -> (N, Dict([species(N)[i] => rand(1:n) for i in 1:richness(N)]))
  return n_random
end

"""
    each_species_its_module(N::T) where {T<:AbstractEcologicalNetwork}

Returns a dictionary in which each species is its own module. This is used as a
starting point for `lp` and `salp` internally. This is often a very poor
starting point for `brim`, and should probably not be used on its own.
"""
function each_species_its_module(N::T) where {T<:AbstractEcologicalNetwork}
  L = Dict([species(N)[i]=>i for i in 1:richness(N)])
  return (N, L)
end


text:
# Latafat, Patrinos, "Asymmetric forward–backward–adjoint splitting for
# solving monotone inclusions involving three operators", Computational
# Optimization and Applications, vol. 68, no. 1, pp. 57-93 (2017).
#
# Latafat, Patrinos, "Primal-dual proximal algorithms for structured convex
# optimization: a unifying framework", In Large-Scale and Distributed
# Optimization, Giselsson and Rantzer, Eds. Springer International Publishing,
# pp. 97–120 (2018).
#
# Chambolle, Pock, "A First-Order Primal-Dual Algorithm for Convex Problems
# with Applications to Imaging", Journal of Mathematical Imaging and Vision,
# vol. 40, no. 1, pp. 120-145 (2011).
#
# Condat, "A primal–dual splitting method for convex optimization
# involving Lipschitzian, proximable and linear composite terms",
# Journal of Optimization Theory and Applications, vol. 158, no. 2,
# pp 460-479 (2013).
#
# Vũ, "A splitting algorithm for dual monotone inclusions involving
# cocoercive operators", Advances in Computational Mathematics, vol. 38, no. 3,
# pp. 667-681 (2013).

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalCore: Zero, IndZero, convex_conjugate
using LinearAlgebra
using Printf

"""
    AFBAIteration(; <keyword-arguments>)

Iterator implementing the asymmetric forward-backward-adjoint algorithm (AFBA, see [1]).

This iterator solves convex optimization problems of the form

    minimize f(x) + g(x) + (h □ l)(L x),

where `f` is smooth, `g` and `h` are possibly nonsmooth and `l` is strongly
convex. Symbol `□` denotes the infimal convolution, and `L` is a linear mapping.

Points `x0` and `y0` are the initial primal and dual iterates, respectively.
If unspecified, functions `f`, `g`, and `h` default to the identically zero
function, `l` defaults to the indicator of the set `{0}`, and `L` defaults to
the identity. Important keyword arguments, in case `f` and `l` are set, are
the Lipschitz constants `beta_f` and `beta_l` (see below).

The iterator implements Algorithm 3 of [1] with constant stepsize (α_n=λ)
for several prominant special cases:
1) θ = 2          ==>   Corresponds to the Vu-Condat Algorithm [3, 4].
2) θ = 1, μ=1
3) θ = 0, μ=1
4) θ ∈ [0,∞), μ=0

See [2, Section 5.2] and [1, Figure 1] for stepsize conditions, special cases,
and relation to other algorithms.

See also: [`AFBA`](@ref).

# Arguments
- `x0`: initial primal point.
- `y0`: initial dual point.
- `f=Zero()`: smooth objective term.
- `g=Zero()`: proximable objective term.
- `h=Zero()`: proximable objective term.
- `l=IndZero()`: strongly convex function.
- `L=I`: linear operator (e.g. a matrix).
- `beta_f=0`: Lipschitz constant of the gradient of `f`.
- `beta_l=0`: Lipschitz constant of the gradient of `l` conjugate.
- `theta=1`: nonnegative algorithm parameter.
- `mu=1`: algorithm parameter in the range [0,1].
- `gamma1`: primal stepsize (see [1] for the default choice).
- `gamma2`: dual stepsize (see [1] for the default choice).

# References
1. Latafat, Patrinos, "Asymmetric forward-backward-adjoint splitting for solving monotone inclusions involving three operators", Computational Optimization and Applications, vol. 68, no. 1, pp. 57-93 (2017).
2. Latafat, Patrinos, "Primal-dual proximal algorithms for structured convex optimization: a unifying framework", In Large-Scale and Distributed Optimization, Giselsson and Rantzer, Eds. Springer International Publishing, pp. 97-120 (2018).
3. Condat, "A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms", Journal of Optimization Theory and Applications, vol. 158, no. 2, pp 460-479 (2013).
4. Vũ, "A splitting algorithm for dual monotone inclusions involving cocoercive operators", Advances in Computational Mathematics, vol. 38, no. 3, pp. 667-681 (2013).
"""
Base.@kwdef struct AFBAIteration{R,Tx,Ty,Tf,Tg,Th,Tl,TL,Tbetaf,Tbetal,Ttheta,Tmu,Tlambda}
    f::Tf = Zero()
    g::Tg = Zero()
    h::Th = Zero()
    l::Tl = IndZero()
    L::TL = if isa(h, Zero) L = 0 * I else I end
    x0::Tx
    y0::Ty
    beta_f::Tbetaf = if isa(f, Zero)
        real(eltype(x0))(0)
    else
        error("argument beta_f must be specified together with f")
    end
    beta_l::Tbetal = if isa(l, IndZero)
        real(eltype(x0))(0)
    else
        error("argument beta_l must be specified together with l")
    end
    theta::Ttheta = real(eltype(x0))(1)
    mu::Tmu = real(eltype(x0))(1)
    lambda::Tlambda = real(eltype(x0))(1)
    gamma::Tuple{R, R} = if lambda != 1
        error("if lambda != 1, then you need to provide stepsizes manually")
    else
        T = real(eltype(x0))
        AFBA_default_stepsizes(L, h, T(theta), T(mu), T(beta_f), T(beta_l))
    end
end

Base.IteratorSize(::Type{<:AFBAIteration}) = Base.IsInfinite()

"""
    VuCondatIteration(; <keyword-arguments>)

Iterator implementing the Vũ-Condat primal-dual algorithm [1, 2].

This iterator solves convex optimization problems of the form

    minimize f(x) + g(x) + (h □ l)(L x),

where `f` is smooth, `g` and `h` are possibly nonsmooth and `l` is strongly
convex. Symbol `□` denotes the infimal convolution, and `L` is a linear mapping.

This iteration is equivalent to [`AFBAIteration`](@ref) with `theta=2`;
for all other arguments see [`AFBAIteration`](@ref).

See also: [`AFBAIteration`](@ref), [`VuCondat`](@ref).

# References
1. Condat, "A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms", Journal of Optimization Theory and Applications, vol. 158, no. 2, pp 460-479 (2013).
2. Vũ, "A splitting algorithm for dual monotone inclusions involving cocoercive operators", Advances in Computational Mathematics, vol. 38, no. 3, pp. 667-681 (2013).
"""
VuCondatIteration(; kwargs...) = AFBAIteration(kwargs..., theta=2)

"""
    ChambollePockIteration(; <keyword-arguments>)

Iterator implementing the Chambolle-Pock primal-dual algorithm [1].

This iterator solves convex optimization problems of the form

    minimize g(x) + h(L x),

where `g` and `h` are possibly nonsmooth, and `L` is a linear mapping.

See also: [`AFBAIteration`](@ref), [`ChambollePock`](@ref).

This iteration is equivalent to [`AFBAIteration`](@ref) with `theta=2`, `f=Zero()`, `l=IndZero()`;
for all other arguments see [`AFBAIteration`](@ref).

# References
1. Chambolle, Pock, "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging", Journal of Mathematical Imaging and Vision, vol. 40, no. 1, pp. 120-145 (2011).
"""
ChambollePockIteration(; kwargs...) = AFBAIteration(kwargs..., theta=2, f=Zero(), l=IndZero())

Base.@kwdef struct AFBAState{Tx,Ty}
    x::Tx
    y::Ty
    xbar::Tx = similar(x)
    ybar::Ty = similar(y)
    gradf::Tx = similar(x)
    gradl::Ty = similar(y)
    FPR_x::Tx = similar(x)
    FPR_y::Ty = similar(y)
    temp_x::Tx = similar(x)
    temp_y::Ty = similar(y)
end

function Base.iterate(iter::AFBAIteration, state::AFBAState = AFBAState(x=copy(iter.x0), y=copy(iter.y0)))
    # perform xbar-update step
    gradient!(state.gradf, iter.f, state.x)
    mul!(state.temp_x, iter.L', state.y)
    state.temp_x .+= state.gradf
    state.temp_x .*= -iter.gamma[1]
    state.temp_x .+= state.x
    prox!(state.xbar, iter.g, state.temp_x, iter.gamma[1])

    # perform ybar-update step
    gradient!(state.gradl, convex_conjugate(iter.l), state.y)
    state.temp_x .= iter.theta .* state.xbar .+ (1 - iter.theta) .* state.x
    mul!(state.temp_y, iter.L, state.temp_x)
    state.temp_y .-= state.gradl
    state.temp_y .*= iter.gamma[2]
    state.temp_y .+= state.y
    prox!(state.ybar, convex_conjugate(iter.h), state.temp_y, iter.gamma[2])

    # the residues
    state.FPR_x .= state.xbar .- state.x
    state.FPR_y .= state.ybar .- state.y

    # perform x-update step
    state.temp_y .= (iter.mu * (2 - iter.theta) * iter.gamma[1]) .* state.FPR_y
    mul!(state.temp_x, iter.L', state.temp_y)
    state.x .+= iter.lambda .* (state.FPR_x .- state.temp_x)

    # perform y-update step
    state.temp_x .= ((1 - iter.mu) * (2 - iter.theta) * iter.gamma[2]) .* state.FPR_x
    mul!(state.temp_y, iter.L, state.temp_x)
    state.y .+= iter.lambda .* (state.FPR_y .+ state.temp_y)

    return state, state
end

default_stopping_criterion(tol, ::AFBAIteration, state::AFBAState) = norm(state.FPR_x, Inf) + norm(state.FPR_y, Inf) <= tol
default_solution(::AFBAIteration, state::AFBAState) = (state.xbar, state.ybar)
default_display(it, ::AFBAIteration, state::AFBAState) = @printf("%6d | %7.4e\n", it, norm(state.FPR_x, Inf) + norm(state.FPR_y, Inf))

"""
    AFBA(; <keyword-arguments>)

Constructs the asymmetric forward-backward-adjoint algorithm (AFBA, see [1]).

This algorithm solves convex optimization problems of the form

    minimize f(x) + g(x) + (h □ l)(L x),

where `f` is smooth, `g` and `h` are possibly nonsmooth and `l` is strongly
convex. Symbol `□` denotes the infimal convolution, and `L` is a linear mapping.

The returned object has type `IterativeAlgorithm{AFBAIteration}`,
and can be called with the problem's arguments to trigger its solution.

See also: [`AFBAIteration`](@ref), [`IterativeAlgorithm`](@ref).

# Arguments
- `maxit::Int=10_000`: maximum number of iteration
- `tol::1e-5`: tolerance for the default stopping criterion
- `stop::Function`: termination condition, `stop(::T, state)` should return `true` when to stop the iteration
- `solution::Function`: solution mapping, `solution(::T, state)` should return the identified solution
- `verbose::Bool=false`: whether the algorithm state should be displayed
- `freq::Int=100`: every how many iterations to display the algorithm state
- `display::Function`: display function, `display(::Int, ::T, state)` should display a summary of the iteration state
- `kwargs...`: additional keyword arguments to pass on to the `AFBAIteration` constructor upon call

# References
1. Latafat, Patrinos, "Asymmetric forward-backward-adjoint splitting for solving monotone inclusions involving three operators", Computational Optimization and Applications, vol. 68, no. 1, pp. 57-93 (2017).
2. Latafat, Patrinos, "Primal-dual proximal algorithms for structured convex optimization: a unifying framework", In Large-Scale and Distributed Optimization, Giselsson and Rantzer, Eds. Springer International Publishing, pp. 97-120 (2018).
"""
AFBA(;
    maxit=10_000,
    tol=1e-5,
    stop=(iter, state) -> default_stopping_criterion(tol, iter, state),
    solution=default_solution,
    verbose=false,
    freq=100,
    display=default_display,
    kwargs...
) = IterativeAlgorithm(AFBAIteration; maxit, stop, solution, verbose, freq, display, kwargs...)

"""
    VuCondat(; <keyword-arguments>)

Constructs the Vũ-Condat primal-dual algorithm [1, 2].

This algorithm solves convex optimization problems of the form

    minimize f(x) + g(x) + (h □ l)(L x),

where `f` is smooth, `g` and `h` are possibly nonsmooth and `l` is strongly
convex. Symbol `□` denotes the infimal convolution, and `L` is a linear mapping.

The returned object has type `IterativeAlgorithm{AFBAIteration}`,
and can be called with the problem's arguments to trigger its solution.

See also: [`VuCondatIteration`](@ref), [`AFBAIteration`](@ref), [`IterativeAlgorithm`](@ref).

# Arguments
- `maxit::Int=10_000`: maximum number of iteration
- `tol::1e-5`: tolerance for the default stopping criterion
- `stop::Function`: termination condition, `stop(::T, state)` should return `true` when to stop the iteration
- `solution::Function`: solution mapping, `solution(::T, state)` should return the identified solution
- `verbose::Bool=false`: whether the algorithm state should be displayed
- `freq::Int=100`: every how many iterations to display the algorithm state
- `display::Function`: display function, `display(::Int, ::T, state)` should display a summary of the iteration state
- `kwargs...`: additional keyword arguments to pass on to the `AFBAIteration` constructor upon call

# References
1. Condat, "A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms", Journal of Optimization Theory and Applications, vol. 158, no. 2, pp 460-479 (2013).
2. Vũ, "A splitting algorithm for dual monotone inclusions involving cocoercive operators", Advances in Computational Mathematics, vol. 38, no. 3, pp. 667-681 (2013).
"""
VuCondat(; kwargs...) = AFBA(; kwargs..., theta=2)

"""
    ChambollePock(; <keyword-arguments>)

Constructs the Chambolle-Pock primal-dual algorithm [1].

This algorithm solves convex optimization problems of the form

    minimize g(x) + h(L x),

where `g` and `h` are possibly nonsmooth, and `L` is a linear mapping.

The returned object has type `IterativeAlgorithm{AFBAIteration}`,
and can be called with the problem's arguments to trigger its solution.

See also: [`ChambollePockIteration`](@ref), [`AFBAIteration`](@ref), [`IterativeAlgorithm`](@ref).

# Arguments
- `maxit::Int=10_000`: maximum number of iteration
- `tol::1e-5`: tolerance for the default stopping criterion
- `stop::Function`: termination condition, `stop(::T, state)` should return `true` when to stop the iteration
- `solution::Function`: solution mapping, `solution(::T, state)` should return the identified solution
- `verbose::Bool=false`: whether the algorithm state should be displayed
- `freq::Int=100`: every how many iterations to display the algorithm state
- `display::Function`: display function, `display(::Int, ::T, state)` should display a summary of the iteration state
- `kwargs...`: additional keyword arguments to pass on to the `AFBAIteration` constructor upon call

# References
1. Chambolle, Pock, "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging", Journal of Mathematical Imaging and Vision, vol. 40, no. 1, pp. 120-145 (2011).
"""
ChambollePock(; kwargs...) = AFBA(; kwargs..., f=Zero(), l=IndZero(), theta=2)

function AFBA_default_stepsizes(L, h::Zero, theta::R, mu::R, beta_f::R, beta_l::R) where R
    return R(1.99) / beta_f, R(1)
end

function AFBA_default_stepsizes(L, h, theta::R, mu::R, beta_f::R, beta_l::R) where R
    par = R(5) # scaling parameter for comparing Lipschitz constants and \|L\|
    par2 = R(100)   # scaling parameter for α
    alpha = R(1)
    nmL = R(opnorm(L))

    if theta ≈ 2 # default stepsize for Vu-Condat
        if nmL > par * max(beta_l, beta_f)
            alpha = R(1)
        elseif beta_f > par * beta_l
            alpha = par2 * nmL / beta_f
        elseif beta_l > par * beta_f
            alpha = beta_l / (par2 * nmL)
        end
        gamma1 = R(1) / (beta_f / 2 + nmL / alpha)
        gamma2 = R(0.99) / (beta_l / 2 + nmL * alpha)
    elseif theta ≈ 1 && mu ≈ 1 # SPCA
        if nmL > par2 * beta_l # for the case beta_f = 0
            alpha = R(1)
        elseif beta_l > par * beta_f
            alpha = beta_l / (par2 * nmL)
        end
        gamma1 = beta_f > 0 ? R(1.99) / beta_f : R(1) / (nmL / alpha)
        gamma2 = R(0.99) / (beta_l / 2 + gamma1 * nmL^2)
    elseif theta ≈ 0 && mu ≈ 1 # PPCA
        temp = R(3)
        if beta_f ≈ 0
            nmL *= sqrt(temp)
            if nmL > par * beta_l
                alpha = R(1)
            else
                alpha = beta_l / (par2 * nmL)
            end
            gamma1 = R(1) / (beta_f / 2 + nmL / alpha)
            gamma2 = R(0.99) / (beta_l / 2 + nmL * alpha)
        else
            if nmL > par * max(beta_l, beta_f)
                alpha = R(1)
            elseif beta_f > par * beta_l
                alpha = par2 * nmL / beta_f
            elseif beta_l > par * beta_f
                alpha = beta_l / (par2 * nmL)
            end
            xi = 1 + 2 * nmL / (nmL + alpha * beta_f / 2)
            gamma1 = R(1) / (beta_f / 2 + nmL / alpha)
            gamma2 = R(0.99) / (beta_l / 2 + xi * nmL * alpha)
        end
    elseif mu ≈ 0 # SDCA & PDCA
        temp = theta^2 - 3 * theta + 3
        if beta_l ≈ 0
            nmL *= sqrt(temp)
            if nmL > par * beta_f
                alpha = R(1)
            else
                alpha = par2 * nmL / beta_f
            end
            gamma1 = R(1) / (beta_f / 2 + nmL / alpha)
            gamma2 = R(0.99) / (beta_l / 2 + nmL * alpha)
        else
            if nmL > par * max(beta_l, beta_f)
                alpha = R(1)
            elseif beta_f > par * beta_l
                alpha = par2 * nmL / beta_f
            elseif beta_l > par * beta_f
                alpha = beta_l / (par2 * nmL)
            end
            eta = 1 + (temp - 1) * alpha * nmL / (alpha * nmL + beta_l / 2)
            gamma1 = R(1) / (beta_f / 2 + eta * nmL / alpha)
            gamma2 = R(0.99) / (beta_l / 2 + nmL * alpha)
        end
    elseif theta ≈ 0 && mu ≈ 0.5 # PPDCA
        if beta_l ≈ 0 || beta_f ≈ 0
            if nmL > par * max(beta_l, beta_f)
                alpha = R(1)
            elseif beta_f > par * beta_l
                alpha = par2 * nmL / beta_f
            elseif beta_l > par * beta_f
                alpha = beta_l / (par2 * nmL)
            end
        else
            alpha = sqrt(beta_l / beta_f) / 2
        end
        gamma1 = R(1) / (beta_f / 2 + nmL / alpha)
        gamma2 = R(0.99) / (beta_l / 2 + nmL * alpha)
    else
        error("this choice of theta and mu is not supported!")
    end

    return gamma1, gamma2
end


text:
module util

function load_input(day::Int, strip::Bool = true)::String
    s = open("inputs/day_$day.txt") do f
        read(f, String)
    end
    if strip
        Base.strip(s)
    else
        s
    end
end

end


text:
" Surface area scaling rules "
abstract type AbstractShape end

struct Isomorph <: AbstractShape end

@columns struct V0morph{Mo} <: AbstractShape
    Vd::Mo  | 4.0 | mol | Gamma(2.0, 0.2) | [0.0, 1000.0] | _ | "reference"
end

@columns struct V1morph{Mo} <: AbstractShape
    Vd::Mo  | 4.0 | mol | Gamma(2.0, 0.2) | [0.0, 1000.0] | _ | "reference"
end

@columns struct V1V0morph{Mo,B} <: AbstractShape
    Vd::Mo   | 4.0 | mol | Gamma(2.0, 0.2) | [0.0, 1000.0] | _ | "reference"
    Vmax::Mo | 4.0 | mol | Gamma(2.0, 0.2) | [0.0, 1000.0] | _ | "reference"
    β::B     | 4.0 | mol | Gamma(2.0, 0.2) | [0.0, 10.0] | _ | "reference"
end

@columns struct Plantmorph{Mo} <: AbstractShape
    M_Vref::Mo     | 0.1 | mol | Gamma(2.0, 0.2) | [0.00001, 20.0]  | true | "Scaling reference"
    M_Vscaling::Mo | 1.0 | mol | Gamma(2.0, 0.2) | [0.0001, 2000.0] | true | "Scaling mass"
end

update_shape!(o, u) = set_shape!(o, shape_correction(shape_pars(o), u.V))

shape_correction(f::Nothing, V) = 1
shape_correction(f::Isomorph, V) = 1
shape_correction(f::V0morph, V) = (V / f.Vd)^(-2//3)
shape_correction(f::V1morph, V) = (V / f.Vd)^(1//3)
shape_correction(f::V1V0morph, V) = (V / f.Vd)^(1//3 - (V/f.Vmax)^f.β)
shape_correction(f::Plantmorph, V) = (V / f.M_Vref)^(-V / f.M_Vscaling)


text:
export TransformerPowerAnalysis, transformer_power_analysis
export total_power, winding_power, r_core, input_power, output_power
export efficiency

"""
    TransformerPowerAnalysis

Object to store the result of a transformer_power_analysis().

**Fields**
- `transformer`             -- from input
- `converter`               -- type of converter.
- `flux_density`            -- peak flux density in Tesla
- `voltage`                 -- (primary, secondary) Vpp
- `current`                 -- (promary, secondary) Ipp
- `ac_winding_resistance`   -- (primary, secondary) Ω
- `dc_winding_resistance`   -- (primary, secondary) Ω
- `equilivent_resistance`   -- (primary, secondary) Ω
- `v_core`                  -- voltage, after drop accross input resistance, referenced to 1 turn
- `core_specific_power`     -- power dissipated in core in W/m^3
- `core_total_power`        -- power dissipated in core in W
"""
struct TransformerPowerAnalysis
  transformer :: Transformer
  converter :: Converter
  flux_density ::Float64
  voltage :: Tuple # (primary, secondary)
  current :: Tuple # (primary, secondary)
  ac_winding_resistance :: Tuple
  dc_winding_resistance :: Tuple
  equilivent_resistance :: Tuple
  v_core  :: Float64 # volatge referenced to 1T corrected for input resistance loss
  core_specific_power :: Float64
  core_total_power :: Float64
end

"""
    TransformerPowerAnalysis(t::Transformer,
                             c::Converter)

Determine the losses of a transformer.
"""
function transformer_power_analysis(t::Transformer, c::Converter)
  n = turns(t)[2]/turns(t)[1]
  seconds = duty(c)/frequency(c)
  r = equilivent_resistance(t,c) #winding_resistance(t, frequency(c))
  # initial conditions
  v1 = v_in(c)
  i2 = i_out(c)
  i1 = -i2 * n
  core_total_power = 0.0
  core_specific_power = 0.0
  v3 = 0.0
  v2 = 0.0
  flux_density = 0.0
  i1_old = 0.0
  for j in 1:10
    # estimate v2
    v3 = v1 - i1*r[1] # v3 referenced to primary
    v3 = v3>0.0 ? v3 : 0.0 # avoid log of negative number. unlikely? delete line?
    v2 = n*v3 + i2*r[2]
    # estimate average core power dissipation
    flux_density = v3*seconds/(2.0*turns(t)[1]*effective_area(t)) # fix it one more time
    core_specific_power = specific_power_loss(t, flux_density, frequency(c))
    core_total_power = effective_volume(t) * core_specific_power
    # estimate i1
    i1_old = i1
    i1 = new_i1(c, core_total_power, i2, v3, n)
    # println(abs(i1_old-i1),"    ",i1)
    if abs(i1_old-i1)<eps()*10
      break
    end
   end
  v3 /= turns(t)[1] # reference v3 to 1 turn
  TransformerPowerAnalysis(t, c, flux_density, (v1, v2), (i1, i2),
    winding_resistance(t, frequency(c)), winding_resistance(t,0.0),
    r,
    v3, core_specific_power,
      core_total_power)
end
function equilivent_resistance(t::Transformer, c::PushPull)
  winding_resistance(t, frequency(c))
end
function equilivent_resistance(t::Transformer, c::Forward)
  (1.0-duty(c)).*winding_resistance(t, frequency(c)) .+ duty(c).*winding_resistance(t, 0.0)
end
new_i1(c::PushPull, core_total_power::Float64, i2::Float64, v3::Float64, n::Float64) = core_total_power/v3 - i2*n
new_i1(c::Forward, core_total_power::Float64, i2::Float64, v3::Float64, n::Float64) = core_total_power/(v3*duty(c)) - i2*n

"""
    total_power(tpa::TransformerPowerAnalysis)

Total power dissipation of transformer.  Sum of winding and core losses.
"""
total_power(tpa::TransformerPowerAnalysis) = sum(winding_power(tpa)) + core_total_power(tpa)

"""
    winding_power(tpa::TransformerPowerAnalysis)

Power dissipation due to windign resistance.  Tuple (primary, secondary).
"""
winding_power(tpa::TransformerPowerAnalysis) = winding_power(tpa, converter(tpa))
function winding_power(tpa::TransformerPowerAnalysis, c::PushPull)
  current(tpa).^2 .* equilivent_resistance(tpa)
end
function winding_power(tpa::TransformerPowerAnalysis, c::Forward)
  duty(c).*current(tpa).^2 .* equilivent_resistance(tpa)
end

"""
    r_core(tpa::TransformerPowerAnalysis)

Resistance, referenced to 1 turn, for spice model to match core power dissipation.
"""
r_core(tpa::TransformerPowerAnalysis) = r_core(tpa, converter(tpa))
function r_core(tpa::TransformerPowerAnalysis, c::PushPull)
  v_core(tpa)^2/core_total_power(tpa)
end
function r_core(tpa::TransformerPowerAnalysis, c::Forward)
  duty(c)*v_core(tpa)^2/core_total_power(tpa)
end
"""
    power_error(ta::TransformerPowerAnalysis)

returns the difference between total_power and the power loss computed from the
voltages and currents.
"""
power_error(ta::TransformerPowerAnalysis) = total_power(ta) - (input_power(ta) + output_power(ta))

"""
    input_power(tpa::TransformerPowerAnalysis)

Average power delivered to the primary of the transformer.
"""
input_power(tpa::TransformerPowerAnalysis) = input_power(tpa, converter(tpa))
input_power(tpa::TransformerPowerAnalysis, c::PushPull) = voltage(tpa)[1]*current(tpa)[1]
input_power(tpa::TransformerPowerAnalysis, c::Forward) = voltage(tpa)[1]*current(tpa)[1]*duty(c)

"""
    output_power(tpa::TransformerPowerAnalysis)

Average power from the secondary of the transformer.
"""
output_power(tpa::TransformerPowerAnalysis) = output_power(tpa, converter(tpa))
output_power(tpa::TransformerPowerAnalysis, c::PushPull) = voltage(tpa)[2]*current(tpa)[2]
output_power(tpa::TransformerPowerAnalysis, c::Forward) = voltage(tpa)[2]*current(tpa)[2]*duty(c)

"""
    efficiency(tpa::TransformerPowerAnalysis)

Efficiency (-Pout/Pin) of transformer.
"""
efficiency(tpa::TransformerPowerAnalysis) = -output_power(tpa)/input_power(tpa)

core_name(tpa::TransformerPowerAnalysis) = core_name(transformer(tpa))
ferrite_name(tpa::TransformerPowerAnalysis) = ferrite_name(transformer(tpa))


text:
using UAPIC
using Test

@testset " Test particles input file read " begin

mesh = Mesh( 0, 4π, 128, 0, 2π, 64 )

particles = read_particles("particles.dat", mesh)

@test particles.nbpart == 204800

end

@testset " Test Particles-MeshFields interaction " begin

     xmin, xmax = 0.0, 20.0
     ymin, ymax = 0.0, 20.0
     nx,   ny   = 20, 20

     mesh = Mesh( xmin, xmax, nx, ymin, ymax, ny )

     dx, dy = mesh.dx, mesh.dy

     fields = MeshFields( mesh )

     nbpart = 121
     w      = 1/nbpart

     particles = Particles( nbpart, w )

     k = 1
     for i = 5:nx-5, j = 5:ny-5

         particles.x[1,k] = (i-0.5)*dx
         particles.x[2,k] = (j-0.5)*dx

         k += 1

     end

     compute_rho_m6!( fields, particles )

     @test integrate( fields.ρ, mesh ) ≈ 0.0 atol = 1e-4

     for i = 1:nx+1, j = 1:ny+1

         fields.e[1,i,j] = (i-1)*dx
         fields.e[2,i,j] = (j-1)*dy

     end

     interpol_eb_m6!( particles, fields )

     err_x, err_y = 0.0, 0.0

     for k = 1:nbpart

         xp = particles.x[1,k]
         yp = particles.x[2,k]
         ex = particles.e[1,k]
         ey = particles.e[2,k]

         err_x += abs( particles.e[1,k] - xp )
         err_y += abs( particles.e[2,k] - yp )

     end

     err_x /= nbpart
     err_y /= nbpart

     @test err_x ≈ 0.0 atol = 1e-6
     @test err_y ≈ 0.0 atol = 1e-6

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 058a3333-0567-43b7-ac5f-1f6688325a08
begin
	"""
	Mark an instance of a custom struct as immutable. The resulting object is also an `AbstractDict`, where the keys are the struct fields (converted to strings).
	"""
	struct ImmutableMarker{T} <: AbstractDict{String,Any}
		source::T
	end

	function Base.getindex(ldict::ImmutableMarker, key::String)
	    Base.getfield(ldict.source, Symbol(key))
	end
	# disabled because it's immutable!
	# Base.setindex!(ldict::ImmutableMarker, args...) = Base.setindex!(ldict.source, args...)
	# Base.delete!(ldict::ImmutableMarker, args...) = Base.delete!(ldict.source, args...)
	Base.keys(ldict::ImmutableMarker{T}) where T = String.(fieldnames(T))
	# Base.values(ldict::ImmutableMarker) = Base.values(ldict.source)
	Base.length(ldict::ImmutableMarker) = nfields(ldict.source)
	Base.iterate(ldict::ImmutableMarker) = Base.iterate(ldict, 1)
	function Base.iterate(ldict::ImmutableMarker{T}, i) where T
		a = ldict.source

		if i <= nfields(a)
			name = fieldname(T, i)
			(String(name) => getfield(a, name), i + 1)
		end
	end

end

# ╔═╡ 273c7c85-8178-44a7-99f0-581754aeb8c8
begin
	"""
	Mark a vector as being append-only: let Firebasey know that it can diff this array simply by comparing lengths, without looking at its contents.

	It was made specifically for logs: Logs are always appended, OR the whole log stream is reset. AppendonlyMarker is like SubArray (a view into another array) except we agree to only ever append to the source array. This way, firebase can just look at the index and diff based on that.
	"""
	struct AppendonlyMarker{T} <: AbstractVector{T}
		mutable_source::Vector{T}
		length_at_time_of_creation::Int
	end
	AppendonlyMarker(arr) = AppendonlyMarker(arr, length(arr))

	# Poor mans vector-proxy
	# I think this is enough for Pluto to show, and for msgpack to pack
	function Base.size(arr::AppendonlyMarker)
		return (arr.length_at_time_of_creation,)
	end
	Base.getindex(arr::AppendonlyMarker, index::Int) = arr.mutable_source[index]
	Base.iterate(arr::AppendonlyMarker, args...) = Base.iterate(arr.mutable_source, args...)
end

# ╔═╡ ef7032d1-a666-48a6-a56e-df175f5ed832
md"""
## ImmutableMarker
"""

# ╔═╡ 183cef1f-bfe9-42cd-8239-49e9ed00a7b6
md"""
## AppendonlyMarker(s)

Example of how to solve performance problems with Firebasey:
We make a new type with a specific diff function.
It might be very specific per problem, but that's fine for performance problems (I think).
It also keeps the performance solutions as separate modules/packages to whatever it is you're actually modeling.
"""

# ╔═╡ 2284ae12-5b8c-4542-81fa-c4d34f2483e7
# @skip_as_script @test length([AppendonlyMarker([1,2,3], 1)...]) == 1

# ╔═╡ 971709de-074e-49cf-8bd4-9c675b037dfd
md"## `@skip_as_script`"

# ╔═╡ 9fce9aa9-d3c6-4134-9692-8a8756fa3cff
function is_inside_pluto(m::Module)
	if isdefined(m, :PlutoForceDisplay)
		return m.PlutoForceDisplay
	else
		isdefined(m, :PlutoRunner) && parentmodule(m) == Main
	end
end

# ╔═╡ db9167c4-7ba7-42e1-949b-0ad18e2d7b25
"""
	@skip_as_script expression

Marks a expression as Pluto-only, which means that it won't be executed when running outside Pluto. Do not use this for your own projects.
"""
macro skip_as_script(ex)
	if is_inside_pluto(__module__)
		esc(ex)
	else
		nothing
	end
end

# ╔═╡ 55975e53-f70f-4b70-96d2-b144f74e7cde
@skip_as_script struct A
	x
	y
	z
end

# ╔═╡ d7e0de85-5cb2-4036-a2e3-ca416ea83737
@skip_as_script id1 = ImmutableMarker(A(1,"asdf",3))

# ╔═╡ 08350326-526e-4c34-ab27-df9fbf69243e
@skip_as_script id2 = ImmutableMarker(A(1,"asdf",4))

# ╔═╡ aa6192e8-410f-4924-8250-4775e21b1590
@skip_as_script id1d, id2d = Dict(id1), Dict(id2)

# ╔═╡ 35d3bcd7-af51-466a-b4c4-cc055e74d01d
@skip_as_script appendonly_1, appendonly_2 = let
	array_1 = [1,2,3,4]
	appendonly_1 = AppendonlyMarker(array_1)
	push!(array_1, 5)
	appendonly_2 = AppendonlyMarker(array_1)

	appendonly_1, appendonly_2
end;

# ╔═╡ 1017f6cc-58ac-4c7b-a6d0-a03f5e387f1b
@skip_as_script appendonly_1_large, appendonly_2_large = let
	large_array_1 = [
		Dict{String,Any}(
			"x" => 1,
			"y" => [1,2,3,4],
			"z" => "hi",
		)
		for i in 1:10000
	];
	appendonly_1 = AppendonlyMarker(large_array_1)
	push!(large_array_1, Dict("x" => 5))
	appendonly_2 = AppendonlyMarker(large_array_1)

	appendonly_1, appendonly_2
end;

# ╔═╡ 33facb62-af28-4e94-946d-545637f320e9
"The opposite of `@skip_as_script`"
macro only_as_script(ex) is_inside_pluto(__module__) ? nothing : esc(ex) end

# ╔═╡ 092c4b11-8b75-446f-b3ad-01fa858daebb
# Only define this in Pluto - assume we are `using Test` otherwise
begin
	@skip_as_script begin
		import Pkg
		Pkg.activate(mktempdir())
		Pkg.add(Pkg.PackageSpec(name="PlutoTest"))
		using PlutoTest
	end
	# Do nothing inside pluto (so we don't need to have Test as dependency)
	# test/Firebasey is `using Test` before including this file
	@only_as_script begin
		if !isdefined(@__MODULE__, Symbol("@test"))
			macro test(e...) nothing; end
			macro test_throws(e...) nothing; end
			macro test_broken(e...) nothing; end
			macro testset(e...) nothing; end
		end
	end
end

# ╔═╡ 06492e8d-4500-4efe-80ee-55bf1ee2348c
@skip_as_script @test length([AppendonlyMarker([1,2,3])...]) == 3

# ╔═╡ dc5cd268-9cfb-49bf-87fb-5b7db4fa6e3c
md"## Import Firebasey when running inside notebook"

# ╔═╡ 0c2f23d8-8e98-47b7-9c4f-5daa70a6c7fb
# OH how I wish I would put in the time to refactor with fromFile or SOEMTGHINLAS LDKJ JULIA WHY ARE YOU LIKE THIS GROW UP
@skip_as_script Firebasey = let
	wrapper_module = Module()
	Core.eval(wrapper_module, :(module Firebasey
			include("Firebasey.jl")
		end
	))
	wrapper_module.Firebasey
end

# ╔═╡ 2903d17e-c6fd-4cea-8585-4db26a00b0e7
function Firebasey.diff(a::AppendonlyMarker, b::AppendonlyMarker)
	if a.mutable_source !== b.mutable_source
		[Firebasey.ReplacePatch([], b)]
	else
		if a.length_at_time_of_creation > b.length_at_time_of_creation
			throw(ErrorException("Not really supposed to diff AppendonlyMarker with the original being longer than the next version (you know, 'append only' and al)"))
		end

		map(a.length_at_time_of_creation+1:b.length_at_time_of_creation) do index
			Firebasey.AddPatch([index], b.mutable_source[index])
		end
	end
end

# ╔═╡ 129dee79-61c0-4524-9bef-388837f035bb
function Firebasey.diff(a::ImmutableMarker, b::ImmutableMarker)
	if a.source !== b.source
		Firebasey.diff(Dict(a), Dict(b))
		# Firebasey.JSONPatch[Firebasey.ReplacePatch([], b)]
	else
		Firebasey.JSONPatch[]
	end
end

# ╔═╡ 138d2cc2-59ba-4f76-bf66-ecdb98cf4fd5
@skip_as_script Firebasey.diff(id1, id2)

# ╔═╡ 8537488d-2ff9-42b7-8bfc-72d43fca713f
@skip_as_script @test Firebasey.diff(appendonly_1, appendonly_2) == [Firebasey.AddPatch([5], 5)]

# ╔═╡ 70179239-357a-424d-bac3-3a1431aff536
var"@track" = Firebasey.var"@track"

# ╔═╡ a5f43f47-6189-413f-95a0-d98f927bb7ce
@skip_as_script @track for _ in 1:1000 Firebasey.diff(id1, id1) end

# ╔═╡ ab5089cc-fec8-43b9-9aa4-d6fa96e231e0
@skip_as_script @track for _ in 1:1000 Firebasey.diff(id1d, id1d) end

# ╔═╡ a84dcdc3-e9ed-4bf5-9bec-c9cbfc267c17
@skip_as_script @track for _ in 1:1000 Firebasey.diff(id1, id2) end

# ╔═╡ f696bb85-0bbd-43c9-99ea-533816bc8e0d
@skip_as_script @track for _ in 1:1000 Firebasey.diff(id1d, id2d) end

# ╔═╡ 37fe8c10-09f0-4f72-8cfd-9ce044c78c13
@skip_as_script @track for _ in 1:1000 Firebasey.diff(appendonly_1_large, appendonly_2_large) end

# ╔═╡ Cell order:
# ╟─ef7032d1-a666-48a6-a56e-df175f5ed832
# ╠═058a3333-0567-43b7-ac5f-1f6688325a08
# ╠═129dee79-61c0-4524-9bef-388837f035bb
# ╠═55975e53-f70f-4b70-96d2-b144f74e7cde
# ╠═d7e0de85-5cb2-4036-a2e3-ca416ea83737
# ╠═08350326-526e-4c34-ab27-df9fbf69243e
# ╠═138d2cc2-59ba-4f76-bf66-ecdb98cf4fd5
# ╠═aa6192e8-410f-4924-8250-4775e21b1590
# ╠═a5f43f47-6189-413f-95a0-d98f927bb7ce
# ╠═ab5089cc-fec8-43b9-9aa4-d6fa96e231e0
# ╠═a84dcdc3-e9ed-4bf5-9bec-c9cbfc267c17
# ╠═f696bb85-0bbd-43c9-99ea-533816bc8e0d
# ╟─183cef1f-bfe9-42cd-8239-49e9ed00a7b6
# ╠═273c7c85-8178-44a7-99f0-581754aeb8c8
# ╠═2903d17e-c6fd-4cea-8585-4db26a00b0e7
# ╠═35d3bcd7-af51-466a-b4c4-cc055e74d01d
# ╠═1017f6cc-58ac-4c7b-a6d0-a03f5e387f1b
# ╠═06492e8d-4500-4efe-80ee-55bf1ee2348c
# ╠═2284ae12-5b8c-4542-81fa-c4d34f2483e7
# ╠═8537488d-2ff9-42b7-8bfc-72d43fca713f
# ╠═37fe8c10-09f0-4f72-8cfd-9ce044c78c13
# ╟─971709de-074e-49cf-8bd4-9c675b037dfd
# ╟─9fce9aa9-d3c6-4134-9692-8a8756fa3cff
# ╟─db9167c4-7ba7-42e1-949b-0ad18e2d7b25
# ╟─33facb62-af28-4e94-946d-545637f320e9
# ╟─dc5cd268-9cfb-49bf-87fb-5b7db4fa6e3c
# ╠═0c2f23d8-8e98-47b7-9c4f-5daa70a6c7fb
# ╠═70179239-357a-424d-bac3-3a1431aff536
# ╠═092c4b11-8b75-446f-b3ad-01fa858daebb


text:
using OrdinaryDiffEq, ModelingToolkit, Test, SparseArrays

N = 3
xyd_brusselator = range(0,stop=1,length=N)
brusselator_f(x, y, t) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.
limit(a, N) = ModelingToolkit.ifelse(a == N+1, 1, ModelingToolkit.ifelse(a == 0, N, a))
function brusselator_2d_loop(du, u, p, t)
  A, B, alpha, dx = p
  alpha = alpha/dx^2
  @inbounds for I in CartesianIndices((N, N))
    i, j = Tuple(I)
    x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
    ip1, im1, jp1, jm1 = limit(i+1, N), limit(i-1, N), limit(j+1, N), limit(j-1, N)
    du[i,j,1] = alpha*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) +
                B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t)
    du[i,j,2] = alpha*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) +
                A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
    end
end

# Test with tuple parameters
p = (3.4, 1., 10., step(xyd_brusselator))

function init_brusselator_2d(xyd)
  N = length(xyd)
  u = zeros(N, N, 2)
  for I in CartesianIndices((N, N))
    x = xyd[I[1]]
    y = xyd[I[2]]
    u[I,1] = 22*(y*(1-y))^(3/2)
    u[I,2] = 27*(x*(1-x))^(3/2)
  end
  u
end

u0 = init_brusselator_2d(xyd_brusselator)
prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
                                        u0,(0.,11.5),p)
sys = modelingtoolkitize(prob_ode_brusselator_2d)

# test sparse jacobian pattern only.
prob = ODEProblem(sys, u0, (0, 11.5), sparse=true, jac=false)
JP = prob.f.jac_prototype
@test findnz(Symbolics.jacobian_sparsity(map(x->x.rhs, equations(sys)), states(sys)))[1:2] == findnz(JP)[1:2]

# test sparse jacobian
prob = ODEProblem(sys, u0, (0, 11.5), sparse=true, jac=true)
@test findnz(calculate_jacobian(sys, sparse=true))[1:2] == findnz(prob.f.jac_prototype)[1:2]

# test when not sparse
prob = ODEProblem(sys, u0, (0, 11.5), sparse=false, jac=true)
@test prob.f.jac_prototype == nothing

prob = ODEProblem(sys, u0, (0, 11.5), sparse=false, jac=false)
@test prob.f.jac_prototype == nothing

# test when u0 is nothing
f = DiffEqBase.ODEFunction(sys, u0=nothing, sparse=true, jac=true)
@test findnz(f.jac_prototype)[1:2] == findnz(JP)[1:2]
@test eltype(f.jac_prototype) == Float64

f = DiffEqBase.ODEFunction(sys, u0=nothing, sparse=true, jac=false)
@test findnz(f.jac_prototype)[1:2] == findnz(JP)[1:2]
@test eltype(f.jac_prototype) == Float64

# test when u0 is not Float64
u0 = similar(init_brusselator_2d(xyd_brusselator), Float32)
prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
                                        u0,(0.,11.5),p)
sys = modelingtoolkitize(prob_ode_brusselator_2d)

prob = ODEProblem(sys, u0, (0, 11.5), sparse=true, jac=false)
@test eltype(prob.f.jac_prototype) == Float32

prob = ODEProblem(sys, u0, (0, 11.5), sparse=true, jac=true)
@test eltype(prob.f.jac_prototype) == Float32

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
using Luxor, Random

Random.seed!(2021)

function meshes(s)
  Δ = s/600 # scale factor, original design was 600units
  scale(Δ)
  # background
  squircle(O, 285, 285, :clip, rt=0.2)
  sethue("slateblue4")
  paint()

  # build point array
  pts = polysample(box(O, 600, 600, vertices=true), 16)
  for pt in ngon(O, 100, 3, π/6, vertices=true)
    for θ in 0:π/8:6π
      for i in 160:40:265
        push!(pts, pt  + polar(i, θ))
      end
    end
  end

  # make triangulation
  tris = polytriangulate(pts)

  # draw triangles
  setline(.5)
  for tri in tris
    sethue("grey")
    poly(tri, :stroke, close=true)
  end

  # clipping for julia circles
  @layer begin
    for pt in ngon(O + (0, 35), 150, 3, π/6, vertices=true)
      circlepath(pt, 110, :path)
    end
    sethue("white")
    setline(3)
    strokepreserve()
    clip()
    setline(1)
    for tri in tris
      sethue([Luxor.julia_red, Luxor.julia_green, Luxor.julia_purple][rand(1:end)])
      poly(tri, :fillpreserve, close=true)
      sethue("white")
      strokepath()
    end
    sethue("white")
    circle.(pts, 3, :fill)
    clipreset()
  end

  # draw nodes
  sethue("black")
  circle.(pts, 3, :fill)

  #outline
  setline(4)
  sethue("black")
  squircle(O, 285, 285, :stroke, rt=0.2)
end

function logo(s, fname)
  Drawing(s, s, fname)
  origin()
  meshes(s)
  finish()
  preview()
end

function logotext(w, h, fname)
  Drawing(w, h, fname)
  origin()
  table = Table([h], [h, w - h])
  @layer begin
    translate(table[1])
    meshes(h)
  end
  @layer begin
    translate(table[2])
    background("white")
    sethue("black")
    # find all fonts available on Linux with `fc-list | -f 2 -d ":"`
    fontface("Julius Sans One")
    fontsize(h/2.5)
    text("Meshes.jl", halign=:center, valign=:middle)
  end
  @layer begin
    translate(table[1])
    meshes(h)
  end
  finish()
  preview()
end

for ext in [".svg", ".png"]
  logo(240, "../assets/logo"*ext)
  logotext(600, 200, "../assets/logo-text"*ext)
end