text:
#   Author: V. Vitale
#   Feb 2022
#   TDVP

include("mps.jl")
include("mpo.jl")
include("contractions.jl")

using TensorOperations
using KrylovKit
using LinearAlgebra

function tdvp!(psi::MPS, W::MPO,dt::Complex,is_hermitian::Bool; krylovdim=10,sweeps=1,tol=1e-15,chimax=128)

    right_orthogonalize!(psi)

    d = dims(psi)[2][2]
    L = construct_L(psi, W)
    R = construct_R(psi, W)

    #global psi,L,R
    Energy=0
    for sweep in 1:Int(sweeps)
        for i in 1:psi.N-1
            psi.data[i],psi.data[i+1] = evolve_right(psi.data[i],psi.data[i+1],W.data[i],W.data[i+1],
                                                     L[i], R[i+1], dt/2, krylovdim,is_hermitian,tol,chimax)
            if i!=psi.N-1
                L[i+1] = contract_from_left(L[i], psi.data[i], W.data[i])
                psi.data[i+1] = local_step(psi.data[i+1], W.data[i+1],
                                           L[i+1], R[i+1],dt/2 ,krylovdim,is_hermitian,tol)
            end
        end
        for i in psi.N:-1:2
            psi.data[i-1],psi.data[i] = evolve_left(psi.data[i-1], psi.data[i], W.data[i-1],  W.data[i],
                                                    L[i-1], R[i], dt/2, krylovdim,is_hermitian,tol,chimax)
            if i!=2
                R[i-1] = contract_from_right(R[i], psi.data[i], W.data[i])
                psi.data[i-1] = local_step(psi.data[i-1], W.data[i-1],
                                           L[i-1], R[i-1],dt/2,krylovdim,is_hermitian,tol)
            end
        end
    end
    right_orthogonalize!(psi)
    return psi
end

function evolve_right(AL::Array, AR::Array, WL::Array, WR::Array, E::Array, F::Array, dt::Complex, krylovdim::Int,is_hermitian::Bool,tol::Float64,chimax::Int)

    sAL = size(AL)
    sAR = size(AR)
    @tensor A[:] := AL[-1,-2,1]*AR[1,-3,-4]
    A = reshape(A,(sAL[1],sAL[2]*sAR[2],sAR[3]))
    sWL = size(WL)
    sWR = size(WR)

    @tensor M[:] := WL[-1,1,-3,-5]*WR[1,-2,-4,-6]
    M = reshape(M,(sWL[1],sWR[2],sWL[3]*sWR[3],sWL[4]*sWR[4]))
    function H_lin(v)
        @tensor temp_1[ :] := E[ 1, -2, -3] *  v[ 1, -1, -4]
        @tensor temp_2[:] :=  temp_1[-1,-2,-3,1] * F[1,-4,-5]
        @tensor temp[:] := temp_2[1, 2, -1, 3,-3] * M[ 2 3 ; 1 -2 ]
        return temp
    end

    #en,V,info = eigsolve( H_lin , A ,  1  , :SR ; issymmetric = true  )
    V,info = exponentiate( H_lin , -dt ,  A ; ishermitian = is_hermitian, tol=tol ,krylovdim=krylovdim)

    V=reshape(V,(sAL[1]*sAL[2],sAR[2]*sAR[3]))
    U,S,V = svd(V,full=false)
    V=V'
    #S=S/norm(S)
    indices = findall(norm(S) .-cumsum(S.^2) .< tol)
    if length(indices)>0
        chi = indices[1]+1
    else
        chi = size(S)[1]
    end
    if chi>chimax
        chi=chimax
    end
    if size(S)[1] > chi
        U = U[:,1:chi]
        S = S[1:chi]
        V =  V[1:chi,:]
    end

    #S /= norm(S)

    AL = reshape( U , ( sAL[1], sAL[2], :) )
    #"ij,jl,slk->sik"
    S=diagm(S)
    @tensor AR[:] := S[-1,1 ] * V[ 1,-2 ]
    AR=reshape(AR,(:,sAR[2],sAR[3]))
    return AL, AR
end

function evolve_left(AL::Array, AR::Array, WL::Array, WR::Array, E::Array, F::Array, dt::Complex, krylovdim::Int,is_hermitian::Bool,tol::Float64,chimax::Int)

    sAL = size(AL)
    sAR = size(AR)

    @tensor A[:] :=AL[-1,-2,1]*AR[1,-3,-4]
    A = reshape(A,(sAL[1],sAL[2]*sAR[2],sAR[3]))

    sWL = size(WL)
    sWR = size(WR)

    @tensor M[:] := WL[-1,1,-3,-5]*WR[1,-2,-4,-6]
    M = reshape(M,(sWL[1],sWR[2],sWL[3]*sWR[3],sWL[4]*sWR[4]))

    function H_lin(v)
            @tensor temp_1[ :] := E[ 1, -2, -3] *  v[ 1, -1, -4]
            @tensor temp_2[:] :=  temp_1[-1,-2,-3,1] * F[1,-4,-5]
            @tensor temp[:] := temp_2[1, 2, -1, 3,-3] * M[ 2 3 ; 1 -2 ]
        return temp
    end

    V,info = exponentiate( H_lin , -dt ,  A ; ishermitian = is_hermitian, tol=tol ,krylovdim=krylovdim)

    sV=size(V)
    V=reshape(V,(sAL[1]*sAL[2],sAR[2]*sAR[3]))

    U,S,V = svd(V,full=false)
    V=V'
    #S /=norm(S)
    indices = findall(norm(S) .-cumsum(S.^2) .< tol)
    if length(indices)>0
        chi = indices[1]+1
    else
        chi = size(S)[1]
    end
    if chi>chimax
        chi=chimax
    end
    if size(S)[1] > chi
        U = U[:,1:chi]
        S = S[1:chi]
        V =  V[1:chi,:]
    end

    #S /= norm(S)

    AR = reshape(V,(:,sAR[2],sAR[3]))

    S=diagm(S)
    @tensor AL[:] :=  U[ -1,1 ] * S[ 1,-2 ]
    AL = reshape( AL , ( sAL[1], sAL[2], :) )
    return AL, AR
end

function local_step(A::Array, M::Array, E::Array, F::Array, dt::Complex, krylovdim::Int,is_hermitian::Bool,tol::Float64)

    sA = size(A)

        function H_lin(v)
        @tensor temp_1[ :] := E[ 1, -2, -3] *  v[ 1, -1, -4]
        @tensor temp_2[:] :=  temp_1[-1,-2,-3,1] * F[1,-4,-5]
        @tensor temp[:] := temp_2[1, 2, -1, 3,-3] * M[ 2 3 ; 1 -2 ]
        return temp
    end

    V,info = exponentiate( H_lin , dt ,  A ; ishermitian = is_hermitian, tol=tol ,krylovdim=krylovdim)

    A=reshape(V,(sA[1],sA[2],sA[3]))
    return A
end

max_stars_repo_name: aanniikkeettr/Tensor-Networks

max_issues_repo_name: aanniikkeettr/Tensor-Networks

max_forks_repo_name: aanniikkeettr/Tensor-Networks


text:
export ExponentialLinearQuadratic

import Distributions: pdf, logpdf, ContinuousUnivariateDistribution
import Base: eltype
struct ExponentialLinearQuadratic{ A <: AbstractApproximationMethod, T <: Real } <: ContinuousUnivariateDistribution
    approximation :: A
    a :: T
    b :: T
    c :: T
    d :: T
end

Base.eltype(::Type{ <: ExponentialLinearQuadratic{A, T} }) where { A, T } = T
Base.eltype(::ExponentialLinearQuadratic{A, T})            where { A, T } = T

ExponentialLinearQuadratic(approximation, a::Real, b::Real, c::Real, d::Real)             = ExponentialLinearQuadratic(approximation, promote(a, b, c, d)...)
ExponentialLinearQuadratic(approximation, a::Integer, b::Integer, c::Integer, d::Integer) = ExponentialLinearQuadratic(approximation, float(a), float(b), float(c), float(d))

mean_cov(dist::ExponentialLinearQuadratic) = mean_var(dist)

function mean_var(dist::ExponentialLinearQuadratic)
    # This is equivalent to (x) -> pdf(dist, x) * exp(x^2 / 2)
    adjusted_pdf = let a = dist.a, b = dist.b, c = dist.c, d = dist.d
        (x) -> exp(-(a * x - x^2 + b * exp(c * x + d * x^2 / 2)) / 2)
    end
    return approximate_meancov(dist.approximation, adjusted_pdf, NormalMeanVariance(zero(eltype(dist)), one(eltype(dist))))
end

mean_invcov(dist::ExponentialLinearQuadratic)      = mean_cov(dist) .|> (identity, inv)
mean_std(dist::ExponentialLinearQuadratic)         = mean_var(dist) .|> (identity, sqrt)
weightedmean_cov(dist::ExponentialLinearQuadratic) = weightedmean_var(dist)
weightedmean_std(dist::ExponentialLinearQuadratic) = weightedmean_var(dist) .|> (identity, sqrt)

function weightedmean_var(dist::ExponentialLinearQuadratic)
    m, v = mean_cov(dist)
    return (inv(v) * m, v)
end

function weightedmean_invcov(dist::ExponentialLinearQuadratic)
    m, w = mean_invcov(dist)
    return (w * m, w)
end

Distributions.pdf(dist::ExponentialLinearQuadratic, x::Real)    = exp(logpdf(dist, x))
Distributions.logpdf(dist::ExponentialLinearQuadratic, x::Real) = -(dist.a * x + dist.b * exp(dist.c * x + dist.d * x ^ 2 / 2)) / 2
Distributions.mean(dist::ExponentialLinearQuadratic)            = mean_var(dist)[1]
Distributions.var(dist::ExponentialLinearQuadratic)             = mean_var(dist)[2]
Distributions.std(dist::ExponentialLinearQuadratic)             = mean_std(dist)[2]
Distributions.cov(dist::ExponentialLinearQuadratic)             = var(dist)

invcov(dist::ExponentialLinearQuadratic)       = mean_invcov(dist)[2]
precision(dist::ExponentialLinearQuadratic)    = mean_invcov(dist)[2]
weightedmean(dist::ExponentialLinearQuadratic) = weightedmean_invcov(dist)[1]

prod_analytical_rule(::Type{ <: UnivariateNormalDistributionsFamily }, ::Type{ <: ExponentialLinearQuadratic }) = ProdAnalyticalRuleAvailable()
prod_analytical_rule(::Type{ <: ExponentialLinearQuadratic }, ::Type{ <: UnivariateNormalDistributionsFamily }) = ProdAnalyticalRuleAvailable()

function prod(::ProdAnalytical, left::UnivariateNormalDistributionsFamily, right::ExponentialLinearQuadratic)
    mean, variance = approximate_meancov(right.approximation, (z) -> pdf(right, z), left)
    return NormalMeanVariance(mean, variance)
end

function prod(::ProdAnalytical, left::ExponentialLinearQuadratic, right::UnivariateNormalDistributionsFamily)
    mean, variance = approximate_meancov(left.approximation, (z) -> pdf(left, z), right)
    return NormalMeanVariance(mean, variance)
end

ext: jl


text:
import GLPK

@testset "DERs" begin
    include("fixedLoad.jl")
    include("deferrableLoad.jl")
    include("curtailableLoad.jl")
    include("shiftableLoad.jl")
    include("battery.jl")
    include("thermalLoad.jl")
end


text:
using Flux

m = Flux.Chain(
    Dense(500,32,relu),
    Dense(32,1,relu)
)

loss(x,y) = Flux.mse(m(x),y)

ps = Flux.param(m)

x = rand(500)
y = rand(50)

opt = ADAM()

dat = [(x,y)]

Flux.train!(loss,ps,dat,opt)

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
const PInv{Style, Typ} = Applied{Style, typeof(pinv), <:Tuple{Typ}}
const Inv{Style, Typ} = Applied{Style, typeof(inv), <:Tuple{Typ}}

Inv(A) = applied(inv, A)
PInv(A) = applied(pinv, A)

const InvOrPInv = Union{PInv, Inv}

parent(A::InvOrPInv) = first(A.args)

pinv(A::PInv) = parent(A)
function inv(A::PInv)
    checksquare(parent(A))
    parent(A)
end

inv(A::Inv) = parent(A)
pinv(A::Inv) = inv(A)

ndims(A::InvOrPInv) = ndims(parent(A))

size(A::InvOrPInv) = reverse(size(parent(A)))
axes(A::InvOrPInv) = reverse(axes(parent(A)))
size(A::InvOrPInv, k) = size(A)[k]
axes(A::InvOrPInv, k) = axes(A)[k]
eltype(A::InvOrPInv) = Base.promote_op(inv, eltype(parent(A)))

# Use ArrayLayouts.ldiv instead of \
struct LdivStyle <: ApplyStyle end

ApplyStyle(::typeof(\), ::Type{A}, ::Type{B}) where {A<:AbstractArray,B<:AbstractArray} = LdivStyle()

axes(M::Applied{Style,typeof(\)}) where Style = ldivaxes(M.args...)
axes(M::Applied{Style,typeof(\)}, p::Int)  where Style = axes(M)[p]
size(M::Applied{Style,typeof(\)}) where Style = length.(axes(M))

@inline eltype(M::Applied{Style,typeof(\)}) where Style = eltype(Ldiv(M.args...))
@inline ndims(M::Applied{Style,typeof(\)}) where Style = ndims(last(M.args))

check_applied_axes(A::Applied{<:Any,typeof(\)}) = check_ldiv_axes(A.args...)

######
# PInv/Inv
########

