text:
export group_algebra

################################################################################
#
#  Basic field access
#
################################################################################

base_ring(A::AlgGrp{T}) where {T} = A.base_ring::parent_type(T)

Generic.dim(A::AlgGrp) = size(multiplication_table(A, copy = false), 1)

elem_type(::Type{AlgGrp{T, S, R}}) where {T, S, R} = AlgGrpElem{T, AlgGrp{T, S, R}}

order_type(::AlgGrp{fmpq, S, R}) where { S, R } = AlgAssAbsOrd{AlgGrp{fmpq, S, R}, elem_type(AlgGrp{fmpq, S, R})}
order_type(::Type{AlgGrp{fmpq, S, R}}) where { S, R } = AlgAssAbsOrd{AlgGrp{fmpq, S, R}, elem_type(AlgGrp{fmpq, S, R})}

order_type(::AlgGrp{T, S, R}) where { T <: NumFieldElem, S, R } = AlgAssRelOrd{T, frac_ideal_type(order_type(parent_type(T)))}
order_type(::Type{AlgGrp{T, S, R}}) where { T <: NumFieldElem, S, R } = AlgAssRelOrd{T, frac_ideal_type(order_type(parent_type(T)))}

@doc Markdown.doc"""
    group(A::AlgGrp) -> Group

> Returns the group defining $A$.
"""
group(A::AlgGrp) = A.group

has_one(A::AlgGrp) = true

function (A::AlgGrp{T, S, R})(c::Array{T, 1}) where {T, S, R}
  length(c) != dim(A) && error("Dimensions don't match.")
  return AlgGrpElem{T, typeof(A)}(A, c)
end

@doc Markdown.doc"""
    multiplication_table(A::AlgGrp; copy::Bool = true) -> Array{RingElem, 2}

> Given an group algebra $A$ this function returns the multiplication table of
> $A$: If the function returns $M$ and the basis of $A$ is $g_1,\dots, g_n$ then
> it holds $g_i \cdot g_j = g_{M[i, j]}$.
"""
function multiplication_table(A::AlgGrp; copy::Bool = true)
  if copy
    return deepcopy(A.mult_table)
  else
    return A.mult_table
  end
end

################################################################################
#
#  Construction
#
################################################################################

@doc Markdown.doc"""
    group_algebra(K::Ring, G; op = *) -> AlgGrp

> Returns the group ring $K[G]$.
> $G$ may be any set and `op` a group operation on $G$.
"""
group_algebra(K::Ring, G; op = *) = AlgGrp(K, G, op = op)

group_algebra(K::Ring, G::GrpAbFinGen) = AlgGrp(K, G)

function group_algebra(K::Field, G; op = *)
  A = AlgGrp(K, G, op = op)
  if iszero(characteristic(K))
    A.issemisimple = 1
  else
    A.issemisimple = isone(gcd(characteristic(K), order(G))) ? 1 : 2
  end
  return A
end

function group_algebra(K::Field, G::GrpAbFinGen)
  A = group_algebra(K, G, op = +)
  A.iscommutative = true
  return A
end

################################################################################
#
#  Commutativity
#
################################################################################

iscommutative_known(A::AlgGrp) = (A.iscommutative != 0)

@doc Markdown.doc"""
    iscommutative(A::AlgGrp) -> Bool

> Returns `true` if $A$ is a commutative ring and `false` otherwise.
"""
function iscommutative(A::AlgGrp)
  if iscommutative_known(A)
    return A.iscommutative == 1
  end
  for i in 1:dim(A)
    for j in 1:dim(A)
      if multiplication_table(A, copy = false)[i, j] != multiplication_table(A, copy = false)[j, i]
        A.iscommutative = 2
        return false
      end
    end
  end
  A.iscommutative = 1
  return true
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, A::AlgGrp)
  compact = get(io, :compact, false)
  if compact
    print(io, "Group algebra of dimension ", dim(A), " over ", base_ring(A))
  else
    print(io, "Group algebra of group\n", group(A), "\nover\n", base_ring(A))
  end
end

################################################################################
#
#  Deepcopy
#
################################################################################

# TODO: This is broken. I have to copy everything carefully by hand.
#function Base.deepcopy_internal(A::AlgGrp, dict::IdDict)
#  B = AlgGrp(base_ring(A), group(A))
#  for x in fieldnames(typeof(A))
#    if x != :base_ring && x != :group && isdefined(A, x)
#      setfield!(B, x, Base.deepcopy_internal(getfield(A, x), dict))
#    end
#  end
#  B.base_ring = A.base_ring
#  B.group = A.group
#  return B
#end

################################################################################
#
#  Equality
#
################################################################################

@doc Markdown.doc"""
    ==(A::AlgGrp, B::AlgGrp) -> Bool

> Returns `true` if $A$ and $B$ are equal and `false` otherwise.
"""
function ==(A::AlgGrp{T}, B::AlgGrp{T}) where {T}
  return base_ring(A) == base_ring(B) && group(A) == group(B)
end

###############################################################################
#
#  Trace Matrix
#
###############################################################################

function _assure_trace_basis(A::AlgGrp{T}) where {T}
  if !isdefined(A, :trace_basis_elem)
    A.trace_basis_elem = Vector{T}(undef, dim(A))
    A.trace_basis_elem[1] = base_ring(A)(dim(A))
    for i=2:length(A.trace_basis_elem)
      A.trace_basis_elem[i] = zero(base_ring(A))
    end
  end
  return nothing
end

@doc Markdown.doc"""
    trace_matrix(A::AlgGrp) -> MatElem

> Returns a matrix $M$ over the base ring of $A$ such that
> $M_{i, j} = \mathrm{tr}(b_i \cdot b_j)$, where $b_1, \dots, b_n$ is the
> basis of $A$.
"""
function trace_matrix(A::AlgGrp)
  _assure_trace_basis(A)
  F = base_ring(A)
  n = dim(A)
  M = zero_matrix(F, n, n)
  for i = 1:n
    M[i,i] = tr(A[i]^2)
  end
  for i = 1:n
    for j = i+1:n
      x = tr(A[i]*A[j])
      M[i,j] = x
      M[j,i] = x
    end
  end
  return M
end

################################################################################
#
#  Center
#
################################################################################

@doc Markdown.doc"""
    center(A::AlgGrp) -> AlgAss, AbsAlgAssMor

> Returns the center $C$ of $A$ and the inclusion $C \to A$.
"""
function center(A::AlgGrp{T}) where {T}
  if iscommutative(A)
    B, mB = AlgAss(A)
    return B, mB
  end

  if isdefined(A, :center)
    return A.center::Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}
  end

  B, mB = AlgAss(A)
  C, mC = center(B)
  mD = compose_and_squash(mB, mC)
  A.center = C, mD
  return C, mD
end

################################################################################
#
#  Conversion to AlgAss
#
################################################################################

function AlgAss(A::AlgGrp{T, S, R}) where {T, S, R}
  K = base_ring(A)
  mult = Array{T, 3}(undef, dim(A), dim(A), dim(A))
  B = basis(A)
  d = dim(A)
  for i in 1:d
    for j in 1:d
      l = multiplication_table(A, copy = false)[i, j]
      for k in 1:d
        if k == l
          mult[i, j, k] = one(K)
        else
          mult[i, j, k] = zero(K)
        end
      end
    end
  end
  B = AlgAss(K, mult, one(A).coeffs)
  B.iscommutative = A.iscommutative
  B.issimple = A.issimple
  B.issemisimple = A.issemisimple
  return B, hom(B, A, identity_matrix(K, dim(A)), identity_matrix(K, dim(A)))
end

################################################################################
#
#  Misc
#
################################################################################

Base.enumerate(G::Generic.PermGroup) = enumerate(AllPerms(G.n))

################################################################################
#
#  Generators
#
################################################################################

# Helper function for gens, changes mid in place
function _merge_elts_in_gens!(left::Vector{Tuple{Int, Int}}, mid::Vector{Tuple{Int, Int}}, right::Vector{Tuple{Int, Int}})
  nl = length(left)
  if length(mid) == 0
    mid = deepcopy(left)
  elseif nl != 0
    if left[nl][1] == mid[1][1]
      t = popfirst!(mid)
      prepend!(mid, left)
      tt = (t[1], mid[nl][2] + t[2])
      mid[nl] = tt
    else
      prepend!(mid, left)
    end
  end
  nm = length(mid)
  if nm == 0
    return deepcopy(right)
  end
  if length(right) == 0
    return mid
  end
  if mid[nm][1] == right[1][1]
    t = pop!(mid)
    append!(mid, right)
    tt = (t[1], mid[nm][2] + t[2])
    mid[nm] = tt
  else
    append!(mid, right)
  end
  return mid
end

@doc Markdown.doc"""
    gens(A::AlgGrp, return_full_basis::Type{Val{T}} = Val{false})
      -> Vector{AlgGrpElem}

> Returns a subset of `basis(A)`, which generates $A$ as an algebra over
> `base_ring(A)`.
> If `return_full_basis` is set to `Val{true}`, the function also returns a
> `Vector{AbsAlgAssElem}` containing a full basis consisting of monomials in
> the generators and a `Vector{Vector{Tuple{Int, Int}}}` containing the
> information on how these monomials are built. E. g.: If the function returns
> `g`, `full_basis` and `v`, then we have
> `full_basis[i] = prod( g[j]^k for (j, k) in v[i] )`.
"""
function gens(A::AlgGrp, return_full_basis::Type{Val{T}} = Val{false}) where T
  G = group(A)
  group_gens = gens(G)

  return_full_basis == Val{true} ? nothing : return map(A, group_gens)

  full_group = elem_type(G)[ id(G) ]
  elts_in_gens = Vector{Tuple{Int, Int}}[ Tuple{Int, Int}[] ]
  constructed_elements = falses(BigInt(order(G)))
  constructed_elements[A.group_to_base[id(G)]] = true
  new_elements = Set{Int}()
  for i = 1:length(group_gens)
    g = group_gens[i]
    push!(full_group, g)
    push!(elts_in_gens, Tuple{Int, Int}[ (i, 1) ])
    constructed_elements[A.group_to_base[g]] = true
    push!(new_elements, length(full_group))
  end
  k = 1 + length(group_gens) # == number of constructed elements, i. e. #{ i | constructed_elements[i] == true }

  while k != dim(A) || !isempty(new_elements)
    i = pop!(new_elements)
    g = full_group[i]

    n = length(full_group)
    for r = 1:n
      s = op(g, full_group[r])
      for l = 1:n
        if !iscommutative(A)
          t = op(full_group[l], s)
        else
          t = s
        end
        if constructed_elements[A.group_to_base[t]]
          continue
        end
        constructed_elements[A.group_to_base[t]] = true
        k += 1
        push!(full_group, t)
        coord = _merge_elts_in_gens!(elts_in_gens[l], deepcopy(elts_in_gens[i]), elts_in_gens[r])
        push!(elts_in_gens, coord)
        push!(new_elements, length(full_group))
        if iscommutative(A)
          break
        end
        k == dim(A) ? break : nothing
      end
      k == dim(A) ? break : nothing
    end
  end

  return map(A, group_gens), map(A, full_group), elts_in_gens
end

################################################################################
#
#  Isomorphisms to number fields
#
################################################################################

# Assumes that Gal(K/k) == group(A), where k = base_field(K) and that group(A) is
# abelian.
# Returns a k-linear map from K to A and one from A to K
function _find_isomorphism(K::Union{ AnticNumberField, NfRel{nf_elem} }, A::AlgGrp)
  G = group(A)
  aut = automorphisms(K)

  aut_dict = Dict{elem_type(K), Int}()
  n = length(aut)
  identity = 0
  for i = 1:n
    b = aut[i].prim_img
    aut_dict[b] = i
    if b == gen(K)
      identity = i
    end
  end

  op_array = zeros(Int, n, n)
  for i = 1:n
    for j = i:n
      if i == identity
        k = j
      elseif j == identity
        k = i
      else
        b = aut[i](aut[j].prim_img)
        k = aut_dict[b]
      end
      op_array[i, j] = k
      # It is assumed, that the group is abelian
      op_array[j, i] = k
    end
  end

  function op(i::Int, j::Int)
    return op_array[i, j]
  end

  auttoG, Gtoaut = find_isomorphism([ i for i in 1:n ], op, G)

  alpha = normal_basis(K)
  basis_alpha = Vector{elem_type(K)}(undef, dim(A))
  for (i, g) in enumerate(G)
    f = aut[Gtoaut[g]]
    basis_alpha[A.group_to_base[g]] = f(alpha)
  end

  M = zero_matrix(base_field(K), degree(K), degree(K))
  for i = 1:degree(K)
    a = basis_alpha[i]
    for j = 1:degree(K)
      M[i, j] = coeff(a, j - 1)
    end
  end

  invM = inv(M)

  function KtoA(x::Union{ nf_elem, NfRelElem })
    t = zero_matrix(base_field(K), 1, degree(K))
    for i = 1:degree(K)
      t[1, i] = coeff(x, i - 1)
    end
    y = t*invM
    return A([ y[1, i] for i = 1:degree(K) ])
  end

  function AtoK(x::AlgGrpElem)
    t = matrix(base_field(K), 1, degree(K), coeffs(x))
    y = t*M
    return K(parent(K.pol)([ y[1, i] for i = 1:degree(K) ]))
  end

  return KtoA, AtoK
end


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

@testset "scalar_dist" begin
end

ext: jl


text:
export batchNormNN,getbatchNormNN,initTheta

"""
batchNormNN Neural Network block

 Y_k+1 = layer{k}(theta{k},Y_k)
"""
#TODO: Can probably optimize some functions using the knowledge that we only have a norm and an AFS layer.
mutable struct batchNormNN{T,TQ <: Union{Array{T,2},UniformScaling{Int}}} <: AbstractMeganetElement{T}
    layers   ::Tuple{normLayer{T}, AffineScalingLayer{T}}
    outTimes ::Array{Int,1}
    Q        ::TQ
end

function getbatchNormNN(layers::Tuple{normLayer{T}, AffineScalingLayer{T}},outTimes=eye(Int,length(layers))[:,end],Q=I) where {T <: Number}
	nt   = length(layers)
    nout = nFeatOut(layers[1])

    for k=2:nt
        if nFeatIn(layers[k]) != nout
            error("Dim. of input features of block $k does not match dim. of output features of block $(k-1)");
        end
        nout = nFeatOut(layers[k])
    end
	return batchNormNN(layers,outTimes,Q);
end

# ---------- counting thetas, input and output features -----
function nTheta(this::batchNormNN)
    n::Int = 0;
    for k=1:length(this.layers)
        n += nTheta(this.layers[k]);
    end
    return n
end
nFeatIn(this::batchNormNN)   = nFeatIn(this.layers[1])
nFeatOut(this::batchNormNN) = nFeatOut(this.layers[end])

function nDataOut(this::batchNormNN)
    n=0;
    for k=1:length(this.layers)
        n = n+this.outTimes[k]* nFeatOut(this.layers[k]);
    end
end

function initTheta(this::batchNormNN{T}) where {T <: Number}
    theta = zeros(T,0)
    for k=1:length(this.layers)
        theta = [theta; vec(initTheta(this.layers[k]))]
    end
    return convert(Array{T},theta)
end

# --------- forward problem ----------
function apply(this::batchNormNN{T},theta::Array{T},Y::Array{T,2},tmp::Array,doDerivative=true) where {T<:Number}
    nex = div(length(Y),nFeatIn(this))::Int
    nt = length(this.layers)

    if isempty(tmp) #TODO Will have to make sure size of Y doesnt change
        tmp = Array{Any}(nt+1,2)
    end

    if doDerivative
        if isassigned(tmp,1,1)
            tmp11 = tmp[1,1]
            tmp11 .= Y
        else
            tmp[1,1] = copy(Y)
        end
    end

    Ydata::Array{T,2} = zeros(T,0,nex)
    cnt = 0
    for i=1:nt
        ni = nTheta(this.layers[i])::Int

        Yd::Array{T,2}, Y, tmp[i,2] = apply(this.layers[i],theta[cnt+(1:ni)],Y,doDerivative)
        if this.outTimes[i]==1
            Ydata = [Ydata; this.Q*Yd]
        end
        if doDerivative
            if isassigned(tmp,i+1,1)
                tmp1 = tmp[i+1,1]
                tmp1 .= Y
            else
                tmp[i+1,1] = copy(Y)
            end
        end
        cnt = cnt + ni
    end

    return Ydata,Y,tmp
end

# -------- Jacobian matvecs --------
function JYmv(this::batchNormNN{T},dY::Array{T},theta::Array{T},Y::Array{T},tmp) where {T <: Number}
    nex = div(length(Y),nFeatIn(this))
    nt = length(this.layers)
    cnt = 0
    dYdata = zeros(T,0,nex)
    for i=1:nt
        ni = nTheta(this.layers[i])
        dY = JYmv(this.layers[i],dY,theta[cnt+(1:ni)],tmp[i,1],tmp[i,2])[2]
        if this.outTimes[i]==1
            dYdata = [dYdata; this.Q*dY]
        end
        cnt = cnt+ni
    end
    return dYdata, dY
end

function  Jmv(this::batchNormNN{T},dtheta::Array{T},dY::Array{T},theta::Array{T},Y::Array{T},tmp) where {T <: Number}
    nex = div(length(Y),nFeatIn(this))
    nt = length(this.layers);
    if isempty(dY)
        dY = 0*Y
    end

    dYdata = zeros(T,0,nex)
    cnt = 0
    for i=1:nt
        ni = nTheta(this.layers[i])
        dY = Jmv(this.layers[i],dtheta[cnt+(1:ni)],dY,theta[cnt+(1:ni)],
                tmp[i,1],tmp[i,2])[2]
        if this.outTimes[i]==1
            dYdata = [dYdata; this.Q*dY]
        end
        cnt = cnt+ni
    end
    return dYdata,dY
end

# -------- Jacobian' matvecs --------
function JYTmv(this::batchNormNN{T},Wdata::Array{T},W::Array{T},theta::Array{T},Y::Array{T},tmp) where {T <: Number}

    nex = div(length(Y),nFeatIn(this));
    if !isempty(Wdata)
        Wdata = reshape(Wdata,:,nex);
    end
    if isempty(W)
        W = zero(T)
    elseif length(W)>1
        W     = reshape(W,:,nex)
    end
    nt = length(this.layers)

    cnt = 0; cnt2 = 0;
    for i=nt:-1:1
        ni = nTheta(this.layers[i])
        if this.outTimes[i]==1
            nn = nFeatOut(this.layers[i])
            W = W + this.Q'*Wdata[end-cnt2-nn+1:end-cnt2,:]
            cnt2 = cnt2 + nn
        end
        W  = JYTmv(this.layers[i], W,(T)[],theta[end-cnt-ni+1:end-cnt],
                    tmp[i,1],tmp[i,2])
        cnt = cnt+ni
    end
    return W
end

function JthetaTmv(this::batchNormNN{T},Wdata::Array{T},W::Array{T},theta::Array{T},Y::Array{T},tmp) where {T <: Number}
	return JTmv(this,Wdata,W,theta,Y,tmp)[1];
end

function JTmv(this::batchNormNN,Wdata::Array{T},W::Array{T},theta::Array{T},Y::Array{T},tmp) where {T <: Number}

    nex = div(length(Y),nFeatIn(this))

    if size(Wdata,1)>0
        Wdata = reshape(Wdata,:,nex)
    end
    if length(W)==0
        W = zeros(T,nFeatOut(this),nex)
    elseif length(W)>1
        W     = reshape(W,:,nex)
    end

    dtheta = 0*theta
    nt = length(this.layers)

    cnt = 0; cnt2 = 0
    for i=nt:-1:1
        if this.outTimes[i]==1
            nn = nFeatOut(this.layers[i])
            W += this.Q'*Wdata[end-cnt2-nn+1:end-cnt2,:]
            cnt2 = cnt2 + nn
        end

        ni     = nTheta(this.layers[i])

        dmbi,W = JTmv(this.layers[i],W,zeros(T,0),theta[end-cnt-ni+1:end-cnt],tmp[i,1],tmp[i,2])
        dtheta[end-cnt-ni+1:end-cnt]  = dmbi
        cnt = cnt+ni
    end
    return  vec(dtheta), vec(W)

end


text:
@symbol_func function phi_CS(cur_reactor::AbstractReactor)
  cur_phi = cur_reactor.k

  cur_phi *= cur_reactor.B_0

  cur_phi *= sqrt( h_CS(cur_reactor) )

  cur_phi *= sqrt( G_LO(cur_reactor) )

  cur_phi *= sqrt( K_LP(cur_reactor) )

  cur_phi *= sqrt( cur_reactor.R_0 )

  cur_phi /= sqrt( 1e-7 )

  cur_phi
end


text:
module TestRDA
    using Base.Test
    using DataFrames
    using RData
    using Compat

    # check for Float64 NA
    @test !RData.isna_float64(reinterpret(UInt64, 1.0))
    @test !RData.isna_float64(reinterpret(UInt64, NaN))
    @test !RData.isna_float64(reinterpret(UInt64, Inf))
    @test !RData.isna_float64(reinterpret(UInt64, -Inf))
    @test RData.isna_float64(reinterpret(UInt64, RData.R_NA_FLOAT64))
    # check that alternative NA is also recognized (#10)
    @test RData.isna_float64(reinterpret(UInt64, RData.R_NA_FLOAT64 | ((Base.significand_mask(Float64) + 1) >> 1)))

    testdir = dirname(@__FILE__)

    df = DataFrame(num = [1.1, 2.2])
    @test isequal(sexp2julia(load("$testdir/data/minimal.rda",convert=false)["df"]), df)
    @test isequal(load("$testdir/data/minimal.rda",convert=true)["df"], df)
    @test isequal(load("$testdir/data/minimal_ascii.rda")["df"], df)

    df[:int] = Int32[1, 2]
    df[:logi] = [true, false]
    df[:chr] = ["ab", "c"]
    df[:factor] = pool(df[:chr])
    df[:cplx] = Complex128[1.1+0.5im, 1.0im]
    @test isequal(sexp2julia(load("$testdir/data/types.rda",convert=false)["df"]), df)
    @test isequal(sexp2julia(load("$testdir/data/types_ascii.rda",convert=false)["df"]), df)

    df[2, :] = NA
    append!(df, df[2, :])
    df[3, :num] = NaN
    df[:, :cplx] = @data [NA, @compat(Complex128(1,NaN)), NaN]
    @test isequal(sexp2julia(load("$testdir/data/NAs.rda",convert=false)["df"]), df)
    # ASCII format saves NaN as NA
    df[3, :num] = NA
    df[:, :cplx] = @data [NA, NA, NA]
    @test isequal(sexp2julia(load("$testdir/data/NAs_ascii.rda",convert=false)["df"]), df)

    rda_names = names(sexp2julia(load("$testdir/data/names.rda",convert=false)["df"]))
    expected_names = [:_end, :x!, :x1, :_B_C_, :x, :x_1]
    @test rda_names == expected_names
    rda_names = names(sexp2julia(load("$testdir/data/names_ascii.rda",convert=false)["df"]))
    @test rda_names == [:_end, :x!, :x1, :_B_C_, :x, :x_1]

    rda_envs = load("$testdir/data/envs.rda",convert=false)

    rda_pairlists = load("$testdir/data/pairlists.rda",convert=false)

    rda_closures = load("$testdir/data/closures.rda",convert=false)

    rda_cmpfuns = load("$testdir/data/cmpfun.rda",convert=false)
end


text:
using SchumakerSpline
using Test

gridx = collect(1:10)
gridy = collect(1:10)
grid = gridx * gridy'

schum = Schumaker2d(gridx, gridy, grid; extrapolation = (Linear, Linear))

abs(schum(5,5) - 25) < 100 * eps()
abs(schum(5.5,5.5) - 30.25) < 100 * eps()


text:
include("QuantumOpsTests.jl");

QuantumOpsTests.runtests()


text:
__precompile__()

module DiffEqJump

using DiffEqBase, Compat, RandomNumbers, TreeViews, LinearAlgebra
using DataStructures, PoissonRandom, Random
using FunctionWrappers, Parameters

import DiffEqBase: DiscreteCallback, init, solve, solve!, plot_indices
import Base: size, getindex, setindex!, length, similar, show
import DataStructures: length, update!

import RecursiveArrayTools: recursivecopy!

abstract type AbstractJump end
abstract type AbstractAggregatorAlgorithm end
abstract type AbstractJumpAggregator end

import Base.Threads
@static if VERSION < v"1.3"
  seed_multiplier() = Threads.threadid()
else
  seed_multiplier() = 1
end

include("jumps.jl")
include("massaction_rates.jl")
include("aggregators/aggregators.jl")
include("aggregators/ssajump.jl")
include("aggregators/direct.jl")
include("aggregators/frm.jl")
include("aggregators/sortingdirect.jl")
include("aggregators/nrm.jl")
include("aggregators/bracketing.jl")
include("aggregators/rssa.jl")
include("aggregators/prioritytable.jl")
include("aggregators/directcr.jl")
include("problem.jl")
include("callbacks.jl")
include("solve.jl")
include("extended_jump_array.jl")
include("coupled_array.jl")
include("coupling.jl")
include("SSA_stepper.jl")
include("simple_regular_solve.jl")

export ConstantRateJump, VariableRateJump, RegularJump, MassActionJump,
       JumpSet

export JumpProblem

export SplitCoupledJumpProblem

export Direct, DirectFW, SortingDirect, DirectCR
export BracketData, RSSA
export FRM, FRMFW, NRM

export get_num_majumps, needs_depgraph, needs_vartojumps_map

export init, solve, solve!

export ExtendedJumpArray

end # module


text:
struct StochasticProgram{N, M, S <: NTuple{N, Stage}, SP <: Union{AbstractScenarioProblems, NTuple{M, AbstractScenarioProblems}}}
    stages::S
    scenarioproblems::SP
    generator::Dict{Symbol, Function}
    problemcache::Dict{Symbol, JuMP.Model}
    spsolver::SPSolver

    function StochasticProgram(first_stage_params::Any, second_stage_params::Any, ::Type{S}, solver::SPSolverType, procs::Vector{Int}) where S <: AbstractScenario
        stages = (Stage(1, first_stage_params), Stage(2, second_stage_params))
        scenarioproblems = ScenarioProblems(S, procs)
        SP = typeof(scenarioproblems)
        return new{2, 1, typeof(stages), SP}(stages, scenarioproblems, Dict{Symbol, Function}(), Dict{Symbol, JuMP.Model}(), SPSolver(solver))
    end

    function StochasticProgram(first_stage_params::Any, second_stage_params::Any, scenarios::Vector{<:AbstractScenario}, solver::SPSolverType, procs::Vector{Int})
        stages = (Stage(1, first_stage_params), Stage(2, second_stage_params))
        S = typeof(stages)
        scenarioproblems = ScenarioProblems(scenarios, procs)
        SP = typeof(scenarioproblems)
        return new{2, 1, S, SP}(stages, scenarioproblems, Dict{Symbol, Function}(), Dict{Symbol, JuMP.Model}(), SPSolver(solver))
    end

    function StochasticProgram(stage_params::NTuple{N, Any}, scenario_types::NTuple{M, DataType}, solver::SPSolverType, procs::Vector{Int}) where {N, M}
        N >= 2 || error("Stochastic program needs at least two stages.")
        M == N - 1 || error("Inconsistent number of stages $N and number of scenario types $M")
        stages = ntuple(Val(N)) do i
            Stage(i, stage_params[i])
        end
        S = typeof(stages)
        scenarioproblems = ntuple(Val(M)) do i
            ScenarioProblems(scenarios[i], procs)
        end
        SP = typeof(scenarioproblems)
        return new{N, M, S, SP}(stages, scenarioproblems, Dict{Symbol, Function}(), Dict{Symbol, JuMP.Model}(), SPSolver(solver))
    end

    function StochasticProgram(stage_params::NTuple{N, Any}, scenarios::NTuple{M, Vector{<:AbstractScenario}}, solver::SPSolverType, procs::Vector{Int}) where {N, M}
        N >= 2 || error("Stochastic program needs at least two stages.")
        M == N - 1 || error("Inconsistent number of stages $N and number of scenario types $M")
        stages = ntuple(Val(N)) do i
            Stage(i, stage_params[i])
        end
        S = typeof(stages)
        scenarioproblems = ntuple(Val(M)) do i
            ScenarioProblems(scenarios[i], procs)
        end
        SP = typeof(scenarioproblems)
        return new{N, M, S, SP}(stages, scenarioproblems, Dict{Symbol, Function}(), Dict{Symbol, JuMP.Model}(), SPSolver(solver))
    end