const PInvMatrix{T,Arg} = ApplyMatrix{T,typeof(pinv),<:Tuple{Arg}}
const InvMatrix{T,Arg} = ApplyMatrix{T,typeof(inv),<:Tuple{Arg}}

PInvMatrix(A) = ApplyMatrix(pinv, A)
function InvMatrix(A)
    checksquare(A)
    ApplyMatrix(inv, A)
end

parent(A::PInvMatrix) = first(A.args)
parent(A::InvMatrix) = first(A.args)
axes(A::PInvMatrix) = reverse(axes(parent(A)))
size(A::PInvMatrix) = map(length, axes(A))
inv(A::InvMatrix) = parent(A)
pinv(A::InvMatrix) = parent(A)
pinv(A::PInvMatrix) = parent(A)

abstract type AbstractInvLayout{L} <: MemoryLayout end
struct InvLayout{L} <: AbstractInvLayout{L} end
struct PInvLayout{L} <: AbstractInvLayout{L} end

applylayout(::Type{typeof(inv)}, ::A) where A = InvLayout{A}()
applylayout(::Type{typeof(pinv)}, ::A) where A = PInvLayout{A}()

# Can always  simplify by lowering to \
simplifiable(::Mul{<:AbstractInvLayout}) = Val(true)

copy(M::Mul{<:AbstractInvLayout}) = ArrayLayouts.ldiv(pinv(M.A), M.B)
copy(M::Mul{<:AbstractInvLayout,<:AbstractLazyLayout}) = ArrayLayouts.ldiv(pinv(M.A), M.B)
@inline copy(M::Mul{<:AbstractInvLayout,ApplyLayout{typeof(*)}}) = simplify(M)
copy(L::Ldiv{<:AbstractInvLayout}) = pinv(L.A) * L.B
Ldiv(A::Applied{<:Any,typeof(\)}) = Ldiv(A.args...)

similar(M::Applied{LdivStyle}, ::Type{T}) where T = similar(Ldiv(M), T)
@inline copy(M::Applied{LdivStyle}) = ldiv(arguments(M)...)
@inline copyto!(dest::AbstractArray, M::Applied{LdivStyle}) = copyto!(dest, Ldiv(M))
@inline materialize!(M::Applied{LdivStyle}) = materialize!(Ldiv(M))

@propagate_inbounds getindex(A::Applied{<:Any,typeof(\)}, kj...) = Ldiv(A)[kj...]

###
# * layout
###
@inline _copy_ldiv_mul(A, B₀, B₁...) = apply(*, A \  B₀,  B₁...)
@inline copy(L::Ldiv{<:DiagonalLayout,ApplyLayout{typeof(*)}}) = _copy_ldiv_mul(L.A, arguments(ApplyLayout{typeof(*)}(), L.B)...)
@inline copy(L::Ldiv{<:Any,ApplyLayout{typeof(*)}}) = _copy_ldiv_mul(L.A, arguments(ApplyLayout{typeof(*)}(), L.B)...)
@inline copy(L::Ldiv{<:AbstractLazyLayout,ApplyLayout{typeof(*)}}) = _copy_ldiv_mul(L.A, arguments(ApplyLayout{typeof(*)}(), L.B)...)
@inline copy(L::Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}) = _copy_ldiv_mul(L.A, arguments(ApplyLayout{typeof(*)}(), L.B)...)

@inline _copy_ldiv_ldiv(B, A₀) = A₀ \ B
@inline _copy_ldiv_ldiv(B, A₀, A₁...) = _copy_ldiv_ldiv(A₀ \ B, A₁...)
@inline copy(L::Ldiv{<:ApplyLayout{typeof(*)}}) = _copy_ldiv_ldiv(L.B, arguments(ApplyLayout{typeof(*)}(), L.A)...)
@inline copy(L::Ldiv{<:ApplyLayout{typeof(*)},<:AbstractLazyLayout}) = _copy_ldiv_ldiv(L.B, arguments(ApplyLayout{typeof(*)}(), L.A)...)
@inline copy(L::Ldiv{<:AbstractLazyLayout,<:AbstractLazyLayout}) = lazymaterialize(\, L.A, L.B)
@inline copy(L::Ldiv{<:AbstractLazyLayout}) = lazymaterialize(\, L.A, L.B)
@inline copy(L::Ldiv{<:Any,<:AbstractLazyLayout}) = lazymaterialize(\, L.A, L.B)

function copy(M::Mul{ApplyLayout{typeof(\)}})
    A,B = arguments(ApplyLayout{typeof(\)}(), M.A)
    A \ (B * M.B)
end
###
# Diagonal
###

inv(D::Diagonal{T,<:LazyVector}) where T = Diagonal(inv.(D.diag))

###
# getindex
###

# \ is likely to be specialised
@propagate_inbounds getindex(Ai::InvMatrix{T}, ::Colon, j::Integer) where T = parent(Ai) \ [Zeros{T}(j-1); one(T); Zeros{T}(size(parent(Ai),1)-j)]
@propagate_inbounds getindex(Ai::PInvMatrix{T}, ::Colon, j::Integer) where T = parent(Ai) \ [Zeros{T}(j-1); one(T); Zeros{T}(size(parent(Ai),1)-j)]
getindex(Ai::SubArray{<:Any,2,<:InvMatrix}, ::Colon, j::Integer) = parent(Ai)[:, parentindices(Ai)[2][j]]

@propagate_inbounds getindex(A::PInvMatrix{T}, k::Integer, j::Integer) where T = A[:,j][k]
@propagate_inbounds getindex(A::InvMatrix{T}, k::Integer, j::Integer) where T = A[:,j][k]

getindex(L::ApplyMatrix{<:Any,typeof(\)}, ::Colon, j::Integer) where T = L.args[1] \ L.args[2][:,j]
getindex(L::ApplyMatrix{<:Any,typeof(\)}, k::Integer, j::Integer) where T = L[:,j][k]


text:
module Univariate

#パッケージの読み込み
using  DataFrames , StatsBase #, Gadfly

export NumericSummary,SummaryMerge,Univariate

function NumericSummary(INPUT::DataArray,VarNames::Symbol)  #数値型の基本統計量の算出
  SummaryVar=var(INPUT)
  SummaryStd=std(INPUT)
  SummaryQuartile=quantile( INPUT , [0.00, 0.25, 0.50, 0.75, 1.00])
  SummaryMedian=median(INPUT)
  SummaryMean=mean(INPUT)

  Output=DataFrame(colname=VarNames,Var=SummaryVar,Std=SummaryStd,Mean=SummaryMean, Min=SummaryQuartile[1] ,Quartile1st=SummaryQuartile[2] ,Median=SummaryMedian,Quartile3rd=SummaryQuartile[4],Max=SummaryQuartile[5])

  return Output
end

function NonNumericSummary(INPUT::DataFrame,VarNames::Symbol)  #非数値型の基本統計量の算出
  Output=by(INPUT, VarNames, df -> size(df, 1))

  return Output
end

function SummaryMerge(INPUT::DataFrame, Summary::DataFrame ,FirstMergeFlag::Bool)  #データの集約
  if(FirstMergeFlag ==true)
  {
    Summary=INPUT
    FirstMergeFlag=false
  } else {
    Summary=vcat(Summary,INPUT)
  } end

  return Dict( {"Summary","FirstMergeFlag"},{Summary,FirstMergeFlag})
end

function Univariate(INPUT)
  VarNames=names(INPUT)
  Vartype=eltypes(INPUT)
  SummaryNonNum=DataFrame(colname="",hist=0)
  SummaryNum=DataFrame(colname="",Var=0.0,Std=0.0,Mean=0.0, Min=0.0 ,Quartile1st=0.0 ,Median=0.0,Quartile3rd=0.0,Max=0.0)
  FirstMergeFlagNum=true
  FirstMergeFlagNonNum=true

  for  i = 1 : length(INPUT)
  {
    if( ( Vartype[i]==(Int) || Vartype[i]==(Float32)  || Vartype[i]== (Float64) ) == true )
    {
      #数値型の基本統計量の算出
      work=NumericSummary(INPUT[:,i],VarNames[i])

      #基本統計量の集約
      SummaryNum=SummaryMerge(work,SummaryNum,FirstMergeFlagNum)["Summary"]
      FirstMergeFlagNum=SummaryMerge(work,SummaryNum,FirstMergeFlagNum)["FirstMergeFlag"]
    } else  {
      #非数値型の基本統計量の算出
      work=NonNumericSummary(INPUT,VarNames[i])

      #基本統計量の集約
      SummaryNonNum=SummaryMerge(work,SummaryNonNum,FirstMergeFlagNonNum)["Summary"]
      FirstMergeFlagNonNum=SummaryMerge(work,SummaryNonNum,FirstMergeFlagNonNum)["FirstMergeFlag"]
    } end
  } end

  return Dict({"SummaryNum","SummaryNonNum"},{SummaryNum ,SummaryNonNum})
end

end # module

max_stars_repo_name: wakakusa/Univariate

max_issues_repo_name: wakakusa/Univariate

max_forks_repo_name: wakakusa/Univariate


text:
# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

using Documenter, Polyhedra

makedocs(
    format = :html,
    sitename = "Polyhedra",
    pages = [
        "Index" => "index.md",
        "Installation" => "installation.md",
        "Representation" => "representation.md",
        "Polyhedron" => "polyhedron.md",
        "Containment/Redundancy" => "redundancy.md",
        "Projection/Elimination" => "projection.md",
        "Optimization" => "optimization.md",
        "Utilities" => "utilities.md"
    ],
    # The following ensures that we only include the docstrings from
    # this module for functions define in Base that we overwrite.
    modules = [Polyhedra]
)

deploydocs(
    repo   = "github.com/JuliaPolyhedra/Polyhedra.jl.git",
    target = "build",
    osname = "linux",
    julia  = "1.0",
    deps   = nothing,
    make   = nothing
)


text:
# Legendre Polynomials

@testset "PlmBar and PlBar" begin
    p = PlmBar(4, 0)
    p′ = zeros(15)
    PlmBar!(p′, 4, 0)
    @test p′ == p

    p, dp = PlmBar_d1(4, 0)
    dp′ = zeros(15)
    PlmBar_d1!(p′, dp′, 4, 0)
    @test p′ == p
    @test dp′ == dp

    p = PlBar(4, 0)
    p′ = zeros(5)
    PlBar!(p′, 4, 0)
    @test p′ == p

    p, dp = PlBar_d1(4, 0)
    dp′ = zeros(5)
    PlBar_d1!(p′, dp′, 4, 0)
    @test p′ == p
    @test dp′ == dp
end

@testset "PlmON and PlON" begin
    p = PlmON(4, 0)
    p′ = zeros(15)
    PlmON!(p′, 4, 0)
    @test p′ == p

    p, dp = PlmON_d1(4, 0)
    dp′ = zeros(15)
    PlmON_d1!(p′, dp′, 4, 0)
    @test p′ == p
    @test dp′ == dp

    p = PlON(4, 0)
    p′ = zeros(5)
    PlON!(p′, 4, 0)
    @test p′ == p

    p, dp = PlON_d1(4, 0)
    dp′ = zeros(5)
    PlON_d1!(p′, dp′, 4, 0)
    @test p′ == p
    @test dp′ == dp
end

@testset "PlmSchmidt and PlSchmidt" begin
    p = PlmSchmidt(4, 0)
    p′ = zeros(15)
    PlmSchmidt!(p′, 4, 0)
    @test p′ == p

    p, dp = PlmSchmidt_d1(4, 0)
    dp′ = zeros(15)
    PlmSchmidt_d1!(p′, dp′, 4, 0)
    @test p′ == p
    @test dp′ == dp

    p = PlSchmidt(4, 0)
    p′ = zeros(5)
    PlSchmidt!(p′, 4, 0)
    @test p′ == p

    p, dp = PlSchmidt_d1(4, 0)
    dp′ = zeros(5)
    PlSchmidt_d1!(p′, dp′, 4, 0)
    @test p′ == p
    @test dp′ == dp
end

@testset "PLegendreA and PLegendre" begin
    p = PLegendreA(4, 0)
    p′ = zeros(15)
    PLegendreA!(p′, 4, 0)
    @test p′ == p

    p, dp = PLegendreA_d1(4, 0)
    dp′ = zeros(15)
    PLegendreA_d1!(p′, dp′, 4, 0)
    @test p′ == p
    @test dp′ == dp

    p = PLegendre(4, 0)
    p′ = zeros(5)
    PLegendre!(p′, 4, 0)
    @test p′ == p

    p, dp = PLegendre_d1(4, 0)
    dp′ = zeros(5)
    PLegendre_d1!(p′, dp′, 4, 0)
    @test p′ == p
    @test dp′ == dp
end

@testset "PlmIndex" begin
    index = PlmIndex(4, 3)
    @test index == 14
end


text:
using LabConnections.BeagleBone
import LabConnections.BeagleBone: getdev, write!, closedev, read, initdev, printdev, listdev

using Base.Test