end
TwoStageStochasticProgram{S <: Tuple{Stage, Stage}, SP <: AbstractScenarioProblems} = StochasticProgram{2, 1, S, SP}

# Constructors #
# ========================== #
# Two-stage
# ========================== #
"""
    StochasticProgram(first_stage_params::Any,
                      second_stage_params::Any,
                      solver = JuMP.UnsetSolver(),
                      procs = workers())

Create a new two-stage stochastic program with stage data given by `first_stage_params` and `second_stage_params`. After construction, scenarios of type `Scenario` can be added through `add_scenario!`. Optionally, a capable `solver` can be supplied to later optimize the stochastic program. If multiple Julia processes are available, the resulting stochastic program will automatically be memory-distributed on these processes. This can be avoided by setting `procs = [1]`.
"""
function StochasticProgram(first_stage_params::Any, second_stage_params::Any; solver = JuMP.UnsetSolver(), procs = workers()) where S <: AbstractScenario
    return StochasticProgram(first_stage_params, second_stage_params, Scenario, solver, procs)
end
"""
    StochasticProgram(first_stage_params::Any,
                      second_stage_params::Any,
                      ::Type{S};
                      solver = JuMP.UnsetSolver(),
                      procs = workers()) where S <: AbstractScenario

Create a new two-stage stochastic program with stage data given by `first_stage_params` and `second_stage_params`. After construction, scenarios of type `S` can be added through `add_scenario!`. Optionally, a capable `solver` can be supplied to later optimize the stochastic program. If multiple Julia processes are available, the resulting stochastic program will automatically be memory-distributed on these processes. This can be avoided by setting `procs = [1]`.
"""
function StochasticProgram(first_stage_params::Any, second_stage_params::Any, ::Type{S}; solver = JuMP.UnsetSolver(), procs = workers()) where S <: AbstractScenario
    return StochasticProgram(first_stage_params, second_stage_params, S, solver, procs)
end
"""
    StochasticProgram(; solver = JuMP.UnsetSolver(),
                      procs = workers())

Create a new two-stage stochastic program with scenarios of type `S` and no stage data.
"""
StochasticProgram(; solver = JuMP.UnsetSolver(), procs = workers()) where S <: AbstractScenario = StochasticProgram(nothing, nothing, Scenario; solver = solver, procs = procs)
"""
    StochasticProgram(::Type{S};
                      solver = JuMP.UnsetSolver(),
                      procs = workers()) where S <: AbstractScenario

Create a new two-stage stochastic program with scenarios of type `S` and no stage data.
"""
StochasticProgram(::Type{S}; solver = JuMP.UnsetSolver(), procs = workers()) where S <: AbstractScenario = StochasticProgram(nothing, nothing, S; solver = solver, procs = procs)
"""
    StochasticProgram(first_stage_params::Any,
                      second_stage_params::Any,
                      scenarios::Vector{<:AbstractScenario};
                      solver = JuMP.UnsetSolver(),
                      procs = workers())

Create a new two-stage stochastic program with a given collection of `scenarios`
"""
function StochasticProgram(first_stage_params::Any, second_stage_params::Any, scenarios::Vector{<:AbstractScenario}; solver = JuMP.UnsetSolver(), procs = workers())
    return StochasticProgram(first_stage_params, second_stage_params, scenarios, solver, procs)
end
"""
    StochasticProgram(scenarios::Vector{<:AbstractScenario};
                      solver = JuMP.UnsetSolver(),
                      procs = workers()) where {SD <: AbstractScenario}

Create a new two-stage stochastic program with a given collection of `scenarios` and no stage data.
"""
StochasticProgram(scenarios::Vector{<:AbstractScenario}; solver = JuMP.UnsetSolver(), procs = workers()) = StochasticProgram(nothing, nothing, scenarios; solver = solver, procs = procs)

# Printing #
# ========================== #
function Base.show(io::IO, stochasticprogram::StochasticProgram{N}) where N
    plural(n) = (n==1 ? "" : "s")
    defer_stage(sp, s) = begin
        if deferred_stage(sp,s)
            return " * deferred stage $s"
        else
            xdim = decision_length(sp, s)
            n = nscenarios(sp, s)
            if s == 1
                if N == 2
                    return " * $(xdim) decision variable$(plural(xdim))"
                else
                    return " * Stage $s:\n   * $(xdim) decision variable$(plural(xdim))"
                end
            elseif s == 2 && N == 2
                stype = typename(scenariotype(stochasticprogram))
                return " * $(xdim) recourse variable$(plural(xdim))\n * $(n) scenario$(plural(n)) of type $stype"
            else
                stype = typename(scenariotype(stochasticprogram, s))
                if distributed(stochasticprogram, s)
                    return " * Distributed stage $s:\n   * $(xdim) decision variable$(plural(xdim))\n   * $(n) scenario$(plural(n)) of type $stype"
                else
                    return " * Stage $s:\n   * $(xdim) decision variable$(plural(xdim))\n   * $(n) scenario$(plural(n)) of type $stype"
                end
            end
        end
    end
    stage(sp, s) = begin
        stage_key = Symbol(:stage_,s)
        if has_generator(sp, stage_key)
            return defer_stage(sp,s)
        else
            " * undefined stage $s"
        end
    end
    if N == 2 && distributed(stochasticprogram)
        println(io, "Distributed stochastic program with:")
    else
        println(io, "Stochastic program with:")
    end
    println(io, stage(stochasticprogram, 1))
    for s = 2:N
        println(io, stage(stochasticprogram, s))
    end
    print(io, "Solver is ")
    if isa(spsolver(stochasticprogram), JuMP.UnsetSolver)
        print(io, "default solver")
    else
        print(io, solverstr(spsolver(stochasticprogram)))
    end
end
function Base.print(io::IO, stochasticprogram::StochasticProgram{N}) where N
    print(io, "Stage 1\n")
    print(io, "============== \n")
    print(io, get_stage_one(stochasticprogram))
    for s = 2:N
        print(io, "\nStage $s\n")
        print(io, "============== \n")
        for (id, subproblem) in enumerate(subproblems(stochasticprogram, s))
            @printf(io, "Subproblem %d (p = %.2f):\n", id, probability(scenario(stochasticprogram, id, s)))
            print(io, subproblem)
            print(io, "\n")
        end
    end
end
function Base.print(io::IO, stochasticprogram::StochasticProgram{2})
    print(io, "First-stage \n")
    print(io, "============== \n")
    print(io, get_stage_one(stochasticprogram))
    print(io, "\nSecond-stage \n")
    print(io, "============== \n")
    for (id, subproblem) in enumerate(subproblems(stochasticprogram))
        @printf(io, "Subproblem %d (p = %.2f):\n", id, probability(scenario(stochasticprogram, id)))
        print(io, subproblem)
        print(io, "\n")
    end
end
# ========================== #


text:
using StatsBase

"""
  use the Loyd-Max algorithm to find a quantizer with `M` quantization levels
  for the signal `sig`.
"""
function loydMax(sig::Vector{T}, M::Int64; iterations=100) where T<:Real
  # sort signal to simplify the algorithm
  sig_sort = sort(sig)

  # start with uniform levels
  x = collect(range(minimum(sig), stop=maximum(sig),length=M))
  for i=1:iterations
    # compute decision thresholds
    t = 0.5*(x[1:end-1].+x[2:end])
    # @info "t = $(t)"
    # update representative levels
    num = 0
    denom = 0
    idx = 1
    for j=1:length(sig)
      if idx<M && sig_sort[j]<t[idx]
        num += sig_sort[j]
        denom += 1
      elseif idx<M
        x[idx] = num/denom
        num = 0
        denom = 0
        idx += 1
      else
        num += sig_sort[j]
        denom += 1
      end
      x[M] = num/denom
    end
  end

  return x
end


text:
# Copyright (c) 2021
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

using LinearAlgebra, Random;

"""
    SL0(n::Int64, p::Int64, k::Int64, type::String="Gaussian")

Generate a dataset for compressed sensing.

### Input

- `n       -- Int64: number of rows
- `p`      -- Int64: number of columns
- `k`      -- Int64: number true non-zeros
- `type`   -- String: type to generate ["Gaussian"]

### Output

(A, x, y) for Ax=y
"""
function cs_model(n::Int64, p::Int64, k::Int64, type::String="Gaussian")
    local A, x;
    if type == "Gaussian"
        A = randn(n, p);
        while rank(A) != n
            A = randn(n, p);
        end
        x = zeros(p);
        pp = sign.(randn(k)).*(ones(k)+abs.(randn(k)));
        rp = randperm(p);
        x[rp[1:k]] = pp;
        b = A*x;
    else
        error("Uknown type in cs_model(n,m,type).")
    end

    return A, x, b;
end


text:
module MixtureModels

using Distributions, StatsBase

export fit_em, emstep

include("./univariate_discrete.jl")
include("./univariate_continuous.jl")
include("./fit_em.jl")

emstep(comp::Distribution, samples::AbstractVector; kwargs...) = emstep(comp, samples, uweights(length(samples)); kwargs...)
emstep(comp::Distribution, samples::AbstractVector, weights::AbstractWeights; kwargs...) = emstep(comp, samples, weights, uweights(length(weights)); kwargs...)

end


text:
using Plots
include("N_plots.jl")

# Model and Objective
model, obj = Dynamics.dubinscar_parallelpark
obj_uncon = UnconstrainedObjective(obj)
obj = update_objective(obj,tf=2.,u_min=-2,u_max=2)

# Solver Options
opts = SolverOptions()
opts.verbose = false
opts.cost_tolerance = 1e-6
opts.cost_tolerance_intermediate = 1e-5
opts.constraint_tolerance = 1e-5
opts.resolve_feasible = false
opts.outer_loop_update_type = :default
opts.use_nesterov = true
opts.penalty_scaling = 200
opts.penalty_initial = .1
opts.R_infeasible = 20

# Solve original problem
N = 51
solver = Solver(model, obj, N=N, opts=opts)
n,m,N = get_sizes(solver)
U0 = ones(m,N-1)
res, stats = solve(solver,U0)

# Mintime
opts = SolverOptions()
opts.verbose = false
opts.cost_tolerance = 1e-6
opts.cost_tolerance_intermediate = 5e-5
opts.max_dt = 0.2
opts.min_dt = 1e-3
opts.minimum_time_dt_estimate = obj.tf/(N-1)
opts.constraint_tolerance = 0.0001 # 0.005
opts.R_minimum_time = 2. #15.0 #13.5 # 12.0
opts.constraint_decrease_ratio = .25
opts.penalty_scaling = 10.0
opts.outer_loop_update_type = :individual
opts.iterations = 1000
opts.iterations_outerloop = 30 # 20
opts.square_root = true

ipopt_options = Dict("tol"=>opts.cost_tolerance,"constr_viol_tol"=>opts.constraint_tolerance)

obj_mintime = update_objective(obj,tf=:min)

solver_mintime = Solver(model, obj_mintime, N=N, opts=opts)
results_mintime, stats_mintime = solve(solver_mintime,to_array(res.U))
stats_mintime["iterations"]
evals(solver_mintime,:f) / stats_mintime["iterations"]