@testset "SYS LED Tests" begin

    @testset "Inialization/Termination" begin
        # Initialize three devices
        initdev("sysled", Int32(1))
        @test sum(listdev()) == 1

        initdev("sysled", Int32(3))
        @test sum(listdev()) == 2

        # Attempt to initialize a device which has already been initialized
        @test_throws ErrorException initdev("sysled", Int32(1))
        @test_throws ErrorException initdev("sysled", Int32(3))

        # Attempt to initialize a device with a very high index (no matching channel)
        @test_throws ErrorException initdev("sysled", Int32(1000))

        # Attempt to remove devices which have not been initialized
        @test_throws ErrorException closedev("sysled", Int32(2))
        @test_throws ErrorException closedev("sysled", Int32(4))

        # Remove devices from TOC
        closedev("sysled", Int32(1))
        @test sum(listdev()) == 1

        closedev("sysled", Int32(3))
        @test sum(listdev()) == 0
    end

    @testset "Error Handling" begin

        device = initdev("sysled", Int32(1))

        # Test that an exception is thrown when a faulty ID is given
        @test_throws ErrorException write!(device, "bad_entry")

        # Close device
        closedev("sysled", Int32(1))
    end

    @testset "IO Communication" begin
        # Instanciate all possible leds and perform 10 read/write commands
        device1 = initdev("sysled", Int32(1))
        device2 = initdev("sysled", Int32(2))
        device3 = initdev("sysled", Int32(3))
        device4 = initdev("sysled", Int32(4))

        for i = 1:10
            stateA = "$(i%2)"
            stateB = "$((i+1)%2)"
            write!(device1, stateA)
            @test read(device1) == stateA
            write!(device2, stateB)
            @test read(device2) == stateB
            write!(device3, stateA)
            @test read(device3) == stateA
            write!(device4, stateB)
            @test read(device4) == stateB
            sleep(0.1)
        end

        closedev("sysled", Int32(1))
        closedev("sysled", Int32(2))
        closedev("sysled", Int32(3))
        closedev("sysled", Int32(4))
    end
end


text:
using Documenter, NIPALS_PCA

makedocs(
    modules = [NIPALS_PCA],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Fredrik Pettersson",
    sitename = "NIPALS_PCA.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/Fredrikp-ume/NIPALS_PCA.jl.git",
    push_preview = true
)


text:
# ccall wrapper for calling functions in libraries that might not be available

export @checked, @workspace, @argout, @debug_ccall

"""
    @checked function foo(...)
        rv = ...
        return rv
    end

Macro for wrapping a function definition returning a status code. Two versions of the
function will be generated: `foo`, with the function body wrapped by an invocation of the
`@check` macro (to be implemented by the caller of this macro), and `unsafe_foo` where no
such invocation is present and the status code is returned to the caller.
"""
macro checked(ex)
    # parse the function definition
    @assert Meta.isexpr(ex, :function)
    sig = ex.args[1]
    @assert Meta.isexpr(sig, :call)
    body = ex.args[2]
    @assert Meta.isexpr(body, :block)

    # generate a "safe" version that performs a check
    safe_body = quote
        @check $body
    end
    safe_sig = Expr(:call, sig.args[1], sig.args[2:end]...)
    safe_def = Expr(:function, safe_sig, safe_body)

    # generate a "unsafe" version that returns the error code instead
    unsafe_sig = Expr(:call, Symbol("unsafe_", sig.args[1]), sig.args[2:end]...)
    unsafe_def = Expr(:function, unsafe_sig, body)

    return esc(:($safe_def, $unsafe_def))
end

"""
    @argout call(..., out(output_arg), ...)

Change the behavior of a function call, returning the value of an argument instead.
A common use case is to pass a newly-created Ref and immediately dereference that output:

    @argout(some_getter(Ref{Int}()))[]

If no output argument is specified, `nothing` will be returned.
Multiple output arguments return a tuple.

"""
macro argout(ex)
    Meta.isexpr(ex, :call) || throw(ArgumentError("@argout macro should be applied to a function call"))

    block = quote end

    # look for output arguments (`out(...)`)
    output_vars = []
    args = ex.args[2:end]
    for (i,arg) in enumerate(args)
        if Meta.isexpr(arg, :call) && arg.args[1] == :out
            # allocate a variable
            @gensym output_val
            push!(block.args, :($output_val = $(ex.args[i+1].args[2]))) # strip `output(...)`
            push!(output_vars, output_val)

            # replace the argument
            ex.args[i+1] = output_val
        end
    end

    # generate a return
    push!(block.args, ex)
    if isempty(output_vars)
        push!(block.args, :(nothing))
    elseif length(output_vars) == 1
        push!(block.args, :($(output_vars[1])))
    else
        push!(block.args, :(tuple($(output_vars...))))
    end

    esc(block)
end

"""
    @workspace size=getWorkspaceSize(args...) [eltyp=UInt8] [fallback=nothing] buffer->begin
      useWorkspace(workspace, sizeof(workspace))
    end

Create a GPU workspace vector with element type `eltyp` and size in number of elements (in
the default case of an UInt8 element type this equals to the amount of bytes) determined by
calling `getWorkspaceSize`, and pass it to the closure for use in calling `useWorkspace`.
A fallback workspace size `fallback` can be specified if the regular size would lead to OOM.
Afterwards, the buffer is put back into the memory pool for reuse.

This helper protects against the rare but real issue of `getWorkspaceSize` returning
different results based on the GPU device memory pressure, which might change _after_
initial allocation of the workspace (which can cause a GC collection).
"""
macro workspace(ex...)
    code = ex[end]
    kwargs = ex[1:end-1]

    sz = nothing
    eltyp = :UInt8
    fallback = nothing
    for kwarg in kwargs
        key,val = kwarg.args
        if key == :size
            sz = val
        elseif key == :eltyp
            eltyp = val
        elseif key == :fallback
            fallback = val
        else
            throw(ArgumentError("Unsupported keyword argument '$key'"))
        end
    end

    if sz === nothing
        throw(ArgumentError("@workspace macro needs a size argument"))
    end

    # break down the closure to a let block to prevent JuliaLang/julia#15276
    Meta.isexpr(code, :(->)) || throw(ArgumentError("@workspace macro should be applied to a closure"))
    length(code.args) == 2 || throw(ArgumentError("@workspace closure should take exactly one argument"))
    code_arg = code.args[1]
    code = code.args[2]

    return quote
        sz = $(esc(sz))
        workspace = nothing
        try
          while workspace === nothing || length(workspace) < sz
              workspace = CuArray{$(esc(eltyp))}(undef, sz)
              sz = $(esc(sz))
          end
        catch ex
            $fallback === nothing && rethrow()
            isa(ex, OutOfGPUMemoryError) || rethrow()
            workspace = CuArray{UInt8}(undef, $fallback)
        end

        let $(esc(code_arg)) = workspace
            ret = $(esc(code))
            CUDA.unsafe_free!(workspace)
            ret
        end
    end
end

macro debug_ccall(target, rettyp, argtyps, args...)
    @assert Meta.isexpr(target, :tuple)
    f, lib = target.args

    quote
        # get the call target, as e.g. libcuda() triggers initialization, even though we
        # can't use the result in the ccall expression below as it's supposed to be constant
        $(esc(target))

        print($f, '(')
        for (i, arg) in enumerate(($(map(esc, args)...),))
            i > 1 && print(", ")
            render_arg(stdout, arg)
        end
        print(')')
        rv = ccall($(esc(target)), $(esc(rettyp)), $(esc(argtyps)), $(map(esc, args)...))
        println(" = ", rv)
        for (i, arg) in enumerate(($(map(esc, args)...),))
            if arg isa Base.RefValue
                println(" $i: ", arg[])
            end
        end
        rv
    end
end

render_arg(io, arg) = print(io, arg)
render_arg(io, arg::Union{<:Base.RefValue, AbstractArray}) = summary(io, arg)


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

#VERSION >= v"0.4.0-dev" && __precompile__()

module CompilerTools

# package code goes here
include("debug.jl")
include("helper.jl")
include("ast_walk.jl")
include("lambda.jl")
include("read-write-set.jl")
include("CFGs.jl")
include("liveness.jl")
include("OptFramework.jl")
include("udchains.jl")
include("loops.jl")
include("alias-analysis.jl")
include("tdependence.jl")

importall .OptFramework
export @acc, @noacc

end # module


text:
@plugin struct StaticFiles <: Plugin
    files::Array{String} = []
    dests::Array{String} = []
end

# validate() = ...

function hook(p::StaticFiles, t::Template, pkg_dir::AbstractString)
    pkg = basename(pkg_dir)

    for (file, dest) in zip(p.files, p.dests)
        dest = joinpath(pkg_dir, dest)
        mkpath(dirname(dest))
        cp(file, dest)
    end
end


text:
import Random: rand, seed!
import RandomNumbers: AbstractRNG
using AbstractAlgebra #for skipping outside of skipMat

#####
# Mixmax Random Number Generators from "Spectral Test of the MIXMAX Random Number Generators" (2018) by Narek Martirosyan, Konstantin Savvidy and George Savvidy.
# In particular, this implements the N=17 and N=240 versions.
#
# Skipping ahead is possible, and relatively fast (~140 rand calls with nSkip=1000000) for N=17, thanks to a cache for nSkip < 16*2^30.
# The intent of skipping here is not to retrieve far away independent substreams, but to split a stream in blocks.
#
const M61 = UInt64(2305843009213693951)
const BITS = 61

mutable struct Mixmax{N,SPECIALMUL,SPECIAL} <: AbstractRNG{UInt64}
    v::Vector{UInt64}
    counter::Int
    sumtot::UInt64
    #skipMat::Array{UInt64}(undef, (128,17))
end

export Mixmax, Mixmax17, Mixmax240

function Mixmax{N,SPECIALMUL,SPECIAL}(seed::UInt64 = UInt64(123)) where {N,SPECIALMUL,SPECIAL}
    r = Mixmax{N,SPECIALMUL,SPECIAL}(Vector{UInt64}(undef, N), 1, UInt64(0))
    seed!(r, seed)
    return r
end

function Mixmax17(seed::UInt64 = UInt64(123))
    return Mixmax{17,36,0}(seed)
end

function Mixmax240(seed::UInt64 = UInt64(123))
    return Mixmax{240,32,271828282}(seed)
end

function seed!(r::Mixmax{N,SPECIALMUL,SPECIAL}, seed::UInt64) where {N,SPECIALMUL,SPECIAL}
    sumtot = UInt64(0)
    ovflow = UInt64(0)
    if seed == 0
        throw(DomainError("A non-zero seed is required"))
    end
    l = seed

    for i = 1:N
        l *= UInt64(6364136223846793005)
        l = (l << 32) ⊻ (l >> 32)
        r.v[i] = l & M61
        sumtot += r.v[i]
        if sumtot < r.v[i]
            ovflow + 1
        end
    end
    r.counter = 2
    r.sumtot = modMersenne(modMersenne(sumtot) + (ovflow << 3))
end

#@inline rand(r::Mixmax, ::Type{Float64}) = Float64(rand(r, UInt64)) / M61
@inline rand(r::Mixmax, ::Type{Float64}) = (Float64(rand(r, UInt64) >> 8) + 0.5) /  9007199254740992.0
#actually 61 lowest bits
@inline rand(r::Mixmax, ::Type{UInt64}) = mixmax_get(r)

@inline function mixmax_get(r::Mixmax{N,SPECIALMUL,SPECIAL}) where {N,SPECIALMUL,SPECIAL}
    v = r.v
    if r.counter > N
        iterateRawVec(r)
        r.counter = 3
        return v[2] #first number skipped on purpose
    end
    result = r.v[r.counter]
    r.counter += 1
    result
end

@inline function iterateRawVec(r::Mixmax{N,SPECIALMUL,SPECIAL}) where {N,SPECIALMUL,SPECIAL}
    v = r.v
    temp2 = v[2]
    tempV = r.sumtot
    v[1] = tempV
    sumtotNew = v[1]
    ovflow = UInt64(0)
    tempP = UInt64(0)
    @inbounds for i = 2:N
        tempPO = mulWU(tempP, SPECIALMUL)
        tempP = modadd(tempP, v[i])
        tempV = modMersenne(tempV + tempP + tempPO)
        v[i] = tempV
        sumtotNew += tempV
        if (sumtotNew < tempV)
            ovflow += 1
        end
    end
    temp2 = modMulSpec(r, temp2)
    v[3] = modadd(v[3], temp2)
    sumtotNew += temp2
    if sumtotNew < temp2
        ovflow += 1
    end
    r.sumtot = modMersenne(modMersenne(sumtotNew) + (ovflow << 3))
end

@inline modMulSpec(r::Mixmax{240,SPECIALMUL,SPECIAL}, k::UInt64) where {SPECIALMUL,SPECIAL} =
    fmodmulM61(UInt64(0), UInt64(SPECIAL), k)

@inline function iterateRawVec(r::Mixmax{N,SPECIALMUL,0}) where {N,SPECIALMUL}
    v = r.v
    tempV = r.sumtot
    v[1] = tempV
    sumtotNew = v[1]
    ovflow = UInt64(0)
    tempP = UInt64(0)
    @inbounds for i = 2:N
        tempPO = mulWU(tempP, SPECIALMUL)
        tempP = modadd(tempP, v[i])
        tempV = modMersenne(tempV + tempP + tempPO)
        v[i] = tempV
        sumtotNew += tempV
        if (sumtotNew < tempV)
            ovflow += 1
        end
    end
    r.sumtot = modMersenne(modMersenne(sumtotNew) + (ovflow << 3))
end

@inline function modMersenne(k::UInt64)::UInt64
    return (k & M61) + (k >> BITS)
end

@inline function modadd(foo::UInt64, bar::UInt64)::UInt64
    k = (foo + bar)
    return modMersenne(k)
end

@inline function mulWU(k::UInt64, SPECIALMUL::Int)::UInt64
    return (((k) << (SPECIALMUL) & M61) | ((k) >> (BITS - SPECIALMUL)))
end

@inline function mod128(s::UInt128)::UInt64
    s1 = ((s % UInt64) & M61) + (((s >> 64) % UInt64) * 8) + ((s % UInt64) >> BITS)
    return modMersenne(s1)
end

@inline function fmodmulM61(cum::UInt64, a::UInt64, b::UInt64)::UInt64
    temp = UInt128(a) * UInt128(b) + cum
    return mod128(temp)
end