plot(to_array(results_mintime.U)[1:2,1:solver_mintime.N-1]',linestyle=:solid,color=[1 2],labels=["v" "omega"],linewidth=2)

X0 = rollout(solver,U0)
res_d, stats_d = solve_dircol(solver_mintime,X0,U0,options=ipopt_options)
solve_dircol(solver_mintime,X0,U0,options=ipopt_options)
stats_d["info"]
stats_d["iterations"]
evals(solver_mintime,:f) / stats_d["iterations"]

T = TrajectoryOptimization.total_time(solver,res)
T_min = TrajectoryOptimization.total_time(solver_mintime,results_mintime)
T_d = TrajectoryOptimization.total_time(solver_mintime,res_d)

x = [get_time(solver),get_time(solver_mintime,results_mintime),get_time(solver_mintime,results_mintime)]
# x = [1:N for i= 1:3]
p1 = plot(x[1],to_array(res.U)[1,1:solver_mintime.N-1],labels=["ALTRO (init)" ""],linewidth=2,linestyle=[:solid],color=[:black :black],
    xlabel="Time (s)",ylabel="Linear velocity")
plot!(x[3],res_d.U[1,1:solver_mintime.N-1],label=["DIRCOL (mintime)" ""],linestyle=[:solid],color=[:blue :blue],linewidth=1.5)
plot!(x[2],to_array(results_mintime.U)[1,1:solver_mintime.N-1],label=["ALTRO (mintime)" ""],
    linestyle=[:solid],color=[:darkorange2 :darkorange2],linewidth=2,legend=:none)

p2 = plot(x[1],to_array(res.U)[2,1:solver_mintime.N-1],labels=["ALTRO (init)" ""],linewidth=2,linestyle=[:solid],color=[:black :black],
    xlabel="Time (s)",ylabel="Angular velocity")
plot!(x[3],res_d.U[2,1:solver_mintime.N-1],label=["DIRCOL" ""],linestyle=[:solid],color=[:blue :blue],linewidth=2,legend=:topright)
plot!(x[2],to_array(results_mintime.U)[2,1:solver_mintime.N-1],label=["ALTRO" ""],linestyle=[:solid],color=[:darkorange2 :darkorange2],linewidth=2)
plot(p1,p2,layout=(1,2),size=(500,250),dpi=400)
savefig(joinpath(IMAGE_DIR,"ppark_mintime_control.png"))

plot(to_array(res.X)[1:3,:]',linewidth=2,linestyle=:dash,xlabel="Knot point",ylabel="state",labels=["x" "y" "theta"])
plot!(to_array(results_mintime.X)[1:3,:]',linewidth=2,color=[1 2 3],legend=:topleft,labels="")
plot!(res_d.X[1:3,:]',linewidth=2,color=[1 2 3],legend=:topleft,labels="",style=:dot,dpi=400)
savefig(joinpath(IMAGE_DIR,"ppark_mintime_state.png"))

t = get_time(solver)
t_mintime = get_time(solver_mintime,results_mintime)
plot(t,to_array(res.X)[1:3,:]',linewidth=2,linestyle=:dash,xlabel="time (s)",ylabel="state",labels=["x" "y" "theta"])
plot!(t_mintime,to_array(results_mintime.X)[1:3,1:N-1]',linewidth=2,color=[1 2 3],legend=:topleft,labels="")

p = plot()
plot_trajectory!(res,label="original",width=2,xlabel="x",ylabel="y")
plot_trajectory!(results_mintime,labels="mintime",width=2,legend=:topleft)
plot_trajectory!(res_d,labels="dircol",width=2,legend=:topleft)

# Timing results
U0 = to_array(res.U)
@btime solve($solver_mintime,$U0)
res0, stats0 = solve(solver_mintime,U0)
evals(solver_mintime,:f) / stats0["iterations"]
runtimes = []
iters = Int[]
res_d, stats_d = solve_dircol(solver_mintime,X0,U0,options=ipopt_options)
stats_d["info"]
stats_d["iterations"]
evals(solver_mintime,:f) / stats_d["iterations"]

push!(iters,stats_d["iterations"])
push!(runtimes,stats_d["runtime"])
mean(runtimes)
pop!(iters)
iters


text:
#=
Compute the running median of a sequence of numbers. That is, given a stream of numbers, print out the median of the list so far on each new element.

Recall that the median of an even-numbered list is the average of the two middle numbers.

For example, given the sequence [2, 1, 5, 7, 2, 0, 5], your algorithm should print out:

2
1.5
2
3.5
2
2
2
=#

using Test
include("Solutions/problem33_running_median.jl")

@test print_running_median([2, 1, 5, 7, 2, 0, 5], store = true) == [2, 1.5, 2, 3.5, 2, 2.0, 2]

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
# Axis-aligned bounding box.
# An N-dimensional box requires 2 N-dimensional points to specify the boundary:
#   One point to specify the box origin, and on to specify the opposite (furthest corner)
struct AABB{N,T} <: Face{N,T}
    origin::Point{N,T}
    corner::Point{N,T}
end

const AABB_2D = AABB{2}
const AABB_3D = AABB{3}

# Constructors & Conversions
# -------------------------------------------------------------------------------------------------
AABB(p₁::Point{N,T}, p₂::Point{N,T}) where {N,T} = AABB{N,T}(p₁, p₂)
AABB{N}(p₁::Point{N,T}, p₂::Point{N,T}) where {N,T} = AABB{N,T}(p₁, p₂)

# Base
# -------------------------------------------------------------------------------------------------
# Note: all branches but the coraabb one are pruned by the compiler
function Base.getproperty(aabb::AABB, sym::Symbol)
    if sym === :xmin
        return aabb.origin[1]
    elseif sym === :ymin
        return aabb.origin[2]
    elseif sym === :zmin
        return aabb.origin[3]
    elseif sym === :xmax
        return aabb.corner[1]
    elseif sym === :ymax
        return aabb.corner[2]
    elseif sym === :zmax
        return aabb.corner[3]
    else # fallback to getfield
        return getfield(aabb, sym)
    end
end

# Methods
# -------------------------------------------------------------------------------------------------
@inline width(aabb::AABB) = aabb.xmax - aabb.xmin
@inline height(aabb::AABB) = aabb.ymax - aabb.ymin
@inline depth(aabb::AABB) = aabb.ymax - aabb.ymin
@inline area(aabb::AABB_2D) = height(aabb) * width(aabb)
@inline volume(aabb::AABB_3D) = height(aabb) * width(aabb) * depth(aabb)
@inline Base.in(p::Point_2D, aabb::AABB_2D) = aabb.xmin ≤ p[1] ≤ aabb.xmax &&
                                              aabb.ymin ≤ p[2] ≤ aabb.ymax
@inline Base.in(p::Point_3D, aabb::AABB_3D) = aabb.xmin ≤ p[1] ≤ aabb.xmax &&
                                              aabb.ymin ≤ p[2] ≤ aabb.ymax &&
                                              aabb.zmin ≤ p[3] ≤ aabb.zmax

# DEPRECATED. Leaving because there is potential future use
# # Liang-Barsky line clipping algorithm
# # pₖ = 0	            parallel to the clipping boundaries
# # pₖ = 0 and qₖ < 0	    completely outside the boundary
# # pₖ = 0 and qₖ ≥ 0	    inside the parallel clipping boundary
# # pₖ < 0	            line proceeds from outside to inside
# # pₖ > 0	            line proceeds from inside to outside
# function intersect(l::LineSegment_2D{F}, aabb::AABB_2D{F}) where {F <: AbstractFloat}
#     p₂ = l[2].x - l[1].x
#     p₁ = -p₂
#     p₄ = l[2].y - l[1].y
#     p₃ = -p₄
#
#     q₁ = l[1].x - aabb.xmin
#     q₂ = aabb.xmax - l[1].x
#     q₃ = l[1].y - aabb.ymin
#     q₄ = aabb.ymax - l[1].y
#
#     # Line parallel to clipping window
#     if p₁ == 0 # Vertical line
#         if q₁ < 0 || q₂ < 0 # Outside boundaries
#             return false, SVector(Point_2D{F}(0, 0), Point_2D{F}(0, 0))
#         else # Inside boundaries
#             return true, SVector(Point_2D(l[1].x, aabb.ymin), Point_2D(l[1].x, aabb.ymax))
#         end
#     end
#     if p₃ == 0 # Horizontal line
#         if q₃ < 0 || q₄ < 0 # Outside boundaries
#             return false, SVector(Point_2D{F}(0, 0), Point_2D{F}(0, 0))
#         else # Inside boundaries
#             return true, SVector(Point_2D(aabb.xmin, l[1].y), Point_2D(aabb.xmax, l[1].y))
#         end
#     end
#
#     t₁ = q₁ / p₁
#     t₂ = q₂ / p₂
#     if (p₁ < 0)
#         t_min2 = t₁
#         t_max2 = t₂
#     else
#         t_min2 = t₂
#         t_max2 = t₁
#     end
#
#     t₃ = q₃ / p₃
#     t₄ = q₄ / p₄
#     if (p₃ < 0)
#         t_min3 = t₃
#         t_max3 = t₄
#     else
#         t_min3 = t₄
#         t_max3 = t₃
#     end
#
#     t_start = max(F(0), t_min2, t_min3)
#     t_stop = min(F(1), t_max2, t_max3)
#
#     # Line outside clipping window
#     t_start < t_stop || return false, SVector(Point_2D{F}(0, 0), Point_2D{F}(0, 0))
#
#     return true, SVector(l(t_start), l(t_stop))
# end
#
# Credit to Tavian Barnes (https://tavianator.com/2011/ray_box.html)
# Assumes the line passes all the way through the AABB if it intersects, which is a
# valid assumption for this ray tracing application.
#
# Note that Liang-Barsky (above) performs better on a single-threaded CPU than this
# algorithm. But, the multithreaded performance is approximately the same, and since this is
# algorithm is branchless it tends to perform better than LB on GPUs.
function Base.intersect(l::LineSegment_2D, aabb::AABB_2D)
    𝐮 = l[2] - l[1]
    𝐭₁ = (aabb.origin - l[1]) ./ 𝐮
    𝐭₂ = (aabb.corner - l[1]) ./ 𝐮

    tmin = max(min(𝐭₁[1], 𝐭₂[1]), min(𝐭₁[2], 𝐭₂[2]))
    tmax = min(max(𝐭₁[1], 𝐭₂[1]), max(𝐭₁[2], 𝐭₂[2]))

    return (tmax >= tmin, SVector(l(tmin), l(tmax)))
end

# A random AABB within [0, 1]ᴺ
# What does the distribution of AABBs look like? Is this uniform?
function Base.rand(::Type{AABB{N,T}}) where {N,T}
    coord₁ = rand(T, N)
    coord₂ = rand(T, N)
    return AABB{N,T}(Point{N,T}(min.(coord₁, coord₂)),
                     Point{N,T}(max.(coord₁, coord₂)))
end

# NB random AABB within [0, 1]ᴺ
function Base.rand(::Type{AABB{N,T}}, NB::Int64) where {N,T}
    return [ rand(AABB{N,T}) for i ∈ 1:NB ]
end

# Find the AABB which contains both bb₁ and bb₂
function Base.union(bb₁::AABB{N,T}, bb₂::AABB{N,T}) where {N,T}
    return AABB(Point{N,T}(min.(bb₁.origin, bb₂.origin)),
                Point{N,T}(max.(bb₁.corner, bb₂.corner)))
end

# # Plot
# # -------------------------------------------------------------------------------------------------
# if enable_visualization
#     function convert_arguments(LS::Type{<:LineSegments}, aabb::aabbangle_2D)
#         p₂ = Point_2D(aabb.xmax, aabb.ymin)
#         p₄ = Point_2D(aabb.xmin, aabb.ymax)
#         l₁ = LineSegment_2D(aabb.bl, p₂)
#         l₂ = LineSegment_2D(p₂, aabb.tr)
#         l₃ = LineSegment_2D(aabb.tr, p₄)
#         l₄ = LineSegment_2D(p₄, aabb.bl)
#         lines = [l₁, l₂, l₃, l₄]
#         return convert_arguments(LS, lines)
#     end
#
#     function convert_arguments(LS::Type{<:LineSegments}, R::Vector{<:aabbangle_2D})
#         point_sets = [convert_arguments(LS, aabb) for aabb in R]
#         return convert_arguments(LS, reduce(vcat, [pset[1] for pset ∈ point_sets]))
#     end
#
#     function convert_arguments(M::Type{<:Mesh}, aabb::aabbangle_2D)
#         p₂ = Point_2D(aabb.xmax, aabb.ymin)
#         p₄ = Point_2D(aabb.xmin, aabb.ymax)
#         points = [aabb.bl, p₂, aabb.tr, p₄]
#         faces = [1 2 3;
#                  3 4 1]
#         return convert_arguments(M, points, faces)
#     end
#
#     function convert_arguments(M::Type{<:Mesh}, R::Vector{<:aabbangle_2D})
#         points = reduce(vcat, [[aabb.bl, Point_2D(aabb.xmax, aabb.ymin),
#                                 aabb.tr, Point_2D(aabb.xmin, aabb.ymax)] for aabb ∈ R])
#         faces = zeros(Int64, 2*length(R), 3)
#         j = 0
#         for i in 1:2:2*length(R)
#             faces[i    , :] = [1 2 3] .+ j
#             faces[i + 1, :] = [3 4 1] .+ j
#             j += 4
#         end
#         return convert_arguments(M, points, faces)
#     end
# end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
@deprecate writeMIDIFile(filename::AbstractString, data::MIDIFile) save(File{format"MIDI"}(filename), data)
@deprecate writeMIDIFile(filename::AbstractString, notes::Notes) save(File{format"MIDI"}(filename), notes)
@deprecate readMIDIFile(filename::AbstractString) load(File{format"MIDI"}(filename))


text:
using Gen
using Test

include("regression/static_mala.jl")
include("regression/static_map_optimize.jl")
include("regression/static_mh.jl")
include("regression/static_mh_collapsed.jl")
include("regression/quickstart.jl")
include("regression/dynamic_map_optimize_gibbs.jl")
include("regression/dynamic_mh.jl")
include("regression/dynamic_resimulation_mh.jl")
include("coal/coal.jl")
include("gp_structure/lightweight.jl")
include("gp_structure/incremental.jl")
include("mle/static_mle.jl")

max_stars_repo_name: dendisuhubdy/Gen

max_issues_repo_name: dendisuhubdy/Gen

max_forks_repo_name: dendisuhubdy/Gen


text:
#L2HingeLoss

using LaTeXStrings, LossFunctions, Measures, Latexify

Plots.plot(L2HingeLoss(), c=:olive, leg=:topright, linewidth=2, annotation=(1, 6, L"L\left( a \right) = max\{0,1, - a \}^2", 18))

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
data = vcat([10*randn(2,1).+randn(2,1000) for i in 1:10]...)
@time a = quickshift(data)
Profile.clear_malloc_data()
@time a = quickshift(data)
exit()


text:
@testset "Tropical Curve" begin
    @polyvar x y

    g = Tropical(0) + x^3 + y^3 + log(2)*x*y
    t = TropicalCurve(g)
    n3 = [0.0, -log(2)]
    n2 = [-log(2), 0.0]
    n1 = [log(2), log(2)]
    @test vertices(t)[1] ≈ n1
    @test vertices(t)[2] ≈ n3
    @test vertices(t)[3] ≈ n2

    @test segments(t)[1][1] ≈ n1
    @test segments(t)[1][2] ≈ n3
    @test segments(t)[2][1] ≈ n1
    @test segments(t)[2][2] ≈ n2
    @test segments(t)[3][1] ≈ n3
    @test segments(t)[3][2] ≈ n2

    @test length(halfrays(t)) == 3

    # c = ArchTrop(1 + x + y + x * y + y^2 + x^2 * y).curve
    # @test length(halfrays(c)) == 4
    # @test length(vertices(c)) == 1
end


text:
using HMCExamples

HMCExamples.main_FVGQ_1_joint()


text:
using LinearAlgebra, BenchmarkTools
using Distributed, SharedArrays

@everywhere N = 1000
@everywhere n = 1000
@everywhere W = [rand(n, n) for i = 1:N]
@everywhere A = [W[i]' * W[i] for i = 1:N]
@everywhere b = [rand(n) for i = 1:N]

# Distributed
addprocs(4)
procs()
workers()
nworkers()
myid()
#
# remotecall_fetch(() -> myid(), 3)
# pmap(i -> println("I'm worker $(myid()), working on i=$i"), 1:10)
# @sync @distributed for i in 1:10
#   println("I'm worker $(myid()), working on i=$i")
# end
#
# a = Array{Int}(undef, 10)
# a
# @sync @distributed for i in 1:10
#   println("working on i=$i")
#   a[i] = i^2
# end
# a
#

# @everywhere printsquare(i) = println("working on i=$i: its square it $(i^2)")
# @sync @distributed for i in 1:10
#   printsquare(i)
# end

function solve_dist(x, A, b, N)
    @sync @distributed for i = 1:N
        x[i] = A[i] \ b[i]
    end
end

@everywhere x = SharedArray{Float64}((n, N))

print("run distributed")
solve_dist(x, A, b, N)


text:
module fastICA
    export correlation_matrix, covariance_matrix, eigen_decomposition, sing_val_decomp, whiten, fast_ica
    using LinearAlgebra
    using Statistics

    correlation_matrix(mat::Array{Float64 ,2})::Array{Float64 ,2} = (mat .- mean(mat,dims=1)) ./ std(mat,dims=1)

    covariance_matrix(corrmat::Array{Float64 ,2})::Array{Float64 ,2} = corrmat |> length |> float |> n -> corrmat' * corrmat / (n  - 1.)

    sing_val_decomp(mat::Array{Float64 ,2})::SVD = mat |> correlation_matrix |> svd

    #Returns: a whitened matrix using one of the available whitening techniques(PCA or ZCA)
    #Receives:
    #mat -> the matrix to whiten
    #pca -> optional boolean which is used to determine the whitening technique
    function whiten(mat::Array{Float64 ,2},pca::Bool = true)::Array{Float64 ,2}
        sigma = cov(mat)
        vals,vecs = eigen(sigma)
        #Sort eigvals and eigvecs by DESCENDING order
        vecs = reverse(vecs, dims=2)
        vals = reverse(vals, dims=1)
        pca_w =  Diagonal(1 ./ sqrt.(vals)) * vecs'
        #U,S,_ = svd(sigma)
        #pca_w =  Diagonal(1 ./ sqrt.(S)) * U'
        if(pca)
            return pca_w
        end
        #zca whitening
        return vecs * pca_w
    end

    #auxiliary function for logcosh contrast function
    function contrast_func(alpha::Float64,wx)::Tuple{Array{Float64},Array{Float64}}
       gd= tanh.(alpha * wx)
       gdd= alpha * (1 .- tanh.(alpha * wx).^2)
       return (gd,gdd)
    end

    #Fast Ica deflation algorithm
    function fast_ica_def(maxiter::Int64, nic::Int64,X::Array{Float64,2},tol::Float64,W::Array{Float64,2},alpha::Float64)::Array{Float64 , 2}
        #return W
        retW = zeros(size(W))
        #src : http://www.measurement.sk/2011/Patil.pdf pg 119
        for i = 1:nic
            println("Component n $i")
            wp = W[i,:,]
            #to-do create aux func
            if (i > 1)
                t = zeros(size(wp))
                for u = 1:(i-1)
                   k = sum(wp .* retW[u,:,])
                   t = t + k * retW[u,:,]
                end
                wp = wp - t
            end
            normalize!(wp)

            iter = 0
            chg = 0
            converge = false
            while !converge && iter < maxiter
                iter+=1
                println("iter = $iter")
                wx = wp' * X
                gdx,gddx = contrast_func(alpha,wx)
                xgdx = X .* gdx
                v1 = vec(mapslices(mean, xgdx, dims = 2))
                println(v1)
                v2 = mean(gddx) * wp
                w1 = v1 - v2
                #to-do create aux func
                if (i > 1)
                    t = zeros(size(w1))
                    for u = 1:(i-1)
                        k = sum(w1 .* retW[u,:,])
                        t = t + k * retW[u,:,]
                    end
                    w1 = w1 - t
                end
                normalize!(w1)
                #check for convergence
                chg =  abs.(abs.(sum(w1 .* wp)) .- 1.0)
                println("Tolerance change for iter $iter = $chg")
                wp = w1
                converge = ( chg < tol )
            end
            retW[i,:,] = wp
        end
        println("retW = $retW")
        return retW
    end

   #Returns: the resultant ICA model, an instance of the type ICA
    #Receives:
    #maxiter -> number of iterations for the main loop
    #nic     -> number of independent components
    #X       -> the data matrix, must be whitened(use the whiten function)
    #tol     -> tolerable change of the weights at convergence
    #
    #optional :
    #alpha   -> must be between [1,2]
    #whiten  -> whiten the data
    function fast_ica(maxiter::Int64, nic::Int64,X::Array{Float64,2},tol::Float64, alpha::Float64 = 1. )::Array{Float64,2}
        #validate arguments
        m,n = size(X)
        m > 1 || error("There must be more than one samples, n > 1.")
        maxiter > 1 || error("maxiter must be greater than 1.")
        tol > 0 || error("tol must be positive.")
        alpha >= 1 && alpha <= 2|| error("alpha must be in between 1 and 2")
        comp = min(n,m)
        if(nic > comp)
          nic = comp
        end
        # initialize weights of size n with random values
        #src : https://en.wikipedia.org/wiki/FastICA
        W = randn(comp,comp)

        println("Random matrix = $W")
        return fast_ica_def(maxiter,nic,X,tol,W,alpha)

    end

end

max_stars_repo_name: fica-ps/fica-case-studies

max_issues_repo_name: fica-ps/fica-case-studies

max_forks_repo_name: fica-ps/fica-case-studies


text:
@enum RotationVariables               RCardan123
@enum WStartVariables    WCartesian   WCardan123
@enum JointKind UndefinedJointKind FixKind FixTranslationKind RevoluteKind PrismaticKind FreeMotionKind AbsoluteFreeMotionKind

using Unitful

# FixedJoint is used if an Object3D is rigidly fixed to its parent Object3D
struct FixedJoint <: Modia3D.AbstractJoint
end

# FixedTranslationJoint is used if an Object3D is rigidly fixed to its parent Object3D without any rotation
struct FixedTranslationJoint <: Modia3D.AbstractJoint
end

"""
    joint = FreeMotion(; obj1, obj2, path="", r, rot, v, w)

Return a `joint` that describes the free movement of `obj2::`[`Object3D`](@ref)
with respect to `obj1::`[`Object3D`](@ref). The initial position is `r` and the
initial orientation is `rot` in [Cardan (Tait–Bryan) angles](https://en.wikipedia.org/wiki/Euler_angles#Chained_rotations_equivalence)
(rotation sequence x-y-z). `v` and `w` are the initial cartesian translational
and rotational velocity vectors.
"""
mutable struct FreeMotion <: Modia3D.AbstractJoint
    path::String

    obj1::Modia3D.AbstractObject3D
    obj2::Modia3D.AbstractObject3D

    ndof::Int

    r::SVector{3,Float64}
    rot::SVector{3,Float64}        # cardan angles (rotation sequence x-y-z)

    v::SVector{3,Float64}
    w::SVector{3,Float64}          # angular velocity vector

    a::SVector{3,Float64}
    z::SVector{3,Float64}          # angular acceleration vector

    residue_f::SVector{3,Float64}
    residue_t::SVector{3,Float64}

    function FreeMotion(; obj1::Modia3D.AbstractObject3D,
                          obj2::Modia3D.AbstractObject3D,
                          path::String = "",
                          r::AbstractVector   = Modia3D.ZeroVector3D,
                          rot::AbstractVector = Modia3D.ZeroVector3D,
                          v::AbstractVector   = Modia3D.ZeroVector3D,
                          w::AbstractVector   = Modia3D.ZeroVector3D)

        #(parent,obj,cutJoint) = attach(obj1, obj2)
        #if cutJoint
        #    error("\nError from FreeMotion joint connecting ", Modia3D.fullName(obj1), " with ", Modia3D.fullName(obj2), ":\n",
        #        "    This joint closes a kinematic loop which does not make sense -> remove this joint or correct the kinematic loop!")
        #end
        #
        #if !(obj === obj2)
        #    error("\nError from FreeMotion joint connecting obj1 = ", Modia3D.fullName(obj1), " with obj2 = ", Modia3D.fullName(obj2), ":\n",
        #        "    obj2 has a root, but obj1 has not root. This is currently not supported (exchange obj1 and obj2)!")
        #end

        r   = Modia3D.convertAndStripUnit(SVector{3,Float64}, u"m", r)
        rot = Modia3D.convertAndStripUnit(SVector{3,Float64}, u"rad", rot)
        v   = Modia3D.convertAndStripUnit(SVector{3,Float64}, u"m/s"  , v)
        w   = Modia3D.convertAndStripUnit(SVector{3,Float64}, u"rad/s", w)
        a   = Modia3D.ZeroVector3D
        z   = Modia3D.ZeroVector3D

        residue_f = Modia3D.ZeroVector3D
        residue_t = Modia3D.ZeroVector3D

        obj2.joint      = new(path, obj1, obj2, 6, r, rot, v, w, a, z, residue_f, residue_t)
        obj2.parent     = obj1
        obj2.jointKind  = FreeMotionKind
        obj2.jointIndex = 0
        obj2.ndof       = 6
        obj2.canCollide = true
        obj2.r_rel      = r
        obj2.R_rel      = Rfromrot123(rot)

        push!(obj1.children, obj2)

        return obj2.joint
    end
end


text:
using Random, Statistics
using Oceananigans

Nx = 32       # Number of grid points in x
Ny = 32       # Number of grid points in y
Nz = 32       # Number of grid points in z
Δz = 4.0      # Grid spacing in x, y, z (meters)

grid = RegularRectilinearGrid(size = (Nx, Ny, Nz), extent=(Δz*Nx, Δz*Ny, Δz*Nz),topology = (Periodic, Periodic, Bounded))

buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=2e-4, β=8e-4))

Qʰ = 200  # W m⁻², surface _heat_ flux
ρₒ = 1026 # kg m⁻³, average density at the surface of the world ocean
cᴾ = 3991 # J K⁻¹ s⁻¹, typical heat capacity for seawater

Qᵀ = Qʰ / (ρₒ * cᴾ) # K m⁻¹ s⁻¹, surface _temperature_ flux

dTdz = 0.01 # K m⁻¹

T_bcs = TracerBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵀ),  bottom = BoundaryCondition(Gradient, dTdz))

u₁₀ = 1.00    # m s⁻¹, average wind velocity 10 meters above the ocean
cᴰ = 2.5e-3  # dimensionless drag coefficient
ρₐ = 1.225   # kg m⁻³, average density of air at sea-level

Qᵘ = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) # m² s⁻²

u_bcs = UVelocityBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵘ))