@inline function skip(gen::Mixmax{N,SPECIALMUL,SPECIAL}, s::UInt64) where {N,SPECIALMUL,SPECIAL}
    n = s + gen.counter - 2
    k = n ÷ (N - 1)
    rk = n - k * (N - 1)
    nReal = k * N + rk + 1 #we skip one every 17.nReal is the effective n without skipping
    nMat = nReal ÷ N
    #binary decomp of nMat
    _skipInternal(gen, nMat)
    gen.sumtot = 0
    for i = 1:N
        gen.sumtot = modadd(gen.sumtot, gen.v[i])
    end
    gen.counter = rk + 2
    gen
end

@inline function _skipInternal(gen::Mixmax{N,SPECIALMUL,SPECIAL}, nMat::UInt64) where {N,SPECIALMUL,SPECIAL}
    iStart = ceil(Int, log2(N))
    for i = 0:(iStart-1)
        if (nMat >> i) & UInt64(1) == UInt64(1)
            m = UInt64(1) << i
            for j=1:m
                iterateRawVec(gen)
            end
        end
    end
    F=GF(2^61-1)
    R,y = PolynomialRing(F,"y")
    local polyCoeffs
    if N == 17
        polyCoeffs = polyCoeffs17
    elseif N==240
        polyCoeffs = polyCoeffs240
    else
        throw(DomainError(string("coefficients for characteristic polynomial with N=",N," not found")))
    end
    cp = R(F.(polyCoeffs))
    eArray = zeros(Int, N)
    gCoeffs = zeros(F, 2^iStart)
    gCoeffs[2^(iStart-1)+1] = F(1)
    g = R(gCoeffs)
    iEnd = floor(Int,log2(nMat))
    @inbounds for i = iStart:iEnd
        g = g^2 % cp
        if (nMat >> i) & UInt64(1) == UInt64(1)
            #println("i ",i, " end ",iEnd," nMat ",nMat )
            #eArray = skipMat17[i-4]
             @inbounds for j=1:N
                  eArray[j] = g.coeffs[j].d
             end
            _skipWithCoeffs(gen, eArray)
        end
    end
end

@inline function _skipInternal(gen::Mixmax{17,SPECIALMUL,SPECIAL}, nMat::UInt64) where {SPECIALMUL,SPECIAL}
for i = 0:63
    if (nMat >> i) & UInt64(1) == UInt64(1)
        #println("i ",i, " nMat ",nMat )
        skipTwoPower(gen, i)
    end
end
end

@inline function skipTwoPower(gen::Mixmax{17,SPECIALMUL,SPECIAL}, m::Int64) where {SPECIALMUL,SPECIAL}
    if m <= 4 #2^(4+1) > N
        m2 = UInt64(1) << m
        for i = 1:m2
            iterateRawVec(gen)
        end
    else
        local skipMat = skipMat17
        if m - 4 > length(skipMat)
            throw(DomainError(string("Can not skip, 2^m too large where m=", m)))
        end
        eArray = skipMat17[m-4]
        _skipWithCoeffs(gen, eArray)
    end
end

@inline function _skipWithCoeffs(gen::Mixmax{N,SPECIALMUL,SPECIAL},eArray::Array{Int}) where {N, SPECIALMUL,SPECIAL}
    a2p = zeros(UInt64, N)
    @inbounds for i = 1:N
        @inbounds for j = 1:N
            a2p[j] = fmodmulM61(a2p[j], UInt64(eArray[i]), gen.v[j])
        end
        iterateRawVec(gen)
    end
    gen.v = a2p[1:N]
    gen.sumtot = 0
    for i = 1:N
        gen.sumtot = modadd(gen.sumtot, gen.v[i])
    end
end

function computeMixmaxCharPoly(gen::Mixmax{N,SPECIALMUL,SPECIAL}) where {N,SPECIALMUL,SPECIAL}
    # M=makeAsonovRealMatrix(Int64, 240, 271828282, 2^32+1)
    # M=makeAsonovRealMatrix(Int64, 240, 487013230256099140, 2^51+1)
    # M=makeAsonovRealMatrix(Int64, 256, 487013230256099064, 1)
    # M =makeAsonovRealMatrix(BigInt, 17, 0, 2^36+1)
    M =makeAsonovMatrix(BigInt, N, SPECIAL, 2^SPECIALMUL + 1)
    F=GF(2^61-1)
    R,y = PolynomialRing(F,"y")
    g = charpoly(R, matrix(F,M))
    return g
end

function makeAsonovMatrix(C, N::Int, s::T, m::T) where {T}
    A = zeros(C,N,N)
    for i=1:N
        for j=1:N
            A[i,j] += 1
        end
    end
    for i=2:N
        for j = 2:i
            A[i,j] = (i-j)*m + 2
        end
    end
    A[3,2] += s
    return A
end

include("mixmax17_skip.jl")
include("mixmax240_skip.jl")


text:
using dopingFVM
using Documenter

makedocs(;
    modules=[dopingFVM],
    authors="Felipe A. Díaz",
    repo="https://github.com/TTp95/dopingFVM.jl/blob/{commit}{path}#L{line}",
    sitename="dopingFVM.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://TTp95.github.io/dopingFVM.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/TTp95/dopingFVM.jl",
)


text:
# sample sequences of ergms with different parameters' values from R package ergm
# and test the PseudoLikelihoodScoreDrivenERGM filter
using Utilities,AReg,StaticNets,JLD,MLBase,StatsBase,CSV, RCall
using PyCall; pygui(); using PyPlot

using JLD2,Utilities, GLM
## load R MCMC simulation and estimates and estimate sdergmTest

Nsample = 100
 targetAllTv = true
 dgpType = "AR"
 T = 50
 N = 50
 Nterms = 2
 Nsteps1 ,Nsteps2 = 2,2
 load_fold = "./data/estimatesTest/sdergmTest/gas_MCMC_comparison_estimates/"
 @load(load_fold *"test_Nodes_$(N)_T_$(T)_Sample_$(Nsample)_Ns_" * dgpType * "_$(Nsteps1)_$(Nsteps2)_MPLE_target_$(targetAllTv).jld" ,
     stats_T, changeStats_T,estParSS_T, parDgpT,
     Nsample,T,N,filtPar_T_Nsample,gasParEst,convFlag,pVals_Nsample,scoreAutoc_Nsample,staticEst)

pValTh = 0.05
 tmp = sum( pVals_Nsample.<pValTh,2)
 (-tmp[1] + tmp[2] + 100)/200
#compute the Average MSE
MSE_sample  = zeros(Nterms,3,Nsample)
 for n=1:Nsample
     if sum(isnan(filtPar_T_Nsample[:,:,n]))==0
    ssDiff = ssFiltPar = estParSS_T[n,:,:]' .-parDgpT
    constDiff =   staticEst[n] .-parDgpT
    gasDiff = filtPar_T_Nsample[:,:,n].-parDgpT
    #MSE_sample[:,1,n] = mean(constDiff.^2,2)
    MSE_sample[:,2,n] = mean(ssDiff.^2,2)
    MSE_sample[:,3,n] = mean(gasDiff.^2,2)
    else
        pVals_Nsample[:,n] = 0.1
    end

 end
 avgMSE = meanSq(MSE_sample,3)
 display(avgMSE)

close()
 figure(figsize=(9,7))
 labSize = 26
 tickSize = 20
  titSize = 28
 legSize = 20
 legTex = ["SD-ERGM "; " ERGM "; "DGP"]
 legTex1 = ["Constant ERGM";"SDERGM <MSE> = $(avgMSE[2,1])"; "Sequence of ERGM = $(avgMSE[2,1])"; "DGP"]
 legTex2 = ["Constant ERGM";"SDERGM <MSE> = $(avgMSE[2,2])"; "Sequence of ERGM = $(avgMSE[2,1])"; "DGP"]
 for n=1:Nsample
 indTvPar = BitArray([true,true])
 gasFiltPar  = filtPar_T_Nsample[:,:,n]
 parInd = 1
 subplot(1,2,1);#plot(1:T,ones(T)*staticEst[n][parInd],"-b")
                plot(1:T,gasFiltPar[parInd,:],"r")
                plot(1:T,estParSS_T[n,:,parInd],".b",markersize = 2)
                plot(1:T,parDgpT[parInd,:],"k",linewidth=5)

 parInd = 2
 subplot(1,2,2);#plot(1:T,ones(T)*staticEst[n][parInd],"-b")
                plot(1:T,gasFiltPar[parInd,:],"r")
                plot(1:T,estParSS_T[n,:,parInd],".b",markersize = 2)
                plot(1:T,parDgpT[parInd,:],"k",linewidth=5)
 end
 namePar1 = "Number of Links"
 namePar2 = "GWESP"
 subplot(1,2,1);    title(namePar1,fontsize = titSize); legend(legTex,fontsize = legSize,framealpha = 0.95)
 ylim([-3.5,-1.8])
 xticks(fontsize = tickSize)
 yticks(fontsize = tickSize)
 subplot(1,2,2);  title(namePar2,fontsize = titSize); legend(legTex,fontsize = legSize,framealpha = 0.95)
 ylim([-0.5,1.2])

  xticks(fontsize = tickSize)
  yticks(fontsize = tickSize)
   tight_layout()

 #Plotta info relative ai tests
figure()
  pValTh = 0.05
      subplot(2,2,1); parInd = 1;   plt[:hist]((pVals_Nsample[parInd,:]), bins=logspace(minimum(log10(pVals_Nsample[1,:])),maximum(log10(pVals_Nsample[parInd,:])), 20))
      if maximum((pVals_Nsample[parInd,:])) >pValTh
           axvspan((pValTh),maximum((pVals_Nsample[parInd,:])),color = "r",alpha = 0.1);
      end
           xscale("log");
      subplot(2,2,2);  parInd = 2;   plt[:hist]((pVals_Nsample[parInd,:]), bins=logspace(minimum(log10(pVals_Nsample[2,:])),maximum(log10(pVals_Nsample[parInd,:])), 20))
      if maximum((pVals_Nsample[parInd,:])) >pValTh
           axvspan((pValTh),maximum((pVals_Nsample[parInd,:])),color = "r",alpha = 0.1);
      end
           xscale("log");
      subplot(2,2,3); plot((pVals_Nsample[1,:]),(pVals_Nsample[2,:]),".")
      axvline((pValTh));axhline((pValTh))
      xscale("log");yscale("log")
      xlabel("P-Value " * namePar1)
      ylabel("P-Value " * namePar2)
      tight_layout()

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
#
# Important directories
simdir          = "/remote/ceph2/group/gerda/data/mpik/gerda-simulations/gerda-gems-db/pss-dev"
simdir_old      = "/remote/ceph2/group/gerda/data/simulation/gerda-mage-sim"
mapping_file    = simdir_old*"/UTILS/det-data/ged-mapping.json"
parameters_file = simdir_old*"/UTILS/det-data/ged-parameters.json"

#########################################
# Create mapping dictionary ############
sim_to_channel = open(mapping_file) do f
    mapping  = JSON.parse(f)["mapping"]
    channels = Dict()
    for ged in keys(mapping)
        channels[mapping[ged]["sim"]] = [mapping[ged]["channel"], mapping[ged]["name"]]
    end
    return channels
end;

#########################################
# Create mapping dictionary ############
channel_to_name = open(mapping_file) do f
    mapping  = JSON.parse(f)["mapping"]
    channels = Dict()
    for ged in keys(mapping)
        channels[mapping[ged]["channel"]] = mapping[ged]["name"]
    end
    return channels
end;

#########################################
# Create mapping dictionary ############
parameters = open(parameters_file) do f
    temp = JSON.parse(f)
    temp["GTF45"] = temp["GTF45_2"]
    return temp
end;

#########################################
# Electronics Model #####################
function applyElectronics(pulse; Ts = 10e-9, GBP = 2750e6, tau = 180e-6, Kv = 150e3, Cd = 50e-12, Cf = 0.35e-12, Rf = 500e6)

    wop = GBP / (2 * pi * Kv)
    Cmod = Cf + Cd
    wmod = 1.0 / (Rf * Cmod)
    alfa = Cmod / (Cf * GBP)

    b0 = 1.0 / alfa
    a2 = 1.0
    a1 = 1.0 / alfa + wop + wmod
    a0 = 1.0 / (tau * alfa) + wmod*wop

    # then the transfer function in the *Laplace* s-domain looks like this:
    #                       b0
    #   T(s) = ----------------------------
    #             a2 * s^2 + a1 * s + a0

    # PolynomialRatio needs z-transform paramters: s- and z-domains can be connected by
    # the bilinear transform:
    #        z - 1
    # s = K -------- , K = Ts/2  , Ts - sampling period
    #        z + 1
    #
    # we can then convert T(s) to T(z):
    #              bz2 * z^2 + bz1 * z + bz0
    #   T(z) = -------------------------------
    #              az2 * z^2 + az1 * z + az0
    #

    K = 2/Ts

    az2 = 1.0   # normalized
    az1 = (2*a0 - 2*K^2)/(K^2 + a1*K + a0)
    az0 = (K^2 - a1*K + a0)/(K^2 + a1*K + a0)

    bz2 = b0/(K^2 + a1*K + a0)
    bz1 = 2*b0/(K^2 + a1*K + a0)
    bz0 = b0/(K^2 + a1*K + a0)

    myfilter = PolynomialRatio([bz2, bz1, bz0], [az2, az1, az0])

    filtered = filt(myfilter, vcat([0], diff(pulse)))

end