@inline Qˢ(x, y, t, S, evaporation_rate) = - evaporation_rate * S

evaporation_rate = 1e-3 / 3600

evaporation_bc = BoundaryCondition(Flux, Qˢ, field_dependencies=:S, parameters=evaporation_rate)

S_bcs = TracerBoundaryConditions(grid, top=evaporation_bc)

model = IncompressibleModel(architecture = CPU(),
                            advection = UpwindBiasedFifthOrder(),
                            timestepper = :RungeKutta3,
                            grid = grid,
                            coriolis = FPlane(f=1e-4),
                            buoyancy = buoyancy,
                            closure = AnisotropicMinimumDissipation(),
                            boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs))

Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise

Tᵢ(x, y, z) = 20 + dTdz * z + dTdz * model.grid.Lz * 1e-6 * Ξ(z)

uᵢ(x, y, z) = sqrt(abs(Qᵘ)) * 1e-3 * Ξ(z)

set!(model, u=uᵢ, w=uᵢ, T=Tᵢ, S=35)

wizard = TimeStepWizard(cfl=1.0, Δt=1.0, max_change=1.1, max_Δt=20.0)

simulation = Simulation(model, iteration_interval = 10, stop_time=60*60*2, Δt = wizard)
run!(simulation)

Nsimulation = Simulation(model, iteration_interval = 10, stop_time=60*60*3, Δt = 20.0)

Nsimulation.output_writers[:fields] =
        JLD2OutputWriter(model, model.velocities,
                         schedule = TimeInterval(60),
                         dir = ".",
                         prefix = "velocities",
                         force = true)

run!(Nsimulation)