#########################################
# Extend simulated waveforms ############
function add_baseline_and_extend_tail(wv::RadiationDetectorSignals.RDWaveform{T,U,TV,UV}, n_baseline_samples::Int, total_waveform_length::Int) where {T,U,TV,UV}
    new_signal::Vector{eltype(UV)} = Vector{eltype(UV)}(undef, total_waveform_length)
    new_signal[1:n_baseline_samples] .= zero(eltype(UV))
    if length(wv.value) <= total_waveform_length - n_baseline_samples
        new_signal[n_baseline_samples+1:n_baseline_samples+length(wv.value)] = wv.value
        new_signal[n_baseline_samples+length(wv.value)+1:end] .= wv.value[end]
    else
        new_signal[n_baseline_samples+1:end] = wv.value[1:total_waveform_length - n_baseline_samples]
    end
    new_times = if TV <: AbstractRange
        range( zero(first(wv.time)), step = step(wv.time), length = total_waveform_length )
    else
        error("Not yet definted for timestamps of type `$(TV)`")
    end
    return RDWaveform( new_times, new_signal )
end

#########################################
# Read / Create h5 sim file #############
function readcreateh5(filepath::String, output_path::String; read::Bool=true, overwrite::Bool=false,
    max_refine::Int=4)
    config_filename = filepath
    ssd_config_file = filepath
    #output_path *= split(basename(filepath), ".")[1]*"/"
    if !isdir(output_path)
        mkdir(output_path)
    end
    h5_sim_file = joinpath(output_path, split(basename(config_filename), ".config")[1]*".h5")

    T = Float32
    is_file = isfile(h5_sim_file)

    if overwrite
        is_file = false
    end

    s = if is_file
        if read
                @info("File already exist. Read in starts..")

                s = HDF5.h5open(h5_sim_file, "r") do h5f
                    Simulation(LegendHDF5IO.readdata(h5f, "SSD_Simulation"));
            end;
            SSD.set_charge_drift_model!(s, ADLChargeDriftModel(T = SSD.get_precision_type(s.detector)));
            return s
        else
            @info("File already exist. Read = false")
        end
    else
        @info("File does not exist. Simulation starts..")
        s = Simulation{T}(ssd_config_file)

        calculate_electric_potential!(  s,
                                        convergence_limit = 1e-7,
                                        use_nthreads = Base.Threads.nthreads(),
                                        max_refinements = max_refine,
                                        max_n_iterations = 500000,
                                        init_grid_spacing = (T(1e-3), T(2e-3), T(2e-3)),
                                        refinement_limits = (T(1e-5), T(1e-1), T(1e-5)),
                                        min_grid_spacing  = (T(2e-6), T(1e-1), T(2e-6)),
                                        depletion_handling = true )
        SSD.calculate_electric_field!(s, n_points_in_φ = 72)
        SSD.set_charge_drift_model!(s, ADLChargeDriftModel(T = T))
        SSD.apply_charge_drift_model!(s)
        for contact in s.detector.contacts
            calculate_weighting_potential!( s, contact.id,
                                            convergence_limit = 1e-7,
                                            use_nthreads = Base.Threads.nthreads(),
                                            max_refinements = 4,
                                            max_n_iterations = 500000,
                                            init_grid_spacing = (T(1e-3), T(2e-3), T(2e-3)),
                                            refinement_limits = (T(1e-4), T(1e-1), T(1e-4)),
                                            min_grid_spacing  = (T(2e-6), T(1e-1), T(2e-6)),
                                            n_points_in_φ = 2,
                                            verbose = true )
        end
        # IO
        HDF5.h5open(h5_sim_file, "w") do h5f
            LegendHDF5IO.writedata( h5f, "SSD_Simulation", NamedTuple(s)  )
        end

        return s
    end
end

function read_AE_params(filename::String)
    fit_params = JSON.parsefile(filename)
    tt = Table(detector=[], set=[], energy=[], scale=[], sig=[], mu=[], C=[], l=[], e=[], t=[], scale_err=[], sig_err=[], mu_err=[], C_err=[], l_err=[], e_err=[], t_err=[])
    for det in keys(fit_params) # Detector
        for set in keys(fit_params[det])
            for E in keys(fit_params[det][set])
                if "bat_params" in keys(fit_params[det][set][E])
                    append!(tt, Table(
                            detector=[det], set=[set], energy=[parse(Int64, E)],
                            scale=[fit_params[det][set][E]["bat_params"][1][1]],
                            sig=[fit_params[det][set][E]["bat_params"][1][2]],
                            mu=[fit_params[det][set][E]["bat_params"][1][3]],
                            C=[fit_params[det][set][E]["bat_params"][1][4]],
                            l=[fit_params[det][set][E]["bat_params"][1][5]],
                            e=[fit_params[det][set][E]["bat_params"][1][6]],
                            t=[fit_params[det][set][E]["bat_params"][1][7]],
                            scale_err=[fit_params[det][set][E]["bat_uncert"][1][1]],
                            sig_err=[fit_params[det][set][E]["bat_uncert"][1][2]],
                            mu_err=[fit_params[det][set][E]["bat_uncert"][1][3]],
                            C_err=[fit_params[det][set][E]["bat_uncert"][1][4]],
                            l_err=[fit_params[det][set][E]["bat_uncert"][1][5]],
                            e_err=[fit_params[det][set][E]["bat_uncert"][1][6]],
                            t_err=[fit_params[det][set][E]["bat_uncert"][1][7]]
                        )
                    );
                end
            end
        end
    end
    return tt
end

function get_cut_values(fit_params, ch, set)
    fit_cut_values = []
    fit_survival =[]
    fit_survival_err =[]
    for cut in keys(fit_params[channel_to_name[ch]]["DEP_cuts"][set])
        if !(cut in ["cut_sum", "cut_fit"])
            push!(fit_cut_values, parse(Float64, cut))
            push!(fit_survival, fit_params[channel_to_name[ch]]["DEP_cuts"][set][cut]["survival"])
            push!(fit_survival_err, fit_params[channel_to_name[ch]]["DEP_cuts"][set][cut]["survival_err"])
        end
    end
    sum_cut_values = fit_params[channel_to_name[ch]]["DEP_cuts"][set * "_sum"]["cut_values"]
    sum_survival = fit_params[channel_to_name[ch]]["DEP_cuts"][set * "_sum"]["survival"]
    sum_survival_err = fit_params[channel_to_name[ch]]["DEP_cuts"][set * "_sum"]["survival_err"]
    return fit_cut_values, fit_survival, fit_survival_err, sum_cut_values, sum_survival, sum_survival_err
end
function get_cut(fit_params, ch, set)
    return fit_params[channel_to_name[ch]]["DEP_cuts"][set]["cut_fit"], fit_params[channel_to_name[ch]]["DEP_cuts"][set]["cut_sum"]
end

function find_intersect(samples, threshold::Real, noisefilter = 1)
    if noisefilter < 1 noisefilter = 1 end
    if length(samples) < 1 error("Empty array") end
    i = 1
    x = samples[i]
    counter = 0
    intersect = i
    if x != threshold
        findHigher = x < threshold
        while i < length(samples)
            i += 1
            x = samples[i]
            if (findHigher && (x >= threshold)) || (!findHigher && (x <= threshold))
                if counter == 0
                    intersect = i
                end
                if counter >= noisefilter-1
                    return intersect
                else
                    counter += 1
                end
            else
                counter = 0
            end
        end
        #println("No intersect found")
        return NaN
    end
    return NaN
end

simdir_old      = "/remote/ceph2/group/gerda/data/simulation/gerda-mage-sim"
mapping_file    = "$simdir_old/UTILS/det-data/ged-mapping.json"
parameters_file = "$simdir_old/UTILS/det-data/ged-parameters.json"

#>------- Moving Window Average -----------------------<#
#>-----------------------------------------------------<#

function movingaverage(X::Vector,numofele::Int)
    BackDelta = div(numofele,2)
    ForwardDelta = isodd(numofele) ? div(numofele,2) : div(numofele,2) - 1
    len = length(X)
    Y = similar(X)
    for n = 1:len
        lo = max(1,n - BackDelta)
        hi = min(len,n + ForwardDelta)
        Y[n] = mean(X[lo:hi])
    end
    return Y
end

function movingaverage(X::Vector,numofele::Int,iterations::Int)
    i = 1
    while i <= iterations
        X = movingaverage(X,numofele)
        i += 1
    end
    return X
end

function movingaverage(X::Vector,numofele::Int,BackDelta::Int64,ForwardDelta::Int64)
    len = length(X)
    Y = similar(X)
    for n = 1:len
        lo = max(1,n - BackDelta)
        hi = min(len,n + ForwardDelta)
        Y[n] = mean(X[lo:hi])
    end
    return Y
end

function movingaverage(X::Vector,numofele::Int,BackDelta::Int64,ForwardDelta::Int64,iterations::Int)
    i = 1
    while i <= iterations
        X = movingaverage(X,numofele,BackDelta,ForwardDelta)
        i += 1
    end
    return X
end

function mwa(samples::Vector, window::Integer)
#window length is in number of samples and has to be an odd number

    if !isodd(window)
        error("For this method the window size has to be an odd integer")
    end

    newsamples = zeros(eltype(samples), length(samples))
    p = convert(eltype(window), (window - 1) / 2)

	#padding before and after the original vector
    samples = vcat(ones(eltype(samples), p) .* samples[1],
		samples, ones(eltype(samples), p) .* samples[length(samples)])

    newsamples[1] = sum(samples[1:window])
    for i in 2:length(newsamples)
        newsamples[i] = newsamples[i-1] + samples[i+2*p] - samples[i-1]
    end

    newsamples./window #return after normalization

end

function multi_mwa(samples::Vector, window::Integer, reps::Integer)
    for i in 1:1:reps
        samples = mwa(samples, window)
    end
    return samples
end

function mwa(samples::Vector, window::Integer, reps::Integer)
    for i in 1:1:reps
        samples = mwa(samples, window)
    end
    return samples
end

#>------- Get average maximum -------------------------<#
#>-----------------------------------------------------<#
function get_avg_maximum(pulse, intervall)
    index = findall(x -> x == maximum(pulse), pulse)[1]
    intervall = (index - Int(floor(intervall/2))+1):1:(index + Int(ceil(intervall/2)))
    return mean(pulse[intervall])
end

#>------- Create mapping dictionary -------------------<#
#>-----------------------------------------------------<#
sim_to_channel = open(mapping_file) do f
    mapping  = JSON.parse(f)["mapping"]
    channels = Dict()
    for ged in keys(mapping)
        channels[mapping[ged]["sim"]] = [mapping[ged]["channel"], mapping[ged]["name"]]
    end
    return channels
end;

#>------- Create mapping dictionary -------------------<#
#>-----------------------------------------------------<#
channel_to_name = open(mapping_file) do f
    mapping  = JSON.parse(f)["mapping"]
    channels = Dict()
    for ged in keys(mapping)
        channels[mapping[ged]["channel"]] = mapping[ged]["name"]
    end
    return channels
end;
ctn  = deepcopy(channel_to_name)

#>------- Create mapping dictionary -------------------<#
#>-----------------------------------------------------<#
parameters = open(parameters_file) do f
    temp = JSON.parse(f)
    temp["GTF45"] = temp["GTF45_2"]
    return temp
end;

#>------- Create Coax dictionary ----------------------<#
#>-----------------------------------------------------<#
channel_to_coax = open(mapping_file) do f
    mapping  = JSON.parse(f)["mapping"]
    channels = Dict()
    for ged in keys(mapping)
        if mapping[ged]["channel"] in [8,9,10,27,28,29,36,37,38,39]
            channels[mapping[ged]["channel"]] = true
        else
            channels[mapping[ged]["channel"]] = false
        end
    end
    return channels
end;

#>------- Create BEGe dictionary ----------------------<#
#>-----------------------------------------------------<#
channel_to_bege = open(mapping_file) do f
    mapping  = JSON.parse(f)["mapping"]
    channels = Dict()
    for ged in keys(mapping)
        if mapping[ged]["channel"] in [8,9,10,27,28,29,36,37,38,39]
            channels[mapping[ged]["channel"]] = false
        else
            channels[mapping[ged]["channel"]] = true
        end
    end
    return channels
end;
ctb = deepcopy(channel_to_bege)

function check_parameters(tt, td, GBP, tau, Cd, Cf)
    cut_events = Table( waveform    = [],
                        evtno       = [],
                        detno       = [],
                        det_ch      = [],
                        edep        = []);
    for i in 1:1:size(tt,1)
        pulse = tt[i].waveform.value
        filtered_pulse = applyElectronics(pulse; Ts = 10e-9, GBP = GBP, tau = tau, Kv = 150e3, Cd = Cd, Cf = Cf, Rf = 500e6)
        inter = find_intersect(filtered_pulse, maximum(filtered_pulse)/2, 5)
        if inter + 399 < length(pulse) && inter - 400 > 0
            cut  = filtered_pulse[(inter-400):1:(inter+400-1)]
            time = (0:10:7990)u"ns"
            cut_pulse = RDWaveform(time, cut)

            append!(cut_events, Table(  waveform    = [cut_pulse],
                                        evtno       = [tt[i].evtno],
                                        detno       = [tt[i].detno],
                                        det_ch      = [tt[i].det_ch],
                                        edep        = [tt[i].edep]))
        end
    end
    if size(cut_events,1) > 0
        sp = zeros(length(cut_events[1].waveform.value))
        for i in 1:1:size(cut_events,1)
            sp .+= cut_events[i].waveform.value
        end
        sp ./= size(cut_events,1)
        sp ./= maximum(sp)
        data_sp = td[1].superpulse.value ./ maximum(td[1].superpulse.value)
        #chi2_a = sum((sp[1:400] .- data_sp[1:400]).^2 ./ data_sp[1:400] )
        #chi2_b = sum((sp[401:end] .- data_sp[401:end]).^2 ./ data_sp[401:end] )
        chi2 = sum(abs.((sp .- data_sp).^2 ./ data_sp) )
        #return abs(chi2_a) + abs(chi2_b), [GBP, tau, Cd, Cf]
        return chi2, [GBP, tau, Cd, Cf]
    else
        return false
    end
end

function check_parameters_AE(tt, td, GBP, tau, Cd, Cf, ch; std_AE=2)
    cut_events = Table( waveform    = [],
                        evtno       = [],
                        detno       = [],
                        det_ch      = [],
                        edep        = []);
    for i in 1:1:size(tt,1)
        pulse = tt[i].waveform.value
        filtered_pulse = applyElectronics(pulse; Ts = 10e-9, GBP = GBP, tau = tau, Kv = 150e3, Cd = Cd, Cf = Cf, Rf = 500e6)
        inter = find_intersect(filtered_pulse, maximum(filtered_pulse)/2, 5)
        if inter + 399 < length(pulse) && inter - 400 > 0
            cut  = filtered_pulse[(inter-399):1:(inter+400)]
            time = (0:10:7990)u"ns"
            cut_pulse = RDWaveform(time, cut)

            append!(cut_events, Table(  waveform    = [cut_pulse],
                                        evtno       = [tt[i].evtno],
                                        detno       = [tt[i].detno],
                                        det_ch      = [tt[i].det_ch],
                                        edep        = [tt[i].edep]))
        end
    end
    if size(cut_events,1) > 0
        sp = zeros(length(cut_events[1].waveform.value))
        if channel_to_bege[ch]
            #> With A/E filter
            A = []
            for i in 1:1:size(cut_events, 1)
                diff_pulse = diff(cut_events.waveform[i].value ./sum(cut_events.edep[i]).val) # normalized to energy
                a = maximum(diff_pulse)
                push!(A, a)
            end
            µ = mean(A)
            sig = std(A)
            good = []
            for i in 1:1:size(cut_events,1)
                if A[i] >= µ - std_AE*sig && A[i] <= µ + std_AE*sig
                    sp  .+= cut_events[i].waveform.value
                    push!(good, i)
                end
            end
            sp ./= length(good)
        else
            #> Normal superposition of all pulses
            for i in 1:1:size(cut_events,1)
                sp .+= cut_events[i].waveform.value
            end
            sp ./= size(cut_events,1)
        end
        sp ./= maximum(sp)
        data_sp = td[1].superpulse.value ./ maximum(td[1].superpulse.value)
        diff_data = diff(data_sp)
        diff_sim  = diff(sp)
        #chi2_a = sum((sp[1:400] .- data_sp[1:400]).^2 ./ data_sp[1:400] )
        #chi2_b = sum((sp[401:end] .- data_sp[401:end]).^2 ./ data_sp[401:end] )
        non_zero = findall(x->x != 0, diff_data)
        chi2 = sum( abs.((diff_sim[non_zero] .- diff_data[non_zero]).^2 ./ diff_data[non_zero]) )
        return chi2, [GBP, tau, Cd, Cf]
    else
        return false
    end
end

function plot_param_set(tt, td, params, ch; std_AE=2)
    cut_events = Table( waveform    = [],
                        evtno       = [],
                        detno       = [],
                        det_ch      = [],
                        edep        = []);
    for i in 1:1:size(tt,1)
        pulse = tt[i].waveform.value
        filtered_pulse = applyElectronics(pulse; Ts = 10e-9, GBP = params[1], tau = params[2], Kv = 150e3, Cd = params[3], Cf = params[4], Rf = 500e6)
        inter = find_intersect(filtered_pulse, maximum(filtered_pulse)/2, 5)
        cut  = filtered_pulse[(inter-399):1:(inter+400)]
        time = (0:10:7990)u"ns"
        cut_pulse = RDWaveform(time, cut)

        append!(cut_events, Table(  waveform    = [cut_pulse],
                                    evtno       = [tt[i].evtno],
                                    detno       = [tt[i].detno],
                                    det_ch      = [tt[i].det_ch],
                                    edep        = [tt[i].edep]))
    end

    sp = zeros(length(cut_events[1].waveform.value))
    if channel_to_bege[ch]
        #> With A/E filter
        A = []
        for i in 1:1:size(cut_events, 1)
            diff_pulse = diff(cut_events.waveform[i].value ./sum(cut_events.edep[i]).val) # normalized to energy
            a = maximum(diff_pulse)
            push!(A, a)
        end
        µ = mean(A)
        sig = std(A)
        good = []
        for i in 1:1:size(cut_events,1)
            if A[i] >= µ - std_AE*sig && A[i] <= µ + std_AE*sig
                sp  .+= cut_events[i].waveform.value
                push!(good, i)
            end
        end
        sp ./= length(good)
    else
        #> Normal superposition of all pulses
        for i in 1:1:size(cut_events,1)
            sp .+= cut_events[i].waveform.value
        end
        sp ./= size(cut_events,1)
    end
    sp ./= maximum(sp)
    data_sp = td[1].superpulse.value ./ maximum(td[1].superpulse.value);
    diff_data = diff(data_sp);
    diff_sim  = diff(sp);

    chi2 = sum( abs.((diff_sim .- diff_data).^2 ./ diff_data) )
    p1 = plot(sp, label="Simulation: chi² of diff = "*string( round(chi2, digits=2) ));
    p1 = plot!(data_sp, color=:black, label="Superpulse (data)");
    p1 = plot!(title=channel_to_name[ch]*" | get electronics");
    p2 = plot(diff_data, label="Diff(data)");
    p2 = plot!(diff_sim, label="Diff(sim)");
    display(plot(p1,p2, layout=(2,1), xlim=[300, 500]))

end

function applyElectronics(pulse; Ts = 10e-9, GBP = 2750e6, tau = 180e-6, Kv = 150e3, Cd = 50e-12, Cf = 0.35e-12, Rf = 500e6)

    wop = GBP / (2 * pi * Kv)
    Cmod = Cf + Cd
    wmod = 1.0 / (Rf * Cmod)
    alfa = Cmod / (Cf * GBP)

    b0 = 1.0 / alfa
    a2 = 1.0
    a1 = 1.0 / alfa + wop + wmod
    a0 = 1.0 / (tau * alfa) + wmod*wop

    # then the transfer function in the *Laplace* s-domain looks like this:
    #                       b0
    #   T(s) = ----------------------------
    #             a2 * s^2 + a1 * s + a0

    # PolynomialRatio needs z-transform paramters: s- and z-domains can be connected by
    # the bilinear transform:
    #        z - 1
    # s = K -------- , K = Ts/2  , Ts - sampling period
    #        z + 1
    #
    # we can then convert T(s) to T(z):
    #              bz2 * z^2 + bz1 * z + bz0
    #   T(z) = -------------------------------
    #              az2 * z^2 + az1 * z + az0
    #

    K = 2/Ts

    az2 = 1.0   # normalized
    az1 = (2*a0 - 2*K^2)/(K^2 + a1*K + a0)
    az0 = (K^2 - a1*K + a0)/(K^2 + a1*K + a0)

    bz2 = b0/(K^2 + a1*K + a0)
    bz1 = 2*b0/(K^2 + a1*K + a0)
    bz0 = b0/(K^2 + a1*K + a0)

    myfilter = PolynomialRatio([bz2, bz1, bz0], [az2, az1, az0])

    filtered = filt(myfilter, vcat([0], diff(pulse)))

end

function add_baseline_and_extend_tail(wv::RadiationDetectorSignals.RDWaveform{T,U,TV,UV}, n_baseline_samples::Int, total_waveform_length::Int) where {T,U,TV,UV}
    new_signal::Vector{eltype(UV)} = Vector{eltype(UV)}(undef, total_waveform_length)
    new_signal[1:n_baseline_samples] .= zero(eltype(UV))
    if length(wv.value) <= total_waveform_length - n_baseline_samples
        new_signal[n_baseline_samples+1:n_baseline_samples+length(wv.value)] = wv.value
        new_signal[n_baseline_samples+length(wv.value)+1:end] .= wv.value[end]
    else
        new_signal[n_baseline_samples+1:end] = wv.value[1:total_waveform_length - n_baseline_samples]
    end
    new_times = if TV <: AbstractRange
        range( zero(first(wv.time)), step = step(wv.time), length = total_waveform_length )
    else
        error("Not yet definted for timestamps of type `$(TV)`")
    end
    return RDWaveform( new_times, new_signal )
end


text:
module MultilevelCoordinateSearch

using StaticArrays, DataStructures
using Compat

include("types.jl")
include("init.jl")
include("util.jl")

end # module


text:
# Julia wrapper for header: /usr/include/scip/pub_paramset.h
# Automatically generated using Clang.jl wrap_c