text:
using JLD
using Base.Test
import Base.convert
using MAT
using ForwardDiff
using Polynomials
import ForwardDiff.DiffBase
using Optim
# using Pandas
using JSON

# @everywhere path="/home/hcli64/hcli64751/rats/fit_behaviour_julia/"
@everywhere path="/home/emma/github/fit_behaviour_julia/"

#
@everywhere include(path*"functions/functions_julia.jl")
@everywhere include(path*"functions/parameter_functions.jl")
@everywhere include(path*"functions/utils.jl")
@everywhere rat_num = "1"
@everywhere path_load=path*"scripts_rats/data/"
@everywhere filename="processed_data_rat"*rat_num*".json"

# local
# @everywhere path="/home/genis/fit_behaviour_julia/"
#######bsc#########
#
# @everywhere path_functions=path*"scripts/functions/"
# @everywhere path_data=path*"data/"
# @everywhere path_results=path*"results/"
#
# @everywhere include(path_functions*"parameter_functions.jl")
# @everywhere include(path_functions*"functions_julia.jl")
# @everywhere include(path_functions*"utils.jl")

Nframes=10
Tframe=.05

model_num="6"

@everywhere args=[("c2",[1]),("c4",0),("st",[1,2]),("x0",0),("bias_d",0),("hbias_d",0),("sigma_a",[1])]
@everywhere consts=[("c2",[2]),("c4",[1]),("st",[3,4]),("x0",[1,2]),("bias_d",[1]),("hbias_d",[1]),("sigma_a",0)]

@everywhere z = [0,1,0,0,0,0,0,0]

@everywhere lb=[-10.,0.0,0.0,0.1]
@everywhere ub=[10.,10,10,10.]

@everywhere model_f=Dict(:c2_f=>c2_urgency,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
:hbias_d_f=>history_bias_d_const,:st_f=>st_logistic_bias2,:x0_f=>initial_condition_bias_hbias,
:ro_f=>read_out_perfect,:lambda=>0.0)

x_ini=[ 0.6,4.72,6.38,1.52]
@everywhere using JSON
@everywhere model_num="10"

@everywhere Nstim = 100

@everywhere data = Dict()
open(path_load*filename, "r") do f
    global data
    dicttxt = readstring(f)  # file information to string
    data=JSON.parse(dicttxt)  # parse and transform data
end

choices_10 = data["choices_10"][200+1:200+Nstim,:]
stim_10 = data["stim_10"][200+1:100+Nstim,:]
mu_vec = [-1.0, -0.9, -0.8, -0.7, -0.6, -0.55, -0.5, -0.45, -0.4, -0.3, -0.2, -0.15, -0.1, 0.0,
 0.1, 0.15, 0.2, 0.3, 0.4, 0.45, 0.5, 0.55, 0.6, 0.7, 0.8, 0.9, 1.0, -0.65, 0.65, -0.35, 0.35]
mu_10 = data["coherences_10"][200+1:100+Nstim,:]
sigma_10 = data["sigma_10"][200+1:100+Nstim,:]

@everywhere history_bias=zeros(Nstim)
###### check the gradient ##########

#
# println("computeLL")
# tic()
# PR=ProbRight2(stim,history_bias,Tframe,param2,model_f)
# toc()

println("computeLL2")
tic()
ll= ComputeLL2_mu(stim_10,sigma_10,mu_10,mu_vec,choices_10,history_bias,Tframe, args, x_ini,consts,z,model_f)
toc()

println("grads")
tic()
ll_1,grads=LL_function_grads_mu(stim_10,sigma_10,mu_10,mu_vec,choices_10,history_bias,Tframe, args, x_ini,consts,z,model_f)
toc()

grad_num=zeros(length(x))
hess_num=zeros(length(x))

dxx=0.001
println("grad num")

tic()
for iparam in 1:length(x)
        # println(iparam)
        y=x[:]
        # print(y[iparam]+dxx)
        # y[iparam]=y[iparam]+dxx
        lldx=LL_original=ComputeLL2_mu(stim_10,sigma_10,mu_10,mu_vec,choices_10,history_bias,Tframe, args, x_ini,consts,z,model_f)
        grad_num[iparam]=(lldx-ll)/dxx

        #y[iparam]=y[iparam]-2*dxx
        #lldx_m=LL_original=ComputeLL2(stim,choices,history_bias,Tframe, args, y,consts,z,model_f)
        #hess_num[iparam]=(lldx-2.*ll+lldx_m)/(dxx^2)
end
toc()

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
using VDPTag2Experiments
using Test

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


text:
using olcamat
using Base.Test

function test_matrixio()
    A = [1.0 2.0 3.0 ;
         4.0 5.0 6.0 ]
    path = joinpath(tempdir(), "_olcamat_jl_test.bin")
    olcamat.write_matrix(A, path)
    B = olcamat.read_matrix(path)
    @test A == B
    c1 = olcamat.read_column(path, 1)
    @test c1 == [1.0, 4.0]
    c2 = olcamat.read_column(path, 2)
    @test c2 == [2.0, 5.0]
    c3 = olcamat.read_column(path, 3)
    @test c3 == [3.0, 6.0]
    rm(path)
    @test isfile(path) == false
end

test_matrixio()


text:
using CSV # This is a pacakge we use for loading CSV Files.
using ScikitLearn
using DataFrames
@sk_import model_selection:train_test_split
@sk_import linear_model:LogisticRegression
@sk_import metrics:confusion_matrix
@sk_import metrics:classification_report
@sk_import model_selection:cross_val_score
@sk_import metrics:plot_confusion_matrix
@sk_import model_selection:KFold
pathtodata = joinpath("julia-scripts", "model-zoo", "covid_cleaned.csv")
data = DataFrame(CSV.File(pathtodata))

X = convert(Array, data[!, Not(:covid_res)])
y = convert(Array, data[!, :covid_res])

X_train, X_test, y_train, y_test = train_test_split(X, y; test_size=0.33, random_state=42)

simplelogistic = LogisticRegression(; max_iter=300)
fit!(simplelogistic, X_train, y_train)

y_pred = predict(simplelogistic, X_train)
print(classification_report(y_train, y_pred))
cross_val_score(LogisticRegression(; penalty=:none), X_train, y_train)
using PyPlot
plot_confusion_matrix(simplelogistic, X_train, y_train)
p1 = gcf()

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
export simplify

add_verbose_scope(:Simplify)

@doc Markdown.doc"""
    simplify(K::AnticNumberField; canonical::Bool = false) -> AnticNumberField, NfToNfMor

Tries to find an isomorphic field $L$ given by a "simpler" defining polynomial.
By default, "simple" is defined to be of smaller index, testing is done only
using a LLL-basis of the maximal order.

If `canonical` is set to `true`, then a canonical defining polynomial is found,
where canonical is using the definition of PARI's `polredabs`, which is described in
http://beta.lmfdb.org/knowledge/show/nf.polredabs.

Both version require a LLL reduced basis for the maximal order.
"""
function simplify(K::AnticNumberField; canonical::Bool = false, cached = false)
  Qx, x = PolynomialRing(FlintQQ)
  if degree(K) == 1
    L = number_field(x - 1, cached = cached, check = false)[1]
    return L, hom(L, K, gen(K), check = false)
  end
  if canonical
    a, f1 = polredabs(K)
    f = Qx(f1)
  else
    OK = maximal_order(K)
    @vtime :Simplify 3 if isdefined(OK, :lllO)
      ZK = OK.lllO
    else
      prec = 100 + 25*div(degree(K), 3) + Int(round(log(abs(discriminant(OK)))))
      ZK = _lll_for_simplify(OK, prec = prec)
    end
    @vtime :Simplify 3 a = _simplify(ZK)
    @vtime :Simplify 3 f = Qx(minpoly(representation_matrix(OK(a))))
  end
  L = NumberField(f, cached = cached, check = false)[1]
  m = hom(L, K, a, check = false)
  return L, m
end

function _simplify(O::NfOrd)
  K = nf(O)
  Qx, x = PolynomialRing(FlintQQ)
  Zx = PolynomialRing(FlintZZ, "x", cached = false)[1]
  f = Zx(K.pol*denominator(K.pol))

  a = O(gen(K)*denominator(K.pol), false)
  p, d = _find_prime(f)

  B = basis(O, copy = false)
  #First, we search for elements that are primitive using block systems
  F = FlintFiniteField(p, d, "w", cached = false)[1]
  Ft = PolynomialRing(F, "t", cached = false)[1]
  ap = zero(Ft)
  fit!(ap, degree(K)+1)
  rt = roots(f, F)

  n = degree(K)
  indices = Int[]
  for i = 1:length(B)
    b = _block(B[i].elem_in_nf, rt, ap)
    if length(b) == n
      push!(indices, i)
    end
  end
  #Now, we select the one of smallest T2 norm
  I = t2(a)
  for i = 1:length(indices)
    if isone(denominator(B[indices[i]].elem_in_nf))
      continue
    end
    t2n = t2(B[indices[i]].elem_in_nf)
    if t2n < I
      a = B[indices[i]]
      I = t2n
    end
  end
  return a.elem_in_nf
end

 #a block is a partition of 1:n
 #given by the subfield of parent(a) defined by a
 #the embeddings used are in R
 #K = parent(a)
 # then K has embeddings into the finite field (parent of R[1])
 # given by the roots (in R) of the minpoly of K
 #integers in 1:n are in the same block iff a(R[i]) == a(R[j])
 #the length of such a block (system) is the degree of Q(a):Q, the length
 # of a block is the degree K:Q(a)
 # a is primitive iff the block system has length n
function _block(a::nf_elem, R::Array{fq_nmod, 1}, ap::fq_nmod_poly)
  # TODO:
  # Maybe this _tmp business has to be moved out of this function too
  _R = GF(Int(characteristic(base_ring(ap))), cached = false)
  _Ry, _ = PolynomialRing(_R, "y", cached = false)
  _tmp = _Ry()
  nf_elem_to_gfp_poly!(_tmp, a, false) # ignore denominator
  set_length!(ap, length(_tmp))
  for i in 0:(length(_tmp) - 1)
    setcoeff!(ap, i, base_ring(ap)(_get_coeff_raw(_tmp, i)))
  end
#  ap = Ft(Zx(a*denominator(a)))
  s = fq_nmod[evaluate(ap, x) for x = R]
  b = Vector{Int}[]
  a = BitSet()
  i = 0
  n = length(R)
  while i < n
    i += 1
    if i in a
      continue
    end
    z = s[i]
    push!(b, findall(x->s[x] == z, 1:n))
    for j in b[end]
      push!(a, j)
    end
  end
  return b
end

#given 2 block systems b1, b2 for elements a1, a2, this computes the
#system for Q(a1, a2), the compositum of Q(a1) and Q(a2) as subfields of K
function _meet(b1::Vector{Vector{Int}}, b2::Vector{Vector{Int}})
  b = Vector{Int}[]
  for i=b1
    for j = i
      for h = b2
        if j in h
          s = intersect(i, h)
          if ! (s in b)
            push!(b, s)
          end
        end
      end
    end
  end
  return b
end

function _find_prime(f::fmpz_poly)
  p = 2^20
  d = 1
  while true
    p = next_prime(p)
    R = GF(p, cached=false)
    Rt = PolynomialRing(R, "t", cached = false)[1]
    fR = change_base_ring(R, f, parent = Rt)
    if degree(fR) != degree(f) || !issquarefree(fR)
      continue
    end
    FS = factor_shape(fR)
    d = lcm([x for (x, v) in FS])
    if d < degree(fR)^2
      break
    end
  end
  return p, d
end

function polredabs(K::AnticNumberField)
  #intended to implement
  # http://beta.lmfdb.org/knowledge/show/nf.polredabs
  #as in pari
  #TODO: figure out the separation of T2-norms....
  ZK = lll(maximal_order(K))
  I = index(ZK)^2
  D = discriminant(ZK)
  B = basis(ZK, copy = false)
  Zx = FlintZZ["x"][1]
  f = Zx(K.pol)
  p, d = _find_prime(f)

  F = FlintFiniteField(p, d, "w", cached = false)[1]
  Ft = PolynomialRing(F, "t", cached = false)[1]
  ap = zero(Ft)
  fit!(ap, degree(K)+1)
  rt = roots(f, F)

  n = degree(K)

  b = _block(B[1].elem_in_nf, rt, ap)
  i = 2
  while length(b) < degree(K)
    bb = _block(B[i].elem_in_nf, rt, ap)
    b = _meet(b, bb)
    i += 1
  end
  i -= 1