function SCIPparamGetType(param)
    ccall((:SCIPparamGetType, libscip), SCIP_PARAMTYPE, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetName(param)
    ccall((:SCIPparamGetName, libscip), Cstring, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetDesc(param)
    ccall((:SCIPparamGetDesc, libscip), Cstring, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetData(param)
    ccall((:SCIPparamGetData, libscip), Ptr{SCIP_PARAMDATA}, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamIsAdvanced(param)
    ccall((:SCIPparamIsAdvanced, libscip), UInt32, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamIsFixed(param)
    ccall((:SCIPparamIsFixed, libscip), UInt32, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamSetFixed(param, fixed)
    ccall((:SCIPparamSetFixed, libscip), Cvoid, (Ptr{SCIP_PARAM}, UInt32), param, fixed)
end

function SCIPparamGetBool(param)
    ccall((:SCIPparamGetBool, libscip), UInt32, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetBoolDefault(param)
    ccall((:SCIPparamGetBoolDefault, libscip), UInt32, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetInt(param)
    ccall((:SCIPparamGetInt, libscip), Cint, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetIntMin(param)
    ccall((:SCIPparamGetIntMin, libscip), Cint, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetIntMax(param)
    ccall((:SCIPparamGetIntMax, libscip), Cint, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetIntDefault(param)
    ccall((:SCIPparamGetIntDefault, libscip), Cint, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetLongint(param)
    ccall((:SCIPparamGetLongint, libscip), Clonglong, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetLongintMin(param)
    ccall((:SCIPparamGetLongintMin, libscip), Clonglong, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetLongintMax(param)
    ccall((:SCIPparamGetLongintMax, libscip), Clonglong, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetLongintDefault(param)
    ccall((:SCIPparamGetLongintDefault, libscip), Clonglong, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetReal(param)
    ccall((:SCIPparamGetReal, libscip), Cdouble, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetRealMin(param)
    ccall((:SCIPparamGetRealMin, libscip), Cdouble, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetRealMax(param)
    ccall((:SCIPparamGetRealMax, libscip), Cdouble, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetRealDefault(param)
    ccall((:SCIPparamGetRealDefault, libscip), Cdouble, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetChar(param)
    ccall((:SCIPparamGetChar, libscip), UInt8, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetCharAllowedValues(param)
    ccall((:SCIPparamGetCharAllowedValues, libscip), Cstring, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetCharDefault(param)
    ccall((:SCIPparamGetCharDefault, libscip), UInt8, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetString(param)
    ccall((:SCIPparamGetString, libscip), Cstring, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamGetStringDefault(param)
    ccall((:SCIPparamGetStringDefault, libscip), Cstring, (Ptr{SCIP_PARAM},), param)
end

function SCIPparamIsDefault(param)
    ccall((:SCIPparamIsDefault, libscip), UInt32, (Ptr{SCIP_PARAM},), param)
end


text:
import Base.sign

struct PartialTransposeAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int,Int}
    sys::Int
    dims::Vector

    function PartialTransposeAtom(x::AbstractExpr, sys::Int, dims::Vector)
        if x.size[1] ≠ x.size[2]
            throw(ArgumentError("Only square matrices are supported"))
        end
        if !(1 ≤ sys ≤ length(dims))
            throw(
                ArgumentError(
                    "Invalid system, should between 1 and ",
                    length(dims),
                    "; got ",
                    sys,
                ),
            )
        end
        if x.size[1] ≠ prod(dims)
            throw(
                ArgumentError(
                    "Dimension of system doesn't correspond to dimension of subsystems",
                ),
            )
        end
        children = (x,)
        return new(
            :partialtranspose,
            hash(children),
            children,
            x.size,
            sys,
            dims,
        )
    end
end

function sign(x::PartialTransposeAtom)
    return sign(x.children[1])
end

function curvature(x::PartialTransposeAtom)
    return ConstVexity()
end

function monotonicity(x::PartialTransposeAtom)
    return (Nondecreasing(),)
end

function evaluate(x::PartialTransposeAtom)
    return partialtranspose(evaluate(x.children[1]))
end

"""
    partialtranspose(x, sys::Int, dims::Vector)

Returns the partial transpose of `x` over the `sys`th system, where `dims` is a vector of integers encoding the dimensions of each subsystem.
"""
function partialtranspose(x::AbstractMatrix, sys::Int, dims::Vector)
    if size(x, 1) ≠ size(x, 2)
        throw(ArgumentError("Only square matrices are supported"))
    end
    if !(1 ≤ sys ≤ length(dims))
        throw(
            ArgumentError(
                "Invalid system, should between 1 and $(length(dims)); got $sys",
            ),
        )
    end
    if size(x, 1) ≠ prod(dims)
        throw(
            ArgumentError(
                "Dimension of system doesn't correspond to dimension of subsystems",
            ),
        )
    end
    n = length(dims)
    d = prod(dims)
    s = n - sys + 1
    p = collect(1:2n)
    p[s] = n + s
    p[n+s] = s

    rdims = reverse(dims)
    r = reshape(x, (rdims..., rdims...))
    return reshape(permutedims(r, p), (d, d))
end

"""
    permutedims_matrix(dims, p)

Returns a matrix `M` so that for any vector `v` of length `prod(dims)`,

    M*v == vec(permutedims(reshape(v, dims), p))

"""
function permutedims_matrix(dims, p)
    d = prod(dims)
    n = length(dims)
    return sparse(
        reshape(
            PermutedDimsArray(
                reshape(Matrix(I, d, d), (dims..., dims...)),
                (p..., (n+1:2n)...),
            ),
            (d, d),
        ),
    )
end

# borrowing this from transpose.jl:
# Since everything is vectorized, we simply need to multiply x by a permutation
# matrix such that coeff * vectorized(x) - vectorized(x') = 0
function conic_form!(
    x::PartialTransposeAtom,
    unique_conic_forms::UniqueConicForms,
)
    if !has_conic_form(unique_conic_forms, x)
        objective = conic_form!(x.children[1], unique_conic_forms)

        n = length(x.dims)
        d = prod(x.dims)
        s = n - x.sys + 1
        p = collect(1:2n)
        p[s] = n + s
        p[n+s] = s

        rdims = reverse(x.dims)

        partialtranspose_matrix = permutedims_matrix((rdims..., rdims...), p)

        objective = partialtranspose_matrix * objective
        cache_conic_form!(unique_conic_forms, x, objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

function partialtranspose(x::AbstractExpr, sys::Int, dim::Vector)
    return PartialTransposeAtom(x, sys, dim)
end


text:
ENV["GKSwstype"] = "100" #src

# # Supervised Learning
#
# Machine learning is split into two type of tasks, **supervised** and **unsupervised**
# learning.
#
# - Supervised learning is the family of machine learning tasks that deal with predicting/
#   or classifying an output from an input.
# - Unsupervised learning is the family of machine learning tasks that are used on data
#   that is unlabeled, or does not have an associated output.
#
#
# ## Decision Trees

# ### Bagging

# ### Random Forests

# ## Neural Networks

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
# newton.jl: function to do Newtons method, including calculating the Jacobian

export newton
import Utils.free

include("newton_setup.jl")

@doc """
### NonlinearSolvers.newton

  This function uses Newton's method to reduce the residual.  The Jacobian
  is calculated using one of several available methods.

  This two options that control how the jacobian is calculated and stored are
  `jac_type` and `jac_method`.  `jac_type` = 1, 2, or 3 correspond to a
  dense jacobian matrix, a SparseMatrixCSC matrix that uses direct solves

  jac_type | Meaning
     1        dense matrix
     2        SparseMatrixCSC
     3        Petsc explictly stored matrix (using CSR internally)
     4        Petsc matrix shell (jacobian vector products)

  jac_method | Meaning
     1          use finite differences
     2          use complex step

  All combinations of jac_type and jac_method are supported except
  jac_type = 4 and jac_method = 2.

  The magnitude of perturbation for both finite differences and complex step is
  opts["epsilon"]

  This function supports extensive debugging options. See the document
  describing all options supported by PDESolver for details.

  The initial condition in eqn.q_vec should be in whatever variables
  the residual evaluation uses.

  **Arguments**
    * func  : function that evalutes the residual
    * mesh : mesh to use in evaluating the residual
    * sbp : sbp operator to be used to evaluate the residual
    * eqn : EulerData to use to evaluate the residual
    * opts : options dictonary
    * pmesh : mesh used for preconditioning, defaults to mesh

    **Optional Arguments**

    func must have the signature func(mesh, sbp, eqn, opts, t=0.0)

"""->
function newton(func::Function, mesh::AbstractMesh, sbp::AbstractOperator,
                eqn::AbstractSolutionData, opts, pmesh=mesh, t=0.0)

  jac_type = opts["jac_type"]::Int  # jacobian sparse or dense

  # physicsRhs is defind in residual_evaluation.jl
  rhs_func = physicsRhs
  ctx_residual = (func, )

  pc, lo = getNewtonPCandLO(mesh, sbp, eqn, opts)
  ls = StandardLinearSolver(pc, lo, eqn.comm, opts)

  # allocate jac & rhs, and construct newton_data
  newton_data, rhs_vec = setupNewton(mesh, pmesh, sbp, eqn, opts, ls, alloc_rhs=false)

  newtonInner(newton_data, mesh, sbp, eqn, opts, rhs_func, ls,
              rhs_vec, ctx_residual, t)

  #TODO: return newton_data, don't free it
  free(newton_data)
end

"""
  This function contains the algorithm for Newton's method.  It is intended
  to be used by other functions rather than invoked as a solver directly.
  For example [`newton`](@ref) uses newtonInner to solve steady problems while
  [`crank_nicolson`](@ref) uses it to solve the nonlinear problem arising
  from the time discretization.

  For reference, Newton's method solves the problem R(q) = 0 using
    dR/dq delta q = -R(q)
  where R is the rsidual and q is the solution

  On entry, eqn.q_vec must contain the initial guess for q.  On exit, eqn.q_vec
  will contain the solution to f(q) = 0.  eqn.q will also be consistent with
  eqn.q_vec, as will the send and receive buffers in eqn.shared_data

  On exit, rhs_vec will have the residual corresponding to eqn.q_vec in it,
  with the imaginary part set to zero.

  The [`AbstractLO`](@ref) and [`AbstractPC`](@ref) supplied
  inside the [`LinearSolver`](@ref) object must match the `jac_type`.

  When doing inexact Newton-Krylov, `newonInner` modifies the tolerances
  of the linear solver.  Users calling `newtonInner` repeatedly with the
  same linear solver object should reset the initial tolerances as needed.

  **Inputs:**

   * newton_data: NewtonData object, typically obtained from setupNewton
   * mesh: a mesh object
   * sbp: an SBP operator
   * eqn: a solution data object
   * opts: options dictionary
   * ls: a LinearSolver with the preconditioner and linear operator fully
         initialized.
   * rhs_vec: vector to store R(q) in
   * ctx_residual: extra data required by rhs_func

  The user must supply two functions, one to calculate the residual vector
  (referred to as rhs_vec), and another to compute the Jacobian.

  rhs_func should compute (eqn.q_vec) -> (rhs_vec) and have the signature

    rhs_func(mesh, sbp, eqn, opts, rhs_vec, ctx_residual, t=0.0)

  Note that the ctx_residual passed into newtonInner is passed directly
  to rhs_func.  The contents of ctx_residual can be whatever is needed by
  rhs_func to perform its computation.  See [`physicsRhs`](@ref) for
  the specific requirements on rhs_func and for an example implementation.

  The same ctx_residual passed into newtonInner is passed directly to
  [`calcLinearOperator`](@ref)..

  This function supports jacobian/preconditioner freezing using the
  prefix "newton".

  Aliasing restrictions: None.  In particular, rhs_vec *can* alias eqn.res_vec,
                         and this leads so some efficiency because it avoids
                         needlessly copying data.

"""
function newtonInner(newton_data::NewtonData, mesh::AbstractMesh,
                     sbp::AbstractOperator, eqn::AbstractSolutionData,
                     opts, rhs_func, ls::LinearSolver,
                     rhs_vec, ctx_residual=(),
                     t=0.0;)

  # verbose 4 = when newtonInner is used as an inner method for time marching
  #             or something

  @debug1 println(eqn.params.f, "==== Entered newton")
  @debug1 flush(eqn.params.f)

  myrank = mesh.myrank
  reinitNewtonData(newton_data, mesh, sbp, eqn, opts)
  verbose = newton_data.verbose

  # options
  jac_method = opts["jac_method"]::Int  # finite difference or complex step
  jac_type = opts["jac_type"]::Int  # jacobian sparse or dense
  epsilon = opts["epsilon"]::Float64
#  recalc_prec_freq = opts["recalc_prec_freq"]::Int

  @assert opts["parallel_type"] == 2

  @verbose5 @mpi_master println(BSTDOUT, "\nEntered Newtons Method")
  @verbose5 @mpi_master begin
    println("step_tol = ", newton_data.step_tol)
    println("res_abstol = ", newton_data.res_abstol)
    println("res_reltol = ", newton_data.res_reltol)
  end

  if jac_method == 1  # finite difference
    pert = epsilon
  elseif jac_method == 2  # complex step
    pert = complex(0, epsilon)
  end

  # unpack variables
  m = length(rhs_vec)
  Tsol = typeof(rhs_vec[1])
  res_0 = newton_data.res_0  # residual vector (real)
  res_0_norm = 0.0  # norm of res_0  #TODO: rename this to res_norm
  delta_q_vec = newton_data.delta_q_vec  # newton update

  # evaluating residual at initial condition
  @verbose5 @mpi_master println(BSTDOUT, "evaluating residual at initial condition"); flush(BSTDOUT)
  res_0_norm = rhs_func(mesh, sbp, eqn, opts, rhs_vec, ctx_residual, t)
  recordResNorm(newton_data, res_0_norm)

  # extract the real components to res_0
  for i=1:m
    res_0[i] = real(rhs_vec[i])
  end

  # if the user said to use the first residual for computing relative residuals
  if newton_data.set_rel_norm
    newton_data.res_norm_rel = res_0_norm
  end

  writeFiles(newton_data, mesh, sbp, eqn, opts)
  is_converged = checkConvergence(newton_data)

  # if initial residual satisfies tolerances, return
  if is_converged
    return nothing
  end

  # do Newton's method if not converged
  @verbose5 @mpi_master print(BSTDOUT, "\n")

  #----------------------------------------------------------------------------
  # Start of newton iteration loop
  eqn.params.time.t_newton += @elapsed for i=1:newton_data.itermax

    @verbose5 @mpi_master println(BSTDOUT, "===== newton iteration: ", i)
    newton_data.itr = i

    # recalculate PC and Jacobian if needed
    doRecalculation(newton_data.recalc_policy, i,
                    ls, mesh, sbp, eqn, opts, ctx_residual, t)

    # compute eigs, condition number, etc.
    doMatrixCalculations(newton_data, opts)

    # negate res
    for j=1:m
      res_0[j] = -res_0[j]
    end

    #TODO: DEBUGGING
    @mpi_master println(BSTDOUT, "writing newton_q_vec_$i")
    writeSolutionFiles(mesh, sbp, eqn, opts, "newton_q_vec_$i")
    flush(BSTDOUT)

    # calculate Newton step
    flush(BSTDOUT)
    tsolve = @elapsed linearSolve(ls, res_0, delta_q_vec, verbose)
    eqn.params.time.t_solve += tsolve
    step_norm = norm(delta_q_vec)  #TODO: use calcNorm()
    recordStepNorm(newton_data, step_norm)

    # perform Newton update
    for j=1:m
      eqn.q_vec[j] += newton_data.step_fac*delta_q_vec[j]
    end

    flush(BSTDOUT)
    saveSolutionToMesh(mesh, real(eqn.q_vec))
    fname = string("newton_", i)
    writeVisFiles(mesh, fname)
    flush(BSTDOUT)

    # calculate residual at updated location, used for next iteration rhs
    res_0_norm = rhs_func(mesh, sbp, eqn, opts, rhs_vec, ctx_residual, t)
    recordResNorm(newton_data, res_0_norm)

    flush(BSTDOUT)
    saveSolutionToMesh(mesh, real(eqn.res_vec))
    fname = string("newton_res_", i)
    writeVisFiles(mesh, fname)
    flush(BSTDOUT)

    # extract real component to res_0
    for j=1:m
      res_0[j] = real(rhs_vec[j])
    end

    writeFiles(newton_data, mesh, sbp, eqn, opts)
    is_converged = checkConvergence(newton_data)

    if is_converged
      # remove the imaginary part of rhs_vec before exiting
      for j=1:m
        rhs_vec[j] = real(rhs_vec[j])
      end
      flush(BSTDOUT)

      return nothing
    end

    updateKrylov(newton_data)
    @verbose5 @mpi_master print(BSTDOUT, "\n")

  end  # end loop over newton iterations

  # print final information before exiting (unconverged)
  @mpi_master begin
    println(BSTDOUT, "Warning: Newton iteration did not converge in ", newton_data.itermax, " iterations")
    println(BSTDOUT, "  Final step size: ", newton_data.step_norm_i)
    println(BSTDOUT, "  Final residual: ", res_0_norm)
    println(BSTDOUT, "  Final relative residual: ", res_0_norm/newton_data.res_norm_rel)
    @verbose5 close(newton_data.fconv);
  end
  flush(BSTDOUT)

  # remove the imaginary part of the residual before exiting
  for j=1:m
    rhs_vec[j] = real(rhs_vec[j])
  end
  clearEulerConstants(newton_data.ls)

  return nothing
end               # end of function newton()

"""
  Do a bunch of file IO for Newtons method.  This function is called
  immediately after a residual evaluation is done.  It also calles
  the majorIterationCallback.

  On entry, res_0, eqn.q, eqn.q_vec, and eqn.res should have the values
  from the current iteration of the Newton loop.

  **Inputs**

   * newton_data: the NewtonDAta
   * mesh
   * sbp
   * eqn
   * opts

  **Options Keys**

"""
function writeFiles(newton_data::NewtonData, mesh, sbp, eqn, opts)

  verbose = newton_data.verbose
  myrank = newton_data.myrank
  itr = newton_data.itr

  res_norm = newton_data.res_norm_i
  res_norm_rel = newton_data.res_norm_rel
  step_norm = newton_data.step_norm_i

  @verbose4 @mpi_master begin
    println(BSTDOUT, "residual norm = ", res_norm)
    println(BSTDOUT, "relative residual ", res_norm/res_norm_rel)
  end

  # decide which files to write
  write_q = opts["writeq"]::Bool
  write_rhs = opts["write_rhs"]::Bool
  write_res = opts["write_res"]::Bool
  write_sol = opts["write_sol"]::Bool

  # write starting values for next iteration to file
  if write_sol
    fname = string("q_vec", itr, "_", myrank, ".dat")
    writeSolutionFiles(mesh, sbp, eqn, opts, "q_vec_itr")
  end

  if write_q
    array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec)

    fname = string("q", itr, "_", myrank, ".dat")
    writedlm("q$i_$myrank.dat", eqn.q)
  end

  # write rhs to file
  if write_rhs
    fname = string("rhs", itr, "_", myrank, ".dat")
    writedlm(fname, newton_data.res_0)
  end

  if write_res
    fname = string("res", itr, "_", myrank, ".dat")
    writedlm(fname, eqn.res)
  end

  #TODO: have an option to control this
#    saveSolutionToMesh(mesh, eqn.q_vec)
#    writeVisFiles(mesh, "newton_$i")

  @verbose5 eqn.majorIterationCallback(itr, mesh, sbp, eqn, opts, BSTDOUT)

  # write to convergence file
  @verbose5 @mpi_master begin
    println(newton_data.fconv, itr, " ", res_norm, " ", step_norm)
    println(BSTDOUT, "printed to convergence.dat")
  end

  return nothing
end

"""
  Checks whether the Newton iteration has converged or not and prints
  appropriate messages.

  The residual and step norms for the current iteration must be saved to
  the NewtonData object before this function is called

  **Inputs**

   * newton_data: the NewtonData object

  **Outputs**

   * is_converged: true if converged, false otherwise
"""
function checkConvergence(newton_data::NewtonData)

  verbose = newton_data.verbose
  myrank = newton_data.myrank
  itr = newton_data.itr

  res_norm = newton_data.res_norm_i
  res_norm_rel = newton_data.res_norm_rel
  step_norm = newton_data.step_norm_i

  is_converged = false

  # absolute tolerance
  if res_norm < newton_data.res_abstol
    is_converged = true
    @mpi_master if itr == 0
      println(BSTDOUT, "Initial condition satisfies res_tol with residual norm ", res_norm, " < ", newton_data.res_abstol)
    else
      println(BSTDOUT, "Newton iteration converged with residual norm ", res_norm, " < ", newton_data.res_abstol)
    end
  end

  # relative tolerance
  rel_norm = res_norm/res_norm_rel
  if rel_norm < newton_data.res_reltol
    is_converged = true
    @mpi_master if itr == 0
      println(BSTDOUT, "Initial condition satisfied res_reltol with relative residual ", rel_norm, " < ", newton_data.res_reltol)
    else
      println(BSTDOUT, "Newton iteration converged with relative residual norm ", rel_norm, " < ", newton_data.res_reltol)
    end
  end

  # step tolerance
  if step_norm <= newton_data.step_tol && itr > 0
    println("step tolerance satisfied")
    is_converged = true
    @mpi_master println(BSTDOUT, "Newton iteration converged with step_norm = ", step_norm, " < ", newton_data.step_tol)
    @mpi_master println(BSTDOUT, "Final residual = ", res_norm)
  end

  if is_converged
    @mpi_master if itr == 0
      println(BSTDOUT, "Not entering Newton iteration loop")
    else
      println(BSTDOUT, "Iteration count: ", itr)
    end

    clearEulerConstants(newton_data.ls)
  end

  flush(BSTDOUT)
  return is_converged
end

"""
  Do calculations on the jacobian (eigenvalues, condition number, etc.) if
  requested by options

  **Inputs**

   * newton_data: the NewtonData object
"""
function doMatrixCalculations(newton_data::NewtonData, opts)

  # unpack options
  write_jac = opts["write_jac"]::Bool
  print_cond = opts["print_cond"]::Bool
  print_eigs = opts["print_eigs"]::Bool
  write_eigs = opts["write_eigs"]::Bool
  write_eigdecomp = opts["write_eigdecomp"]::Bool
  jac_type = opts["jac_type"]::Int  # jacobian sparse or dense

  ls = newton_data.ls
  itr = newton_data.itr
  myrank = newton_data.myrank

  # print as determined by options
  if write_jac && jac_type != 4
    jac = getBaseLO(ls.lo).A
    if jac_type == 3
      MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)
      MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)
    end
    writedlm("jacobian$itr.dat", full(jac))
    @mpi_master println(BSTDOUT, "finished printing jacobian"); flush(BSTDOUT)
  end

  # calculate Jacobian condition number
  if print_cond && ( jac_type == 1 || jac_type == 2)
    jac = getBaseLO(ls.lo).A
    println(BSTDOUT, "calculating condition number of jacobian"); flush(BSTDOUT)
    cond_j = cond(full(jac))
    println(BSTDOUT, "Condition number of jacobian = ", cond_j); flush(BSTDOUT);
  end

  # if eigenvalues requested and we can calculate them
  if (( print_eigs || write_eigs) && (jac_type == 1 || jac_type == 2))
    jac = getBaseLO(ls.lo).A
    println(BSTDOUT, "calculating eigenvalues of jacobian")
    flush(BSTDOUT)
    if typeof(jac) <: Array
      eigs_i = reverse(eigvals(jac))
    else
      eigs_i = reverse(eigvals(full(jac)))
    end
    if print_eigs
      println(BSTDOUT, "eigenvalues =")
      for i=1:length(eigs_i)
        println(BSTDOUT, eigs_i[i])
      end
    end
    flush(BSTDOUT)

    if write_eigs
      writedlm("eigs_real$itr.dat", real(eigs_i))
      writedlm("eigs_imag$itr.dat", imag(eigs_i))
    end

  end

  # do the full eigenvalue decomposition
  # if requested and if julia owns the Jacobian matrix
  if write_eigdecomp && ( jac_type == 1 || jac_type == 2)
    jac = getBaseLO(ls.lo).A
    println(BSTDOUT, "doing eigen decomposition"); flush(BSTDOUT)
    # make a dense jacobian so we can get *all* the eigenvalues and vectors
    # the algorithm for sparse matrices can get all - 2 of them, and might
    # have problems for matrices with a wide range of eigenvalues
    jac_dense = full(jac)
    D, V = eig(jac_dense)
    writedlm("eigdecomp_real$itr.dat", real(D))
    writedlm("eigdecomp_imag$itr.dat", imag(D))
#      writedlm("eigdecomp_realvecs$i.dat", real(V))
#      writedlm("eigdecomp_imagvecs$i.dat", imag(V))

    max_val = typemin(Float64)
    min_val = typemax(Float64)
  elseif write_eigdecomp # && we can't calculate it
    println(BSTDERR, "Warning: not performing eigen decomposition for jacobian of type $jac_type")

  end

  return nothing
end


text:
# Tests formulations in prob/opf
@testset "test opf" begin

opf_shunt_ac_objective = [14676.95, 27545.01]
@testset "opf shunt acp - $(i)" for (i,network) in enumerate(networks)

    result = run_opf_shunt(network, ACPPowerModel, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], opf_shunt_ac_objective[i]; atol = 1e0)
end

@testset "opf shunt acr - $(i)" for (i,network) in enumerate(networks)

    result = run_opf_shunt(network, ACRPowerModel, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], opf_shunt_ac_objective[i]; atol = 1e0)
end

opf_shunt_soc_objective = [14666.91, 27510.25]
@testset "opf shunt soc - $(i)" for (i,network) in enumerate(networks)

    result = run_opf_shunt(network, SOCWRPowerModel, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], opf_shunt_soc_objective[i]; atol = 1e0)
end

opf_shunt_dcp_objective = [14642.16, 26982.17]
@testset "opf shunt dcp - $(i)" for (i,network) in enumerate(networks)

    result = run_opf_shunt(network, DCPPowerModel, lp_solver)

    @test isapprox(result["termination_status"], OPTIMAL)
    @test isapprox(result["objective"], opf_shunt_dcp_objective[i]; atol = 1e0)
end

@testset "opf shunt dc - infeasible" begin
    result = run_opf_shunt(network_infeasible, DCPPowerModel, lp_solver)
    @test isapprox(result["termination_status"], INFEASIBLE)
end

opf_cheap_ac_objective = [14676.95, 27542.05]
@testset "opf cheap acp - $(i)" for (i,network) in enumerate(networks)

    result = run_opf_cheap(network, ACPPowerModel, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], opf_cheap_ac_objective[i]; atol = 1e0)
end

@testset "opf cheap acr - $(i)" for (i,network) in enumerate(networks)

    result = run_opf_cheap(network, ACRPowerModel, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], opf_cheap_ac_objective[i]; atol = 1e0)
end

opf_cheap_soc_objective = [14666.81, 27507.30]
@testset "opf cheap soc - $(i)" for (i,network) in enumerate(networks)

    result = run_opf_cheap(network, SOCWRPowerModel, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], opf_cheap_soc_objective[i]; atol = 1e0)
end

# not fully general
# opf_cheap_dcp_objective = [14676.95, 27542.05]
# @testset "opf cheap dcp - $(i)" for (i,network) in enumerate(networks)

#     result = run_opf_cheap(network, ACPPowerModel, lp_solver)

#     @test isapprox(result["termination_status"], OPTIMAL)
#     @test isapprox(result["objective"], opf_cheap_ac_objective[i]; atol = 1e0)
# end

opf_cheap_dc_objective = [14642.16, 26982.17]
@testset "opf cheap_dc dc - $(i)" for (i,network) in enumerate(networks)

    result = run_opf_cheap(network, DCPPowerModel, lp_solver)

    @test isapprox(result["termination_status"], OPTIMAL)
    @test isapprox(result["objective"], opf_cheap_dc_objective[i]; atol = 1e0)
end

@testset "opf shunt dc - infeasible" begin
    result = run_opf_cheap(network_infeasible, DCPPowerModel, lp_solver)
    @test isapprox(result["termination_status"], INFEASIBLE)
end

opf_pg_pf_rect_objective = [122981.99, 9.480588301784407e6]
@testset "opf pg pf rect 5 - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)

    deactivate_rate_a!(network)
    activate_rate_a_violations!(network)

    result = run_opf_cheap_lazy_acr(network, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], opf_pg_pf_rect_objective[i]; atol = 1e0)
end

@testset "opf pg pf rect 5 - infeasible" begin
    network = deepcopy(network_infeasible)

    deactivate_rate_a!(network)
    activate_rate_a_violations!(network)

    result = run_opf_cheap_lazy_acr(network, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_INFEASIBLE)
end

end # close test group


text:
"""
    interval_of_existence(Y::Interval{T}, Z₁::Interval{T}, R::T) where {T<:Real}

Return the interval of existence associated with the polynomial ``Y + (Z_1 - 1)r``.
"""
function interval_of_existence(Y::Interval{T}, Z₁::Interval{T}, R::T) where {T<:Real}
    if inf(Y) < 0 || inf(Z₁) < 0 || R < 0
        return throw(DomainError)
    elseif Y == 0 # exact solution
        if Z₁ ≥ 1
            return Interval(zero(T))
        else
            return Interval(zero(T), R)
        end
    elseif Z₁ ≥ 1
        return emptyinterval(T)
    else
        r = Y/(one(Interval{T}) - Z₁)
        if r ≤ R
            return Interval(sup(r), R)
        else
            return emptyinterval(T)
        end
    end
end

"""
    interval_of_existence(Y::Interval{T}, Z₁::Interval{T}, Z₂::Interval{T}, R::T) where {T<:Real}

Return the interval of existence associated with the polynomial ``Y + (Z_1 - 1)r + Z_2 r^2 / 2``.
"""
function interval_of_existence(Y::Interval{T}, Z₁::Interval{T}, Z₂::Interval{T}, R::T) where {T<:Real}
    if inf(Y) < 0 || inf(Z₁) < 0 || inf(Z₂) < 0 || R < 0
        return throw(DomainError)
    elseif Z₂ == 0
        return interval_of_existence(Y, Z₁, R)
    else
        b = Z₁ - one(Interval{T})
        Δ = b*b - 2Z₂*Y
        if inf(Δ) < 0 # complex roots: d_ = - b - sqrt(complex(Δ)), r₁_ = d_/Z₂, r₂_ = 2Y/d_
            return emptyinterval(T)
        else # real roots
            d = sqrt(Δ) - b
            r₂ = d/Z₂
            if Y == 0 # exact solution
                if Z₁ ≥ 1
                    return Interval(zero(T))
                elseif r₂ ≤ R && r₂ * Z₂ < 1
                    return Interval(zero(T), inf(r₂))
                elseif R * Z₂ < 1
                    return Interval(zero(T), R)
                else
                    Z₂⁻¹ = prevfloat(inf(inv(Z₂)))
                    if Z₂⁻¹ ≤ r₂ && Z₂⁻¹ ≤ R
                        return Interval(zero(T), Z₂⁻¹)
                    else
                        return Interval(zero(T))
                    end
                end
            elseif Z₁ ≥ 1 # Z₁ ≠ 1 (since Y > 0 and Z₁ = 1 implies complex roots)
                return emptyinterval(T)
            else
                r₁ = 2Y/d
                if r₁ ≤ R
                    if r₂ ≤ R && r₂ * Z₂ < 1
                        return Interval(sup(r₁), inf(r₂))
                    elseif R * Z₂ < 1
                        return Interval(sup(r₁), R)
                    else
                        Z₂⁻¹ = prevfloat(inf(inv(Z₂)))
                        if Z₂⁻¹ ≤ r₂ && Z₂⁻¹ ≤ R
                            return Interval(sup(r₁), Z₂⁻¹)
                        elseif r₁ * Z₂ < 1
                            return Interval(sup(r₁))
                        else
                            return emptyinterval(T)
                        end
                    end
                else
                    return emptyinterval(T)
                end
            end
        end
    end
end