#  println("need to use at least the first $i basis elements...")
  pr = 100
  old = precision(BigFloat)
  E = enum_ctx_from_ideal(ideal(ZK, 1), zero_matrix(FlintZZ, 1, 1), prec = pr, TU = BigFloat, TC = BigFloat)
  setprecision(BigFloat, pr)
  while true
    try
      if E.C[end] + 0.0001 == E.C[end]  # very very crude...
        pr *= 2
        continue
      end
      break
    catch e
      if isa(e, InexactError) || isa(e, LowPrecisionLLL) || isa(e, LowPrecisionCholesky)
        pr *= 2
        continue
      end
      rethrow(e)
    end
    setprecision(BigFloat, pr)
    E = enum_ctx_from_ideal(ideal(ZK, 1), zero_matrix(FlintZZ, 1, 1), prec = pr, TU = BigFloat, TC = BigFloat)
  end

  l = zeros(FlintZZ, n)
  l[i] = 1

  scale = 1.0
  enum_ctx_start(E, matrix(FlintZZ, 1, n, l), eps = 1.01)

  a = gen(K)
  all_a = nf_elem[a]
  la = length(a)*BigFloat(E.t_den^2)
  Ec = BigFloat(E.c//E.d)
  eps = BigFloat(E.d)^(1//2)

  found_pe = false
  while !found_pe
    while enum_ctx_next(E)
#      @show E.x
      M = E.x*E.t
      q = elem_from_mat_row(K, M, 1, E.t_den)
      bb = _block(q, rt, ap)
      if length(bb) < n
        continue
      end
      found_pe = true
#  @show    llq = length(q)
#  @show sum(E.C[i,i]*(BigFloat(E.x[1,i]) + E.tail[i])^2 for i=1:E.limit)/BigInt(E.t_den^2)
      lq = Ec - (E.l[1] - E.C[1, 1]*(BigFloat(E.x[1,1]) + E.tail[1])^2) #wrong, but where?
#      @show lq/E.t_den^2

      if lq < la + eps
        if lq > la - eps
          push!(all_a, q)
#          @show "new one", q, minpoly(q), bb
        else
          a = q
          all_a = nf_elem[a]
          if lq/la < 0.8
#            @show "re-init"
            enum_ctx_start(E, E.x, eps = 1.01)  #update upperbound
            Ec = BigFloat(E.c//E.d)
          end
          la = lq
  #        @show Float64(la/E.t_den^2)
        end
      end
    end
    scale *= 2
    enum_ctx_start(E, matrix(FlintZZ, 1, n, l), eps = scale)
    Ec = BigFloat(E.c//E.d)
  end

  setprecision(BigFloat, old)
  all_f = Tuple{nf_elem, fmpq_poly}[(x, minpoly(x)) for x=all_a]
  all_d = fmpq[abs(discriminant(x[2])) for x= all_f]
  m = minimum(all_d)
  L1 = Tuple{nf_elem, fmpq_poly}[]
  for i = 1:length(all_f)
    if all_d[i] == m
      push!(L1, all_f[i])
    end
  end
  L2 = Tuple{nf_elem, fmpq_poly}[minQ(x) for x=L1]
  if length(L2) == 1
    return L2[1]
  end
  L3 = sort(L2, lt = il)

  return L3[1]
end

function Q1Q2(f::PolyElem)
  q1 = parent(f)()
  q2 = parent(f)()
  g = gen(parent(f))
  for j=0:degree(f)
    if isodd(j)
      q2 += coeff(f, j)*g^div(j, 2)
    else
      q1 += coeff(f, j)*g^div(j, 2)
    end
  end
  return q1, q2
end

function minQ(A::Tuple{nf_elem, fmpq_poly})
  a = A[1]
  f = A[2]
  q1, q2 = Q1Q2(f)
  if lead(q1)>0 && lead(q2) > 0
    return (-A[1], f(-gen(parent(f)))*(-1)^degree(f))
  else
    return (A[1], f)
  end
end

function int_cmp(a, b)
  if a==b
    return 0
  end
  if abs(a) == abs(b)
    if a>b
      return 1
    else
      return -1
    end
  end
  return cmp(abs(a), abs(b))
end

function il(F, G)
  f = F[2]
  g = G[2]
  i = degree(f)
  while i>0 && int_cmp(coeff(f, i), coeff(g, i))==0
    i -= 1
  end
  return int_cmp(coeff(f, i), coeff(g, i))<0
end

################################################################################
#
#  LLL for simplify
#
################################################################################

function _lll_for_simplify(M::NfOrd; prec = 100)

  K = nf(M)

  if signature(M)[1] == degree(K)
    #in this case the gram-matrix of the minkowski lattice is the trace-matrix
    #which is exact.
    BM = _lll_gram(ideal(M, 1))[2]
    O1 = NfOrd(K, BM*basis_matrix(M, copy = false))
    O1.ismaximal = M.ismaximal
    if isdefined(M, :index)
      O1.index = M.index
    end
    if isdefined(M, :disc)
      O1.disc = M.disc
    end
    if isdefined(M, :gen_index)
      O1.gen_index = M.gen_index
    end
    M.lllO = O1
    return O1
  end

  if degree(K) == 2 && discriminant(M) < 0
    #in this case the gram-matrix of the minkowski lattice is related to the
    #trace-matrix which is exact.
    #could be extended to CM-fields
    BM = _lll_quad(ideal(M, 1))[2]
    O1 = NfOrd(K, BM*basis_matrix(M, copy = false))
    O1.ismaximal = M.ismaximal
    if isdefined(M, :index)
      O1.index = M.index
    end
    if isdefined(M, :disc)
      O1.disc = M.disc
    end
    if isdefined(M, :gen_index)
      O1.gen_index = M.gen_index
    end
    M.lllO = O1
    return O1
  end

  n = degree(M)
  prec = max(prec, 10*n)
  local g::fmpz_mat
  local d::fmpz_mat
  fl = true
  ctx = Nemo.lll_ctx(0.75, 0.51, :gram)
  while true

    while true
      try
        d = minkowski_gram_mat_scaled(M, prec)
        break
      catch e
        prec = prec*2
      end
    end
    g = identity_matrix(FlintZZ, n)

    prec = div(prec, 2)
    shift!(d, -prec)  #TODO: remove?

    for i=1:n
      fmpz_mat_entry_add_ui!(d, i, i, UInt(nrows(d)))
    end
    @vtime :Simplify 3 ccall((:fmpz_lll, :libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{Nemo.lll_ctx}), d, g, ctx)

    nb = nbits(maximum(abs, g))
    if nb <=  prec
      if nb > div(prec, 2)
        fl = false
      end
      break
    end
    @vprint :Simplify 3 "Still in the loop\n"
    prec *= 4
  end
  ## lattice has lattice disc = order_disc * norm^2
  ## lll needs to yield a basis sth
  ## l[1,1] = |b_i|^2 <= 2^((n-1)/2) disc^(1/n)
  ## and prod(l[i,i]) <= 2^(n(n-1)/2) disc

  disc = abs(discriminant(M))
  di = root(disc, n)+1
  di *= fmpz(2)^(div(n+1,2)) * fmpz(2)^prec
  if cmpindex(d, 1, 1, di) > 0
    fl = false
  else
    pr = prod_diagonal(d)
    if pr > fmpz(2)^(div(n*(n-1), 2)) * disc * fmpz(2)^(n*prec)
      fl = false
    end
  end
  On = NfOrd(K, g*basis_matrix(M, copy = false))
  On.ismaximal = M.ismaximal
  if isdefined(M, :index)
    On.index = M.index
  end
  if isdefined(M, :disc)
    On.disc = M.disc
  end
  if isdefined(M, :gen_index)
    On.gen_index = M.gen_index
  end
  if fl
    M.lllO = On
    return On
  else
    On1 = _lll_for_simplify(On, prec = prec*4)
    M.lllO = On1
    return On1
  end
end


text:
#!/usr/bin/env julia
using Gtk, Gtk.ShortNames

win = Window("Calculate Area", 400, 100)
hbox = Box(:h)
vbox = Box(:v)
label = Label("")
fh = Frame("Height")
fw = Frame("Width")
hent = Entry(); set_gtk_property!(hent, :text, 3)
went = Entry(); set_gtk_property!(went, :text, 4)

push!(win, hbox)
push!(hbox, label)
push!(hbox, vbox)
push!(vbox, fh); push!(fh, hent)
push!(vbox, fw); push!(fw, went)

set_gtk_property!(hbox, :expand, label, true)
set_gtk_property!(label, :use_markup, true)
set_gtk_property!(label, :label, "" * string(*(get_gtk_property(hent, :text, String) |> Meta.parse |> eval, get_gtk_property(went, :text, String) |> Meta.parse |> eval)) * "")
showall(win)

function calarea()
    set_gtk_property!(label, :label, "" * string(*(get_gtk_property(hent, :text, String) |> Meta.parse |> eval, get_gtk_property(went, :text, String) |> Meta.parse |> eval)) * "")

    return nothing
end

signal_connect(x -> calarea(), went, "activate")
signal_connect(x -> calarea(), hent, "activate")

if !isinteractive()
    c = Condition()
    signal_connect(win, :destroy) do widget
        notify(c)
    end
    wait(c)
end

max_stars_repo_name: goropikari/JuliaGUIAppExample

max_issues_repo_name: goropikari/JuliaGUIAppExample

max_forks_repo_name: goropikari/JuliaGUIAppExample


text:
struct HydroDispatchReservoirBudgetUpperBound <: PSI.AbstractHydroReservoirFormulation end

function PSI.DeviceRangeConstraintSpec(
    ::Type{<:PSI.RangeConstraint},
    ::Type{PSI.ActivePowerVariable},
    ::Type{T},
    ::Type{<:HydroDispatchReservoirBudgetUpperBound},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, PSI.AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    if !use_parameters && !use_forecasts
        return PSI.DeviceRangeConstraintSpec(;
            range_constraint_spec = PSI.RangeConstraintSpec(;
                constraint_name = PSI.make_constraint_name(
                    PSI.RangeConstraint,
                    PSI.ActivePowerVariable,
                    T,
                ),
                variable_name = PSI.make_variable_name(PSI.ActivePowerVariable, T),
                limits_func = x -> (min = 0.0, max = PSY.get_active_power(x)),
                constraint_func = PSI.device_range!,
                constraint_struct = PSI.DeviceRangeConstraintInfo,
            ),
        )
    end

    return PSI.DeviceRangeConstraintSpec(;
        timeseries_range_constraint_spec = PSI.TimeSeriesConstraintSpec(
            constraint_name = PSI.make_constraint_name(
                PSI.RangeConstraint,
                PSI.ActivePowerVariable,
                T,
            ),
            variable_name = PSI.make_variable_name(PSI.ActivePowerVariable, T),
            parameter_name = use_parameters ? PSI.ACTIVE_POWER : nothing,
            forecast_label = "max_hourly_hydro_budget",
            multiplier_func = x -> PSY.get_max_active_power(x),
            constraint_func = use_parameters ? PSI.device_timeseries_param_ub! :
                              PSI.device_timeseries_ub!,
        ),
    )
end


text:
using MOCNeutronTransport
using StaticArrays
@testset "Point_3D" begin
    @testset "$T" for T in [Float32, Float64]
        @testset "Constructors" begin
            # 3D single constructor
            p = Point_3D(T(1), T(2), T(3))
            @test p.x == SVector(T.((1, 2, 3)))
            @test typeof(p.x) == typeof(SVector(T.((1, 2, 3))))

            # 2D single constructor
            p = Point_3D(T(1), T(2))
            @test p.x == SVector(T.((1, 2, 0)))
            @test typeof(p.x) == typeof(SVector(T.((1, 2, 0))))

            # 1D single constructor
            p = Point_3D(T(1))
            @test p.x == SVector(T.((1, 0, 0)))
            @test typeof(p.x) == typeof(SVector(T.((1, 0, 0))))

            # 3D tuple constructor
            p = Point_3D( T.((1, 2, 3)) )
            @test p.x == SVector(T.((1, 2, 3)))
            @test typeof(p.x) == typeof(SVector(T.((1, 2, 3))))

            # 2D tuple constructor
            p = Point_3D( T.((1, 2)) )
            @test p.x == SVector(T.((1, 2, 0)))
            @test typeof(p.x) == typeof(SVector(T.((1, 2, 0))))

            # 3D single conversion constructor
            p = Point_3D(Float64, 1, 2, 3)
            @test p.x == SVector(Float64.((1, 2, 3)))
            @test typeof(p.x) == typeof(SVector(Float64.((1, 2, 3))))

            # 2D single conversion constructor
            p = Point_3D(Float64, 1, 2)
            @test p.x == SVector(Float64.((1, 2, 0)))
            @test typeof(p.x) == typeof(SVector(Float64.((1, 2, 0))))

            # 1D single conversion constructor
            p = Point_3D(Float64, 1)
            @test p.x == SVector(Float64.((1, 0, 0)))
            @test typeof(p.x) == typeof(SVector(Float64.((1, 0, 0))))
        end

        @testset "Base" begin
            p = Point_3D(T, 1, 2, 3)

            # zero
            p₀ = zero(p)
            @test all(p₀.x .== T(0))
            @test typeof(p.x[1]) == typeof(T(0))

            # (::Type)
            q = Float64(p)
            @test q.x == SVector(Float64.((1, 2, 3)))
            @test typeof(q.x) == typeof(SVector(Float64.((1, 2, 3))))
        end

        @testset "Operators" begin
            p₁ = Point_3D(T, 1, 2, 0)
            p₂ = Point_3D(T, 2, 4, 6)

            # Point_3D equivalence
            @test p₁ == Point_3D(T, 1, 2, 0)

            # Point_3D isapprox
            p = Point_3D(T, 1, 2 - 10*eps(T), 0 - 10*eps(T))
            @test T(2) ≈ T(2 - 10*eps(T))
            @test isapprox(T(0), T(0 - 10*eps(T)), atol=sqrt(eps(T)))
            @test p ≈ p₁

            # Point_3D addition
            p = p₁ + p₂
            @test p.x == SVector(T.((3,6,6)))
            @test typeof(p.x) == typeof(SVector(T.((3, 6, 6))))

            # Point_3D subtraction
            p = p₁ - p₂
            @test p.x == SVector(T.((-1, -2, -6)))
            @test typeof(p.x) == typeof(SVector(T.((-1, -2, -6))))

            # Cross product
            p₁ = Point_3D(T, 2, 3, 4)
            p₂ = Point_3D(T, 5, 6, 7)
            p = p₁ × p₂
            @test p == Point_3D(T, -3, 6, -3)
            @test typeof(p.x) == typeof(SVector(T.((-3, 6, -3))))

            # Dot product
            @test p₁ ⋅ p₂ == T(10 + 18 + 28)

            # Number addition
            p₁ = Point_3D(T, 1, 2, 3)
            p₂ = Point_3D(T, 2, 4, 6)
            p = p₁ + T(1)
            @test p  == Point_3D(T, 2, 3, 4)
            @test typeof(p.x) == typeof(SVector(T.((2, 3, 4))))

            # Broadcast addition
            parray = [p₁, p₂]
            parray = parray .+ T(1)
            @test parray[1] == p₁ + T(1)
            @test parray[2] == p₂ + T(1)
            @test typeof(parray[1].x) == typeof(SVector(T.((2, 3, 4))))

            # Subtraction
            p = p₁ - T(1)
            @test p == Point_3D(T, 0, 1, 2)
            @test typeof(p.x) == typeof(SVector(T.((0, 1, 2))))

            # Multiplication
            p = T(4)*p₁
            @test p == Point_3D(T, 4, 8, 12)
            @test typeof(p.x) == typeof(SVector(T.((4, 8, 12))))

            # Division
            p = p₁/T(4)
            @test p == Point_3D(T, 1//4, 1//2, 3//4)
            @test typeof(p.x) == typeof(SVector(T.((4, 8, 12))))

            # Unary -
            p = -p₁
            @test p == Point_3D(T, -1, -2, -3)
            @test typeof(p.x) == typeof(SVector(T.((4, 8, 12))))
        end

        @testset "Methods" begin
            p₁ = Point_3D(T, 1, 2, 3)
            p₂ = Point_3D(T, 2, 4, 6)

            # distance
            @test distance(p₁, p₂) == sqrt(T(14))
            @test typeof(distance(p₁, p₂)) == typeof(T(1))

            # norm
            @test norm(p₁) == sqrt(T(14))
            @test norm(p₂) == sqrt(T(4 + 16 + 36))
            @test typeof(norm(p₁)) == typeof(T(1))
        end
    end
end


text:
tests = [
    "transform_stats",
    "wavelet",
    "point_field"
]

println("Running tests:")

for t in tests
    println(" * $t.jl")
    include("$t.jl")
end

# TODO point field tests


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie, ColorSchemes
CairoMakie.activate!() # HIDE
let
    x = LinRange(0,2π,50) # with 50 points in the interval.
    fig = Figure(resolution = (1200, 600), font =:sans, fontsize = 22)
    ax = Axis(fig, xlabel = "x", ylabel = "")
    line1 = lines!(x, sin.(x), color = x,  colormap = :thermal, linewidth = 2)
    line2 = lines!(x, cos.(x), color = sqrt.(x),  colormap = :ice,linewidth = 2)
    line3 = lines!(x, -sin.(x), color = x.^2, colormap = :viridis,linewidth = 2)
    line4 = lines!(x, -cos.(x), color = x/2, colormap = :plasma,linewidth = 2)
    lineas = [line1, line2, line3, line4]
    labels = ["x", "sqrt(x)", "x²", "x/2"]
    cbars = [Colorbar(fig, lineas[i], label = labels[i]) for i in 1:4]
    fig[1, 1] = ax
    fig[1, 2] = cbars[1]
    fig[1, 3] = cbars[2]
    fig[1, 4] = cbars[3]
    fig[1, 5] = cbars[4]
    colgap!(fig.layout, 15)
    fig
    #save("MultipleMulticoloredLineCbarOut.png", fig, px_per_unit = 2)
    save(joinpath(@__DIR__, "output", "MultipleMulticoloredLineCbarOut.png"), fig, px_per_unit = 2) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "ColorSchemes"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using NonstationaryProcesses
using Plots

savedt = 0.01
S = lorenzSim(
    X0 = [0.0, -0.01, 9.0],
    parameter_profile = (constantParameter, constantParameter, constantParameter),
    parameter_profile_parameters = (10.0, 28.0, 8/3), # Sprott's recomendation
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = savedt,
    tmax = 1000.0,
    alg = AutoVern7(Rodas5()),
    solver_opts = Dict(:adaptive => true, :reltol => 1e-12))

P = NonstationaryProcesses.lowpass(S, rampInterval, (1/(2*savedt), 0.0, 200.0, 1000.0), nwindows = 100)()

a = animatespectrum(P, downsample=50, trail=5000, colorgradient=cgrad([:black, :crimson]), phasogram=true, nwindows=100, dpi=100)
gif(a, "./lowpasslorenzsym.mp4", fps=48)


text:
# Triangulations of domains
# ===========================
#
# These examples can be loaded into Julia (Revise.jl recommended)
# and run by calling one of the methods with the optional arguments "Plotter=PyPlot".
# Alternatively, you can download a [jupyter notebook](example_domains.ipynb) created
# from this source.
#---
#
# Set  up environment

using Triangulate
using Test
using Printf

injupyter()=(isdefined(Main, :IJulia) && Main.IJulia.inited)

if injupyter()
    import PyPlot
end

# ### Constrained Delaunay triangulation (CDT) of a domain given by a segment list specifying its boundary.
#
# This is obtained by
# specifying the "p" flag.
function example_domain_cdt(;Plotter=nothing)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([0.0 0.0 ; 1.0 0.0 ; 1.0  1.0 ; 0.6 0.6; 0.0 1.0]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ; 3 4 ; 4 5 ; 5 1 ]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4, 5])
    display(triin)
    (triout, vorout)=triangulate("pQ", triin)
    display(triout)
    plot_in_out(Plotter,triin,triout,title="Domain triangulation")
    @test numberofpoints(triout)>=numberofpoints(triin)
    @test numberofsegments(triout)>=numberofsegments(triin)
    @test numberoftriangles(triout)>0
end
#
injupyter() && example_domain_cdt(Plotter=PyPlot);

# ### Constrained Delaunay triangulation (CDT) of a domain given by a segment list specifying its boundary together with a maximum  area constraint.
#
# This constraint is specfied as a floating
# point number given after the -a flag.
# Be careful to not give it in the exponential format as Triangle would be unable to analyse it.
# Therefore it is dangerous to a number in the string interpolation and it is better to convert
# it to a string before using `@sprintf`.
# Specifying only the maximum area constraint does not prevent very thin
# triangles from occuring at the boundary.
function example_domain_cdt_area(;Plotter=nothing,maxarea=0.05)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([0.0 0.0 ; 1.0 0.0 ; 1.0  1.0 ; 0.6 0.6; 0.0 1.0]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ; 3 4 ; 4 5 ; 5 1 ]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4, 5])
    area=@sprintf("%.15f",maxarea) # Don't use exponential format!
    display(triin)
    (triout, vorout)=triangulate("pa$(area)Q", triin)
    display(triout)
    plot_in_out(Plotter,triin,triout,voronoi=vorout, title="Domain CDT with area constraint")
    @test numberofpoints(triout)>=numberofpoints(triin)
    @test numberofsegments(triout)>=numberofsegments(triin)
    @test numberoftriangles(triout)>0
end
#
injupyter() && example_domain_cdt_area(Plotter=PyPlot,maxarea=0.05);

# ### Boundary conforming  Delaunay triangulation (BCDT) of a domain given by a segment list specifying its boundary
# In addition to the area constraint specify the -D flag
# in order to keep the triangle circumcenters  within the domain.
function example_domain_bcdt_area(;Plotter=nothing,maxarea=0.05)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([0.0 0.0 ; 1.0 0.0 ; 1.0  1.0 ; 0.6 0.6; 0.0 1.0]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ; 3 4 ; 4 5 ; 5 1 ]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4, 5])
    display(triin)
    area=@sprintf("%.15f",maxarea)
    (triout, vorout)=triangulate("pa$(area)DQ", triin)
    display(triout)
    plot_in_out(Plotter,triin,triout,voronoi=vorout, title="Boundary conforming Delaunay triangulation")
    @test numberofpoints(triout)>=numberofpoints(triin)
    @test numberofsegments(triout)>=numberofsegments(triin)
    @test numberoftriangles(triout)>0
end
#
injupyter() && example_domain_bcdt_area(Plotter=PyPlot,maxarea=0.05);

# ### Constrained Delaunay triangulation of a domain with minimum angle condition
#
# The "q" flag  allows to specify a minimum angle
# constraint preventing skinny triangles.
#
# This combination of flags, possibly with an additional "D" flag is recommended
# when creating triangulations for finite element or finite volume methods.
# It the mimimum angle is larger then 28.6 degrees, Triangle's algorithm may
# run into an infinite loop.
function example_domain_qcdt_area(;Plotter=nothing,minangle=20, maxarea=0.05)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([0.0 0.0 ; 1.0 0.0 ; 1.0  1.0 ; 0.6 0.6; 0.0 1.0]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ; 3 4 ; 4 5 ; 5 1 ]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4, 5])
    display(triin)
    area=@sprintf("%.15f",maxarea)
    angle=@sprintf("%.15f",minangle)
    (triout, vorout)=triangulate("pa$(area)q$(angle)", triin)
    display(triout)
    plot_in_out(Plotter,triin,triout,voronoi=vorout, title="Quality triangulation")
    @test numberofpoints(triout)>=numberofpoints(triin)
    @test numberofsegments(triout)>=numberofsegments(triin)
    @test numberoftriangles(triout)>0

end
#
injupyter() && example_domain_qcdt_area(Plotter=PyPlot,maxarea=0.05,minangle=20);

# ### Triangulation of a domain with refinement callback
#
# A maximum area constraint is specified in the `unsuitable` callback
# which is activated via the "u" flag if it has been passed before calling triangulate.
# In addition, the "q" flag  allows to specify a minimum angle
# constraint preventing skinny triangles.
function example_domain_localref(;Plotter=nothing,minangle=20)
    center_x=0.6
    center_y=0.6
    localdist=0.1
    function unsuitable(x1,y1,x2,y2,x3,y3,area)
        bary_x=(x1+x2+x3)/3.0
        bary_y=(y1+y2+y3)/3.0
        dx=bary_x-center_x
        dy=bary_y-center_y
        qdist=dx^2+dy^2
        qdist>1.0e-5 && area>0.1*qdist
    end

    triunsuitable(unsuitable)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([0.0 0.0 ; 1.0 0.0 ; 1.0  1.0 ; 0.6 0.6; 0.0 1.0]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ; 3 4 ; 4 5 ; 5 1 ]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4, 5])
    display(triin)
    angle=@sprintf("%.15f",minangle)
    (triout, vorout)=triangulate("pauq$(angle)Q", triin)
    display(triout)
    plot_in_out(Plotter,triin,triout,voronoi=vorout, title="Quality triangulation with local refinement")
    @test numberofpoints(triout)>=numberofpoints(triin)
    @test numberofsegments(triout)>=numberofsegments(triin)
    @test numberoftriangles(triout)>0
end
#
injupyter() && example_domain_localref(Plotter=PyPlot,minangle=20);

# ### Triangulation of a heterogeneous domain
#
# The segment list specifies its boundary and the inner boundary between subdomains.
# An additional region list is specified which provides "region points" in `regionlist[1,:]`
# and `regionlist[2,:]`.  These kind of mark the subdomains. `regionlist[3,:]` contains an attribute
# which labels the subdomains. `regionlist[4,:]` contains a maximum area value. `size(regionlist,2)`
# is the number of regions.
#
# With the "A" flag, the subdomain labels are spread to all triangles in the corresponding
# subdomains, becoming available in `triangleattributelist[1,:]`.
# With the "a" flag, the area constraints are applied in the corresponding subdomains.
function example_domain_regions(;Plotter=nothing,minangle=20)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([0.0 0.0 ;0.5 0.0; 1.0 0.0 ; 1.0  1.0 ; 0.6 0.6; 0.0 1.0]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ;3 4 ; 4 5 ; 5 6 ; 6 1 ; 2 5]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4, 5, 6, 7])
    triin.regionlist=Matrix{Cdouble}([0.2 0.8; 0.2 0.2; 1 2 ; 0.01 0.05])
    display(triin)
    angle=@sprintf("%.15f",minangle)
    (triout, vorout)=triangulate("paAq$(angle)Q", triin)
    display(triout)

    plot_in_out(Plotter,triin,triout,voronoi=vorout, title="Hetero domain triangulation")

    @test numberofpoints(triout)>=numberofpoints(triin)
    @test numberofsegments(triout)>=numberofsegments(triin)
    @test numberoftriangles(triout)>0
end
#
injupyter() && example_domain_regions(Plotter=PyPlot,minangle=20);

# ### Triangulation of a domain with holes
#
# The segment list specifies its boundary and the boundaries of the holes.
# An additional hole list is specified which provides "hole points" in `holelist[1,:]`
# and `holelist[2,:]`.
function example_domain_holes(;Plotter=nothing,minangle=20,maxarea=0.001)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([0.0 0.0;
                                     1.0 0.0;
                                     1.0 1.0;
                                     0.0 1.0;
                                     0.2 0.2;
                                     0.3 0.2;
                                     0.3 0.3;
                                     0.2 0.3;
                                     0.6 0.6;
                                     0.7 0.6;
                                     0.7 0.7;
                                     0.6 0.7;
                                     ]')
    triin.segmentlist=Matrix{Cint}([1 2; 2 3; 3 4; 4 1; 5 6; 6 7; 7 8; 8 5;  9 10; 10 11; 11 12; 12 9;]')
    triin.segmentmarkerlist=Vector{Int32}([1, 1,1,1, 2,2,2,2, 3,3,3,3])
    triin.holelist=[0.25 0.25; 0.65 0.65;]'
    display(triin)
    area=@sprintf("%.15f",maxarea) # Don't use exponential format!
    angle=@sprintf("%.15f",minangle)
    (triout, vorout)=triangulate("pa$(area)q$(angle)Q", triin)
    display(triout)

    plot_in_out(Plotter,triin,triout,voronoi=vorout, title="Domain with holes")

    @test numberofpoints(triout)>=numberofpoints(triin)
    @test numberofsegments(triout)>=numberofsegments(triin)
    @test numberoftriangles(triout)>0
end
#
injupyter() && example_domain_holes(Plotter=PyPlot,minangle=20,maxarea=0.05);