text:
module CodecZlib

export
    # gzip
    GzipCompressor,
    GzipCompressorStream,
    GzipDecompressor,
    GzipDecompressorStream,

    # zlib
    ZlibCompressor,
    ZlibCompressorStream,
    ZlibDecompressor,
    ZlibDecompressorStream,

    # deflate
    DeflateCompressor,
    DeflateCompressorStream,
    DeflateDecompressor,
    DeflateDecompressorStream

import TranscodingStreams:
    TranscodingStreams,
    TranscodingStream,
    Memory,
    Error,
    initialize,
    finalize,
    splitkwargs
using Libdl

const libzpath = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(libzpath)
    error("CodecZlib.jl is not installed properly, run Pkg.build(\"CodecZlib\") and restart Julia.")
end
include(libzpath)
check_deps()

include("libz.jl")
include("compression.jl")
include("decompression.jl")

end # module


text:
using Documenter:
    DocMeta,
    HTML,
    MathJax3,
    deploydocs,
    makedocs
using PlutoStaticHTML

const NOTEBOOK_DIR = joinpath(pkgdir(PlutoStaticHTML), "docs", "src")
const NOTEBOOK_PATH = joinpath(NOTEBOOK_DIR, "notebook.jl")

"""
    write_notebook()

Write Pluto output to a HTML file.
This avoidings running via the Documenter.jl evaluation, which appears to just hang.
Probably similar cause as https://github.com/JuliaDocs/Documenter.jl/issues/1514.
"""
function write_notebook()
    @info "Running notebook at $NOTEBOOK_PATH"
    append_build_context = true
    html = notebook2html(NOTEBOOK_PATH; append_build_context)
    md_path = joinpath(NOTEBOOK_DIR, "notebook.md")
    md = """
        ```@eval
        # Auto generated file. Do not modify.
        ```

        # Example notebook

        The Julia version and package information below is added after running the notebook by setting `append_build_context=true`.

        ```@raw html
        $html
        ```
        """
    write(md_path, md)
    return nothing
end

write_notebook()

sitename = "PlutoStaticHTML.jl"
pages = [
    "PlutoStaticHTML" => "index.md",
    "Example notebook" => "notebook.md"
]

# Using MathJax3 since Pluto uses that engine too.
mathengine = MathJax3()
prettyurls = get(ENV, "CI", nothing) == "true"
format = HTML(; mathengine, prettyurls)
modules = [PlutoStaticHTML]
strict = true
checkdocs = :none
makedocs(; sitename, pages, format, modules, strict, checkdocs)

repo = "github.com/rikhuijzer/PlutoStaticHTML.jl.git"
push_preview = false
devbranch = "main"
deploydocs(; devbranch, repo, push_preview)

ext: jl


text:
module MatDeforElastIsoModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.DeforModelRedModule: DeforModelRed, DeforModelRed3D, DeforModelRed2DStrain, DeforModelRed2DStress, DeforModelRed2DAxisymm, DeforModelRed1D, nstressstrain, nthermstrain
import FinEtools.MatDeforModule: MatDefor, stress6vto3x3t!, stress3vto2x2t!, stress4vto3x3t!, stress4vto3x3t!
import LinearAlgebra: Transpose, Diagonal, mul!
At_mul_B!(C, A, B) = mul!(C, Transpose(A), B)
A_mul_B!(C, A, B) = mul!(C, A, B)
import LinearAlgebra: eigen, eigvals, norm, cholesky, cross, dot

const mI = Matrix(Diagonal([1.0, 1.0, 1.0, 0.5, 0.5, 0.5]))
const m1 = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0];

"""
    MatDeforElastIso

Linear isotropic elasticity  material.

```
tangentmoduli!::MTAN
update!::MUPD
thermalstrain!::MTHS
```
"""
struct  MatDeforElastIso{MR<:DeforModelRed, MTAN<:Function, MUPD<:Function, MTHS<:Function} <: MatDefor
    mr::Type{MR} # model reduction type
    mass_density::FFlt # mass density
    E::FFlt # Young's modulus
    nu::FFlt # Poisson ratio
    CTE::FFlt # Coefficient of Thermal Expansion
    D::Array{FFlt, 2} # cached matrix of 3D tangent moduli
    tangentmoduli!::MTAN
    update!::MUPD
    thermalstrain!::MTHS
end

function threedD(E::FFlt, nu::FFlt)
    lambda = E * nu / (1 + nu) / (1 - 2*(nu));
    mu = E / 2. / (1+nu);
    D = lambda * m1 * m1' + 2. * mu * mI;
    return D
end

function MatDeforElastIso(mr::Type{DeforModelRed3D},
    mass_density::FFlt, E::FFlt, nu::FFlt, CTE::FFlt)
    return MatDeforElastIso(mr, mass_density, E, nu, CTE, threedD(E, nu),
    tangentmoduli3d!, update3d!, thermalstrain3d!)
end

function MatDeforElastIso(mr::Type{DeforModelRed2DStress},
    mass_density::FFlt, E::FFlt, nu::FFlt, CTE::FFlt)
    return MatDeforElastIso(mr, mass_density, E, nu, CTE, threedD(E, nu),
    tangentmoduli2dstrs!, update2dstrs!, thermalstrain2dstrs!)
end

function MatDeforElastIso(mr::Type{DeforModelRed2DStrain},
    mass_density::FFlt, E::FFlt, nu::FFlt, CTE::FFlt)
    return MatDeforElastIso(mr, mass_density, E, nu, CTE, threedD(E, nu),
    tangentmoduli2dstrn!, update2dstrn!, thermalstrain2dstrn!)
end

function MatDeforElastIso(mr::Type{DeforModelRed2DAxisymm},
    mass_density::FFlt, E::FFlt, nu::FFlt, CTE::FFlt)
    return MatDeforElastIso(mr, mass_density, E, nu, CTE, threedD(E, nu),
    tangentmoduli2daxi!, update2daxi!, thermalstrain2daxi!)
end

function MatDeforElastIso(mr::Type{DeforModelRed1D},
    mass_density::FFlt, E::FFlt, nu::FFlt, CTE::FFlt)
    return MatDeforElastIso(mr, mass_density, E, nu, CTE, threedD(E, nu),
    tangentmoduli1d!, update1d!, thermalstrain1d!)
end

function MatDeforElastIso(mr::Type{MR}, E::FFlt, nu::FFlt) where {MR}
    mass_density = 1.0
    CTE = 0.0
    return MatDeforElastIso(mr, mass_density, E, nu, CTE)
end

################################################################################
# 3-D solid model
################################################################################

"""
    tangentmoduli3d!(self::MatDeforElastIso,
      ms::MatDeforElastIsoMS, D::FFltMat,
      t::FFlt, dt::FFlt, loc::FFltMat, label::FInt)

Calculate the material stiffness matrix.

`D` = matrix of tangent moduli, 6 x 6, supplied as a buffer and overwritten.
"""
function tangentmoduli3d!(self::MatDeforElastIso,  D::FFltMat,  t::FFlt, dt::FFlt, loc::FFltMat, label::FInt)
    copyto!(D,self.D);
    return D
end

"""
    update3d!(self::MatDeforElastIso,  stress::FFltVec, output::FFltVec,
      strain::FFltVec, thstrain::FFltVec=zeros(6), t::FFlt= 0.0, dt::FFlt= 0.0,
      loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)

Update material state.

`strain` = strain vector,
`thstrain` = thermal strain vector,
`t` = current time,
`dt` = current time step,
`loc` = location of the quadrature point in global Cartesian coordinates,
`label` = label of the finite element in which the quadrature point is found.

Output:
`stress` = stress vector, allocated by the caller with a size of the number of
  stress and strain components, `nstressstrain`. The components of the stress vector are
  calculated and stored in the `stress` vector.
`output` =  array which is (if necessary) allocated  in an appropriate size, filled
  with the output quantity, and returned.
"""
function update3d!(self::MatDeforElastIso,  stress::FFltVec, output::FFltVec,  strain::FFltVec, thstrain::FFltVec=zeros(6), t::FFlt= 0.0, dt::FFlt= 0.0,  loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)
    @assert length(stress) == nstressstrain(self.mr)
    A_mul_B!(stress, self.D, strain-thstrain);
    if quantity == :nothing
        #Nothing to be copied to the output array
    elseif quantity == :cauchy || quantity == :Cauchy
        (length(output) >= 6) || (output = zeros(6)) # make sure we can store it
        copyto!(output, stress);
    elseif quantity == :pressure || quantity == :Pressure
        output[1]  =  -sum(stress[1:3])/3.
    elseif quantity == :princCauchy || quantity == :princcauchy
        t = zeros(FFlt,3,3)
        t = stress6vto3x3t!(t,stress);
        ep = eigen(t);
        (length(output) >= 3) || (output = zeros(3)) # make sure we can store it
        copyto!(output,  sort(ep.values, rev=true));
    elseif quantity==:vonMises || quantity==:vonmises || quantity==:von_mises || quantity==:vm
        s1=stress[1]; s2=stress[2]; s3=stress[3];
        s4=stress[4]; s5=stress[5]; s6=stress[6];
        (length(output) >= 1) || (output = zeros(1)) # make sure we can store it
        output[1] = sqrt(1.0/2*((s1-s2)^2+(s1-s3)^2+(s2-s3)^2+6*(s4^2+s5^2+s6^2)))
    end
    return output
end

"""
    thermalstrain3d!(self::MatDeforElastIso, thstrain::FFltMat, dT= 0.0)

Compute thermal strain from the supplied temperature increment.
"""
function thermalstrain3d!(self::MatDeforElastIso, thstrain::FFltVec, dT= 0.0)
    @assert length(thstrain) == nthermstrain(self.mr)
    thstrain[1] = self.CTE*dT
    thstrain[2] = self.CTE*dT
    thstrain[3] = self.CTE*dT
    thstrain[4] = 0.0
    thstrain[5] = 0.0
    thstrain[6] = 0.0
    return thstrain
end

################################################################################
# 2-D plane stress
################################################################################

"""
    tangentmoduli2dstrs!(self::MatDeforElastIso,
      ms::MatDeforElastIsoMS, D::FFltMat,
      t::FFlt, dt::FFlt, loc::FFltMat, label::FInt)

Calculate the material stiffness matrix.

`D` = matrix of tangent moduli, 3 x 3, supplied as a buffer and overwritten.
"""
function tangentmoduli2dstrs!(self::MatDeforElastIso, D::FFltMat,  t::FFlt, dt::FFlt, loc::FFltMat, label::FInt)
    D[1:2, 1:2] = self.D[1:2, 1:2] -  (reshape(self.D[1:2,3], 2, 1) * reshape(self.D[3,1:2], 1, 2))/self.D[3, 3]
    ix=[1, 2, 4];
    for i = 1:3
        D[3,i] = D[i,3] = self.D[4, ix[i]];
    end
    return D
end

"""
    update2dstrs!(self::MatDeforElastIso, stress::FFltVec, output::FFltVec,
      strain::FFltVec, thstrain::FFltVec=zeros(3), t::FFlt= 0.0, dt::FFlt= 0.0,
      loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)

Update material state.

`strain` = strain vector,
`thstrain` = thermal strain vector,
`t` = current time,
`dt` = current time step,
`loc` = location of the quadrature point in global Cartesian coordinates,
`label` = label of the finite element in which the quadrature point is found.

Output:
`stress` = stress vector, allocated by the caller with a size of the number of
  stress and strain components, `nstressstrain`. The components of the stress vector are
  calculated and stored in the `stress` vector.
`output` =  array which is (if necessary) allocated  in an appropriate size, filled
  with the output quantity, and returned.
"""
function update2dstrs!(self::MatDeforElastIso, stress::FFltVec, output::FFltVec,  strain::FFltVec, thstrain::FFltVec=zeros(3), t::FFlt= 0.0, dt::FFlt= 0.0,  loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)
    @assert length(stress) == nstressstrain(self.mr)
    D = zeros(3, 3)
    tangentmoduli2dstrs!(self, D, t, dt, loc, label)
    A_mul_B!(stress, D, strain-thstrain);
    if quantity == :nothing
        #Nothing to be copied to the output array
    elseif quantity == :cauchy || quantity == :Cauchy
        (length(output) >= 3) || (output = zeros(3)) # make sure we can store it
        copyto!(output, stress);
    elseif quantity == :pressure || quantity == :Pressure
        (length(output) >= 1) || (output = zeros(1)) # make sure we can store it
        output[1] = -sum(stress[1:2])/3.
    elseif quantity == :princCauchy || quantity == :princcauchy
        t = zeros(FFlt,2,2)
        t = stress3vto2x2t!(t,stress);
        ep = eigen(t);
        (length(output) >= 2) || (output = zeros(2)) # make sure we can store it
        copyto!(output,  sort(ep.values, rev=true));
    elseif quantity==:vonMises || quantity==:vonmises || quantity==:von_mises || quantity==:vm
        s1=stress[1]; s2=stress[2]; s3=0.0;
        s4=stress[3]; s5=0.0; s6=0.0;
        (length(output) >= 1) || (output = zeros(1)) # make sure we can store it
        output[1] = sqrt(1.0/2*((s1-s2)^2+(s1-s3)^2+(s2-s3)^2+6*(s4^2+s5^2+s6^2)))
    end
    return output
end

"""
    thermalstrain2dstrs!(self::MatDeforElastIso, thstrain::FFltVec, dT= 0.0)

Compute thermal strain from the supplied temperature increment.
"""
function thermalstrain2dstrs!(self::MatDeforElastIso, thstrain::FFltVec, dT= 0.0)
    @assert length(thstrain) == nthermstrain(self.mr)
    thstrain[1] = self.CTE*dT
    thstrain[2] = self.CTE*dT
    thstrain[3] = 0.0
    return thstrain
end

################################################################################
# 2-D plane strain
################################################################################

"""
    tangentmoduli2dstrn!(self::MatDeforElastIso,
      ms::MatDeforElastIsoMS, D::FFltMat,
      t::FFlt, dt::FFlt, loc::FFltMat, label::FInt)

Calculate the material stiffness matrix.

`D` = matrix of tangent moduli, 3 x 3, supplied as a buffer and overwritten.
"""
function tangentmoduli2dstrn!(self::MatDeforElastIso,  D::FFltMat,  t::FFlt, dt::FFlt, loc::FFltMat, label::FInt)
    ix = [1, 2, 4];
    for i = 1:length(ix)
        for j = 1:length(ix)
        D[j,i] = self.D[ix[j], ix[i]];
        end
    end
    return D
end

"""
    update2dstrn!(self::MatDeforElastIso,  stress::FFltVec, output::FFltVec,
      strain::FFltVec, thstrain::FFltVec=zeros(3), t::FFlt= 0.0, dt::FFlt= 0.0,
      loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)

Update material state.

`strain` = strain vector,
`thstrain` = thermal strain vector,
`t` = current time,
`dt` = current time step,
`loc` = location of the quadrature point in global Cartesian coordinates,
`label` = label of the finite element in which the quadrature point is found.

Output:
`stress` = stress vector, allocated by the caller with a size of the number of
  stress and strain components, `nstressstrain`. The components of the stress vector are
  calculated and stored in the `stress` vector.
`output` =  array which is (if necessary) allocated  in an appropriate size, filled
  with the output quantity, and returned.

Note on the principal stresses: The principal stresses are calculated for the
fully three-dimensional  stress state, that is not the "in-plane" maximum and
minimum,  but rather  the  three-dimensional maximum (1) and minimum (3).
The intermediate principal stress is (2).
"""
function update2dstrn!(self::MatDeforElastIso,  stress::FFltVec, output::FFltVec,  strain::FFltVec, thstrain::FFltVec=zeros(4), t::FFlt= 0.0, dt::FFlt= 0.0,  loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)
    @assert length(stress) == nstressstrain(self.mr)
    D = zeros(3, 3)
    tangentmoduli2dstrn!(self, D, t, dt, loc, label)
    A_mul_B!(stress, D, strain-thstrain[1:3]);
    if quantity == :nothing
        #Nothing to be copied to the output array
    elseif quantity == :cauchy || quantity == :Cauchy
        # sigmax, sigmay, tauxy, sigmaz
        # thstrain[4] =The through the thickness thermal strain
        sz = dot(self.D[3, 1:2], strain[1:2]-thstrain[1:2])-self.D[3,3]*thstrain[4];
        (length(output) >= 4) || (output = zeros(4)) # make sure we can store it
        copyto!(output, stress); output[4] = sz
    elseif quantity == :pressure || quantity == :Pressure
        (length(output) >= 1) || (output = zeros(1)) # make sure we can store it
        sz = dot(self.D[3, 1:2], strain[1:2]-thstrain[1:2])-self.D[3,3]*thstrain[4];
        output[1] = -(sum(stress[[1,2]]) + sz)/3.
    elseif quantity == :princCauchy || quantity == :princcauchy
        (length(output) >= 3) || (output = zeros(3)) # make sure we can store it
        t = zeros(FFlt, 3,3)
        sz = dot(self.D[3, 1:2], strain[1:2]-thstrain[1:2])-self.D[3,3]*thstrain[4];
        t = stress4vto3x3t!(t, vcat(stress[1:3], [sz]));
        ep = eigen(t);
        (length(output) >= 3) || (output = zeros(3)) # make sure we can store it
        copyto!(output,  sort(ep.values, rev=true));
    elseif quantity==:vonMises || quantity==:vonmises || quantity==:von_mises || quantity==:vm
        (length(output) >= 1) || (output = zeros(1)) # make sure we can store it
        sz = dot(self.D[3, 1:2], strain[1:2]-thstrain[1:2])-self.D[3,3]*thstrain[4];
        s1=stress[1]; s2=stress[2]; s3=sz;
        s4=stress[3]; s5=0.0; s6=0.0;
        (length(output) >= 1) || (output = zeros(1)) # make sure we can store it
        output[1] = sqrt(1.0/2*((s1-s2)^2+(s1-s3)^2+(s2-s3)^2+6*(s4^2+s5^2+s6^2)))
    end
    return output
end

"""
    thermalstrain2dstrn!(self::MatDeforElastIso, thstrain::FFltVec, dT= 0.0)

Compute thermal strain from the supplied temperature increment.

The thermal strain is evaluated  for the  two normal strains, the shear
strain, and finally for the through the thickness strain.
thstrain = [ex, ey, 0.0, ez].
"""
function thermalstrain2dstrn!(self::MatDeforElastIso, thstrain::FFltVec, dT= 0.0)
    @assert length(thstrain) == nthermstrain(self.mr)
    thstrain[1] = self.CTE*dT
    thstrain[2] = self.CTE*dT
    thstrain[3] = 0.0
    thstrain[4] = self.CTE*dT
    return thstrain
end

################################################################################
# 2-D axially symmetric
################################################################################

"""
    tangentmoduli2daxi!(self::MatDeforElastIso,
      ms::MatDeforElastIsoMS, D::FFltMat,
      t::FFlt, dt::FFlt, loc::FFltMat, label::FInt)

Calculate the material stiffness matrix.

`D` = matrix of tangent moduli, 3 x 3, supplied as a buffer and overwritten.
"""
function tangentmoduli2daxi!(self::MatDeforElastIso,  D::FFltMat,  t::FFlt, dt::FFlt, loc::FFltMat, label::FInt)
    for i = 1:4
        for j = 1:4
        D[j,i] = self.D[i, j];
        end
    end
    return D
end

"""
    update2daxi!(self::MatDeforElastIso,  stress::FFltVec, output::FFltVec,
      strain::FFltVec, thstrain::FFltVec=zeros(3), t::FFlt= 0.0, dt::FFlt= 0.0,
      loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)

Update material state.

`strain` = strain vector,
`thstrain` = thermal strain vector,
`t` = current time,
`dt` = current time step,
`loc` = location of the quadrature point in global Cartesian coordinates,
`label` = label of the finite element in which the quadrature point is found.

Output:
`stress` = stress vector, allocated by the caller with a size of the number of
  stress and strain components, `nstressstrain`. The components of the stress vector are
  calculated and stored in the `stress` vector.
`output` =  array which is (if necessary) allocated  in an appropriate size, filled
  with the output quantity, and returned.
"""
function update2daxi!(self::MatDeforElastIso,  stress::FFltVec, output::FFltVec,  strain::FFltVec, thstrain::FFltVec=zeros(3), t::FFlt= 0.0, dt::FFlt= 0.0,  loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)
    @assert length(stress) == nstressstrain(self.mr)
    D = zeros(4, 4)
    tangentmoduli2daxi!(self, D, t, dt, loc, label)
    A_mul_B!(stress, D, strain-thstrain);
    if quantity == :nothing
        #Nothing to be copied to the output array
    elseif quantity == :cauchy || quantity == :Cauchy
        (length(output) >= 4) || (output = zeros(4)) # make sure we can store it
        copyto!(output, stress);
    elseif quantity == :pressure || quantity == :Pressure
        (length(output) >= 1) || (output = zeros(1)) # make sure we can store it
        output[1] = -sum(stress[[1,2,3]])/3.
    elseif quantity == :princCauchy || quantity == :princcauchy
        t = zeros(FFlt,3,3)
        t = stress4vto3x3t!(t, stress[[1,2,4,3]]);
        ep = eigen(t);
        (length(output) >= 3) || (output = zeros(3)) # make sure we can store it
        copyto!(output,  sort(ep.values, rev=true));
    elseif quantity==:vonMises || quantity==:vonmises || quantity==:von_mises || quantity==:vm
        s1=stress[1]; s2=stress[2]; s3=stress[3];
        s4=stress[4]; s5=0.0; s6=0.0;
        (length(output) >= 1) || (output = zeros(1)) # make sure we can store it
        output[1] = sqrt(1.0/2*((s1-s2)^2+(s1-s3)^2+(s2-s3)^2+6*(s4^2+s5^2+s6^2)))
    end
    return output
end

"""
    thermalstrain2daxi!(self::MatDeforElastIso, thstrain::FFltVec, dT= 0.0)

Compute thermal strain from the supplied temperature increment.

The thermal strain is evaluated  for the  three normal strains and the shear
strain.
"""
function thermalstrain2daxi!(self::MatDeforElastIso, thstrain::FFltVec, dT= 0.0)
    @assert length(thstrain) == nthermstrain(self.mr)
    thstrain[1] = self.CTE*dT
    thstrain[2] = self.CTE*dT
    thstrain[3] = self.CTE*dT
    thstrain[4] = 0.0
    return thstrain
end

################################################################################
# 2-D axially symmetric
################################################################################

"""
    tangentmoduli2daxi!(self::MatDeforElastIso,
      ms::MatDeforElastIsoMS, D::FFltMat,
      t::FFlt, dt::FFlt, loc::FFltMat, label::FInt)

Calculate the material stiffness matrix.

`D` = matrix of tangent moduli, 3 x 3, supplied as a buffer and overwritten.
"""
function tangentmoduli1d!(self::MatDeforElastIso,
    D::FFltMat,
    t::FFlt, dt::FFlt, loc::FFltMat, label::FInt)
    D[1, 1] = self.E;
    return D
end

"""
    update1d!(self::MatDeforElastIso,  stress::FFltVec, output::FFltVec,
      strain::FFltVec, thstrain::FFltVec=zeros(1), t::FFlt= 0.0, dt::FFlt= 0.0,
      loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)

Update material state.

`strain` = strain vector,
`thstrain` = thermal strain vector,
`t` = current time,
`dt` = current time step,
`loc` = location of the quadrature point in global Cartesian coordinates,
`label` = label of the finite element in which the quadrature point is found.

Output:
`stress` = stress vector, allocated by the caller with a size of the number of
  stress and strain components, `nstressstrain`. The components of the stress vector are
  calculated and stored in the `stress` vector.
`output` =  array which is (if necessary) allocated  in an appropriate size, filled
  with the output quantity, and returned.
"""
function update1d!(self::MatDeforElastIso,  stress::FFltVec, output::FFltVec,  strain::FFltVec, thstrain::FFltVec=zeros(1), t::FFlt= 0.0, dt::FFlt= 0.0,  loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)
    @assert length(stress) == nstressstrain(self.mr)
    D = zeros(1, 1)
    tangentmoduli1d!(self, D, t, dt, loc, label)
    A_mul_B!(stress, D, strain-thstrain);
    if quantity == :nothing
        #Nothing to be copied to the output array
    elseif quantity == :cauchy || quantity == :Cauchy
        (length(output) >= 1) || (output = zeros(1)) # make sure we can store it
        copyto!(output, stress);
    elseif quantity == :pressure || quantity == :Pressure
        (length(output) >= 1) || (output = zeros(1)) # make sure we can store it
        output[1] = -sum(stress[[1]])/3.
    elseif quantity == :princCauchy || quantity == :princcauchy
        copyto!(output,  stress[1]);
    elseif quantity==:vonMises || quantity==:vonmises || quantity==:von_mises || quantity==:vm
        s1=stress[1]; s2=0.0; s3=0.0;
        s4=00.0; s5=0.0; s6=0.0;
        (length(output) >= 1) || (output = zeros(1)) # make sure we can store it
        output[1] = sqrt(1.0/2*((s1-s2)^2+(s1-s3)^2+(s2-s3)^2+6*(s4^2+s5^2+s6^2)))
    end
    return output
end

"""
    thermalstrain1d!(self::MatDeforElastIso, thstrain::FFltVec, dT= 0.0)

Compute thermal strain from the supplied temperature increment.

The thermal strain is evaluated  for the axial strain.
"""
function thermalstrain1d!(self::MatDeforElastIso, thstrain::FFltVec, dT= 0.0)
    @assert length(thstrain) == nthermstrain(self.mr)
    thstrain[1] = self.CTE*dT
    return thstrain
end

end


text:
using LinearAlgebra
#import Gaugefields:Verbose_level,Verbose_3,Verbose_2,Verbose_1,println_verbose_level3
using InteractiveUtils
import Gaugefields.Verboseprint_mpi:Verbose_print,println_verbose_level1,println_verbose_level2,println_verbose_level3

#export bicg,bicgstag,shiftedcg,bicgstab_evenodd,reducedshiftedcg,cg

function add!(b,Y,a,X) #b*Y + a*X -> Y
    LinearAlgebra.axpby!(a,X,b,Y) #X*a + Y*b -> Y
end

function add!(Y,a,X) #Y + a*X -> Y
    LinearAlgebra.axpby!(a,X,1,Y) #X*a + Y -> Y
end

function add!(b,Y,a,X,iseven::Bool) #b*Y + a*X -> Y
    LinearAlgebra.axpby!(a,X,b,Y,iseven) #X*a + Y*b -> Y
end

function add!(Y,a,X,iseven::Bool) #Y + a*X -> Y
    LinearAlgebra.axpby!(a,X,1,Y,iseven) #X*a + Y -> Y
end

function bicg(x,A,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) #Ax=b
    #println_verbose_level3()
    println_verbose_level3(verbose,"--------------------------------------")
    println_verbose_level3(verbose,"bicg method")
    #=
    res = deepcopy(b)
    temp1 = similar(x)
    p = similar(x)
    q = similar(x)
    s = similar(x)

    mul!(temp1,A,x)
    add!(res,-1,temp1)
    =#

    temps = get_temporaryvectors_forCG(A)
    res = temps[1]
    substitute_fermion!(res,b)
    temp1 = temps[2]
    mul!(temp1,A,x)
    add!(res,-1,temp1)
    p = temps[4]
    q = temps[5]
    s = temps[6]

    rnorm = real(res⋅res)
    if rnorm < eps
        return
    end
    #println(rnorm)

    mul!(p,A',res)
    c1 = p ⋅ p

    for i=1:maxsteps
        mul!(q,A,p)
        #! ...  c2 = < q | q >
        c2 = q ⋅ q

        alpha = c1 / c2
        #! ...  x   = x   + alpha * p
        add!(x,alpha,p)
        #...  res = res - alpha * q
        add!(res,-alpha,q)
        rnorm = real(res ⋅ res)
        println_verbose_level3(verbose,"$i-th eps: $rnorm")

        if rnorm < eps
            #println("Converged at $i-th step. eps: $rnorm")
            println_verbose_level3(verbose,"Converged at $i-th step. eps: $rnorm")
            println_verbose_level3(verbose,"--------------------------------------")
            return
        end

        mul!(s,A',res)

        #c3 = s * s
        c3 = s ⋅ s

        beta = c3 / c1
        c1 = c3

        add!(beta,p,1,s) #p = beta*p + s

    end

    error("""
    The BICG is not converged! with maxsteps = $(maxsteps)
    residual is $rnorm
    maxsteps should be larger.""")

end

function bicgstab(x,A,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) #Ax=b
    println_verbose_level3(verbose,"--------------------------------------")
    println_verbose_level3(verbose,"bicg-stab method")

    #=
    r = deepcopy(b)
    temp1 = similar(x)
    mul!(temp1,A,x)
    add!(r,-1,temp1)

    rs = deepcopy(r)
    p = deepcopy(r)
    Ap = similar(r)
    s = similar(r)
    t = similar(r)
    =#

    temps = get_temporaryvectors_forCG(A)
    r = temps[1]
    substitute_fermion!(r,b)
    temp1 = temps[2]
    mul!(temp1,A,x)
    add!(r,-1,temp1)

    rs = temps[3]
    substitute_fermion!(rs,r)
    p = temps[4]
    substitute_fermion!(p,r)
    Ap = temps[5]
    s = temps[6]
    t = temps[7]

    rnorm = real(r⋅r)

    if rnorm < eps
        return
    end

    for i=1:maxsteps
        c1 = dot(rs,r)
        mul!(Ap,A,p)
        c2 = dot(rs,Ap)
        α = c1/c2
        #s = r - α*A*p
        add!(0,s,1,r)
        add!(s,-α,Ap)
        mul!(t,A,s)
        d1 = dot(t,s)
        d2 = dot(t,t)
        ω = d1/d2

        #r = (1-ω A)s
        add!(0,r,1,s)
        add!(r,-ω,t)

        #x = x + ωs+ αp
        add!(x,ω,s)
        add!(x,α,p)

        β = (dot(rs,r)/c1)*(α/ω)

        #p = r + β*(1-ωA)*p
        add!(β,p,1,r)
        add!(p,-ω*β,Ap)

        rnorm = real(r ⋅ r)
        println_verbose_level3(verbose,"$i-th eps: $rnorm")

        if rnorm < eps
            println_verbose_level3(verbose,"Converged at $i-th step. eps: $rnorm")
            println_verbose_level3(verbose,"--------------------------------------")
            return
        end

    end

    error("""
    The BICGstab is not converged! with maxsteps = $(maxsteps)
    residual is $rnorm
    maxsteps should be larger.""")

end

function bicgstab_3(x,A,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) #Ax=b
    println_verbose_level3(verbose,"--------------------------------------")
    println_verbose_level3(verbose,"bicg-stab method")
    r = deepcopy(b)
    temp1 = similar(x)
    mul!(temp1,A,x)
    add!(r,-1,temp1)

    rnorm = real(r⋅r)

    if rnorm < eps
        return
    end

    #function test()

    rs = deepcopy(r)
    #rs = deepcopy(b)
    p = deepcopy(r)
    s = similar(r)
    t = similar(r)
    α = 0
    Ap = similar(p)
    As = similar(p)
    ω = 0

    for i=0:maxsteps
        mul!(Ap,A,p)
        c1 = dot(r,rs)
        α = c1/dot(Ap,rs)
        #s = r - alpha*v
        add!(0,s,1,r)
        add!(s,-α,Ap)

        mul!(As,A,s)
        ω = dot(As,s)/dot(As,As)
        add!(x,α,p)
        add!(x,ω,s)

        add!(0,r,1,s) # r = a*r + b*s (a = 0,b = 1)
        add!(r,-ω,t)

        β = (dot(r,rs)/c1)*α/ω

        #p = r + β*(p - ωm*vm)
        add!(β,p,1,r)
        add!(p,-ω*β,Ap)

        rnorm = real(r ⋅ r)
        println_verbose_level3(verbose,"$i-th eps: $rnorm")

        if rnorm < eps
            println_verbose_level3(verbose,"Converged at $i-th step. eps: $rnorm")
            println_verbose_level3(verbose,"--------------------------------------")
            return
        end
    end
    #end

    #= other
    rs = deepcopy(r)
    ρm = 1
    α = 1
    ωm = 1
    vm = similar(x)
    p = similar(x)
    s = similar(x)
    t = similar(x)

    println("rr ",dot(r,r))
    println("rsrs ",dot(rs,rs))
    println("rsr ",dot(rs,r))

    for i=1:maxsteps
        ρ = dot(rs,r)
        println(ρ    )
        β =(ρ/ρm)*(α/ωm)
        #p = r + β*(p - ωm*vm)
        add!(β,p,1,r)
        add!(p,-ωm*β,vm)

        mul!(vm,A,p)
        α = ρ/dot(rs,vm)

        #s = r - alpha*v
        add!(0,s,1,r)
        add!(s,-α,vm)

        mul!(t,A,s)
        ωm = dot(t,s)/dot(r,r)

        add!(x,α,p)
        add!(x,ωm,s)

        add!(0,r,1,s) # r = a*r + b*s (a = 0,b = 1)
        add!(r,-ωm,t)

        rnorm = real(r ⋅ r)
        println_verbose_level3(verbose,"$i-th eps: $rnorm")

        if rnorm < eps
            println_verbose_level3(verbose,"Converged at $i-th step. eps: $rnorm")
            println_verbose_level3(verbose,"--------------------------------------")
            return
        end
    end

    =#

    #= text book
    p = deepcopy(r)
    s = similar(r)
    v = similar(x)
    t = similar(x)
    β = 0.0im
    ω = 0.0im
    ρ2 = 0.0im
    α = 0.0im

    for i=0:maxsteps
        ρ = dot(rs,r)
        if i != 0
            β = α*ρ/(ω*ρ2)

            add!(β,p,1,r)
            add!(p,-ω*β,v) #p = r + beta*(p -omega*v)
        end
        mul!(v,A,p)
        c2 = dot(rs,v)
        α = ρ/c2

        #s = r - α*v
        add!(0,s,1,r)
        add!(s,-α,v)

        rnorm = real(dot(s,s))
        if rnorm < eps
            add!(x,α,p)
            println_verbose_level3(verbose,"Converged at $i-th step. eps: $rnorm")
            println_verbose_level3(verbose,"--------------------------------------")

            return
        end
        mul!(t,A,s)

        ω = dot(t,s)/dot(t,t)
        #r = s - ω*t
        add!(0,r,1,s) # r = a*r + b*s (a = 0,b = 1)
        add!(s,-ω,t)

        add!(x,α,p)
        add!(x,ω,s)
        ρ2 = ρ

    end

    =#
    error("""
    The BICGstab is not converged! with maxsteps = $(maxsteps)
    residual is $rnorm
    maxsteps should be larger.""")

end

function bicgstab_YN(x,A,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) #Ax=b
    println_verbose_level3(verbose,"--------------------------------------")
    println_verbose_level3(verbose,"bicg-stab method")
    r = deepcopy(b)
    temp1 = similar(x)
    mul!(temp1,A,x)
    add!(r,-1,temp1)

    rnorm = real(r⋅r)

    if rnorm < eps
        return
    end
    rs = deepcopy(r)
    p = deepcopy(r)
    s = similar(r)
    Ap = similar(x)
    As = similar(x)
    β = 0.0im
    ω = 0.0im

    for i=0:maxsteps
        if i != 0
            add!(β,p,1,r)
            add!(p,-ω*β,Ap)
        end
        c1 = dot(rs,r)
        mul!(Ap,A,p)
        c2 = dot(rs,Ap)
        α = c1/c2
        add!(0,s,1,r)
        add!(s,-α,Ap)
        mul!(As,A,s)
        s1 = dot(As,s)
        s2 = dot(As,As)

        ω = s1/s2
        add!(x,α,p)
        add!(x,ω,s)
        add!(0,r,1,s)
        add!(r,-ω,As)

        rnorm = real(r ⋅ r)
        println_verbose_level3(verbose,"$i-th eps: $rnorm")

        if rnorm < eps
            println_verbose_level3(verbose,"Converged at $i-th step. eps: $rnorm")
            println_verbose_level3(verbose,"--------------------------------------")
            return
        end

        β = (α/ω)*dot(rs,r)/c1

    end

    error("""
    The BICG is not converged! with maxsteps = $(maxsteps)
    residual is $rnorm
    maxsteps should be larger.""")

end

function bicgstab_evenodd(x,A,b,iseven;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) #Ax=b
    println_verbose_level3(verbose,"--------------------------------------")
    println_verbose_level3(verbose,"bicg-stab even-odd method")
    r = deepcopy(b)
    temp1 = similar(x)
    mul!(temp1,A,x)
    #println(dot(temp1,temp1,iseven))
    add!(r,-1,temp1,iseven)

    rnorm = real(dot(r,r,iseven))

    if rnorm < eps
        return
    end

    rs = deepcopy(r)
    p = deepcopy(r)
    Ap = similar(r)
    s = similar(r)
    t = similar(r)

    for i=1:maxsteps
        c1 = dot(rs,r,iseven)
        mul!(Ap,A,p)
        c2 = dot(rs,Ap,iseven)
        α = c1/c2
        #s = r - α*A*p
        add!(0,s,1,r,iseven)
        add!(s,-α,Ap,iseven)
        mul!(t,A,s)
        d1 = dot(t,s,iseven)
        d2 = dot(t,t,iseven)
        ω = d1/d2

        #r = (1-ω A)s
        add!(0,r,1,s,iseven)
        add!(r,-ω,t,iseven)

        #x = x + ωs+ αp
        add!(x,ω,s,iseven)
        add!(x,α,p,iseven)

        β = (dot(rs,r,iseven)/c1)*(α/ω)

        #p = r + β*(1-ωA)*p
        add!(β,p,1,r,iseven)
        add!(p,-ω*β,Ap,iseven)

        rnorm = real(dot(r,r,iseven))
        println_verbose_level3(verbose,"$i-th eps: $rnorm")

        if rnorm < eps
            println_verbose_level3(verbose,"Converged at $i-th step. eps: $rnorm")
            println_verbose_level3(verbose,"--------------------------------------")
            return
        end

    end

    error("""
    The BICGstab is not converged! with maxsteps = $(maxsteps)
    residual is $rnorm
    maxsteps should be larger.""")

end

function bicgstab_evenodd_YN(x,A,b,iseven;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) #Ax=b
    println_verbose_level3(verbose,"--------------------------------------")
    println_verbose_level3(verbose,"bicg-stab even-odd method")

    r = deepcopy(b)
    temp1 = similar(x)
    mul!(temp1,A,x)
    println(dot(temp1,temp1,iseven))
    add!(r,-1,temp1,iseven)
    rs = deepcopy(r)
    p = deepcopy(r)
    s = similar(r)
    Ap = similar(x)
    As = similar(x)
    β = 0.0im
    ω = 0.0im

    rnorm = real(dot(r,r,iseven))
    if rnorm < eps
        return
    end

    for i=0:maxsteps

        #=
        if i != 0
            add!(β,p,1,r,iseven)
            add!(p,-ω*β,Ap,iseven)
        end
        =#
        c1 = dot(rs,r,iseven)
        mul!(Ap,A,p)
        c2 = dot(rs,Ap,iseven)
        α = c1/c2
        add!(0,s,1,r,iseven)
        add!(s,-α,Ap,iseven)
        mul!(As,A,s)
        s1 = dot(As,s,iseven)
        s2 = dot(As,As,iseven)

        ω = s1/s2
        add!(x,α,p,iseven)
        add!(x,ω,s,iseven)
        add!(0,r,1,s,iseven)
        add!(r,-ω,As,iseven)

        rnorm = real(dot(r,r,iseven))
        println_verbose_level3(verbose,"$i-th eps: $rnorm")

        if rnorm < eps
            println_verbose_level3(verbose,"Converged at $i-th step. eps: $rnorm")
            println_verbose_level3(verbose,"--------------------------------------")
            return
        end

        β = (α/ω)*dot(rs,r,iseven)/c1

        #if i != 0
            add!(β,p,1,r,iseven)
            add!(p,-ω*β,Ap,iseven)
        #end
    end

    error("""
    The BICG is not converged! with maxsteps = $(maxsteps)
    residual is $rnorm
    maxsteps should be larger.""")

end

function cg(x,A,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) #Ax=b
    temps = get_temporaryvectors_forCG(A)

    println_verbose_level3(verbose,"--------------------------------------")
    println_verbose_level3(verbose,"cg method")

    #=
    res = deepcopy(b)
    temp1 = similar(x)
    mul!(temp1,A,x)
    add!(res,-1,temp1)
    q = similar(x)
    p = deepcopy(res)
    =#

    res = temps[1]
    substitute_fermion!(res,b)
    temp1 = temps[2]
    mul!(temp1,A,x)
    add!(res,-1,temp1)
    q = temps[3]
    p = temps[4]
    substitute_fermion!(p,res)

    #p = deepcopy(res)

    rnorm = real(res⋅res)
    #println(rnorm)

    if rnorm < eps
        return
    end

    c1 = p ⋅ p

    for i=1:maxsteps
        mul!(q,A,p)
        c2 = dot(p,q)
        #c2 = p ⋅ q

        α = c1 / c2
        #! ...  x   = x   + alpha * p
        #println("add2")
        add!(x,α,p)
        #...  res = res - alpha * q
        #println("add1")
        add!(res,-α,q)
        c3 = res ⋅ res
        rnorm = real(c3)
        println_verbose_level3(verbose,"$i-th eps: $rnorm")

        if rnorm < eps
            #println("$i eps: $eps rnorm $rnorm")
            println_verbose_level3(verbose,"Converged at $i-th step. eps: $rnorm")
            println_verbose_level3(verbose,"--------------------------------------")
            return
        end

        β = c3 / c1
        c1 = c3

        #println("add3")
        add!(β,p,1,res) #p = beta*p + s

    end

    error("""
    The CG is not converged! with maxsteps = $(maxsteps)
    residual is $rnorm
    maxsteps should be larger.""")

end

function Base.:*(x::Array{T,1},y::Array{T,1}) where T <: Number
    return x'*y
end

function shiftedcg(vec_x,vec_β,x,A,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) #Ax=b

    println_verbose_level3(verbose,"--------------------------------------")
    println_verbose_level3(verbose,"shifted cg method")
    N = length(vec_β)
    temp1 = similar(b)
    r = deepcopy(b)
    p = deepcopy(b)
    q = similar(b)

    vec_r = Array{typeof(r),1}(undef,N)
    vec_p = Array{typeof(p),1}(undef,N)
    for j=1:N
        vec_r[j] = deepcopy(b)
        vec_p[j] = deepcopy(b)
    end

    αm = 1.0
    βm = 0.0

    ρm = ones(ComplexF64,N)
    ρ0 = ones(ComplexF64,N)
    ρp = ones(ComplexF64,N)
    residual = 0

    for i=1:maxsteps
        mul!(q,A,p)

        pAp = p ⋅ q

        rr = dot(r,r)
        αk = rr / pAp

        #! ...  x   = x   + alpha * p
        add!(x,αk,p)

        #...  r = r - alpha * q
        add!(r,-αk,q)

        βk = dot(r,r)/ rr
        add!(βk,p,1,r) #p = beta*p + r

        for j=1:N
            ρkj = ρ0[j]
            if abs(ρkj) < eps
                continue
            end
            ρkmj =ρm[j]
            ρp[j] = ρkj*ρkmj*αm/(ρkmj*αm*(1.0+αk*vec_β[j])+αk*βm*(ρkmj-ρkj))
            αkj = (ρp[j]/ρkj)*αk

            add!(vec_x[j],αkj,vec_p[j])
            βkj = (ρp[j]/ρkj)^2*βk
            add!(βkj,vec_p[j],ρp[j],r)

        end

        ρm[:] = ρ0[:]
        ρ0[:] = ρp[:]
        αm = αk
        βm = βk

        ρMAX = maximum(abs.(ρp))^2
        residual = abs(rr*ρMAX)
        println_verbose_level3(verbose,"$i-th eps: $residual")

        if abs(residual) < eps
            println_verbose_level3(verbose,"Converged at $i-th step. eps: $residual")
            println_verbose_level3(verbose,"--------------------------------------")
            return
        end

    end

    error("""
    The shifted CG is not converged! with maxsteps = $(maxsteps)
    residual is $residual
    maxsteps should be larger.""")

end

function shiftedbicg(vec_x,vec_β,x,A,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) #Ax=b

    println_verbose_level3(verbose,"--------------------------------------")
    println_verbose_level3(verbose,"shifted cg method")
    N = length(vec_β)
    temp1 = similar(b)
    r = deepcopy(b)
    p = deepcopy(b)
    q = similar(b)
    s = similar(b)

    vec_r = Array{typeof(r),1}(undef,N)
    vec_p = Array{typeof(p),1}(undef,N)
    for j=1:N
        vec_r[j] = deepcopy(b)
        vec_p[j] = deepcopy(b)
    end

    αm = 1.0
    βm = 0.0

    ρm = ones(ComplexF64,N)
    ρ0 = ones(ComplexF64,N)
    ρp = ones(ComplexF64,N)
    residual = 0

    mul!(p,A',r)
    c1 = p ⋅ p

    for i=1:maxsteps
        mul!(q,A,p)
        c2 = q ⋅ q
        alpha = c1 / c2
        rr = dot(r,r)
        #! ...  x   = x   + alpha * p
        add!(x,alpha,p)
        #...  res = res - alpha * q
        add!(r,-alpha,q)

        mul!(s,A',r)

        #c3 = s * s
        c3 = s ⋅ s

        beta = c3 / c1
        c1 = c3

        add!(beta,p,1,s) #p = beta*p + s

        αk = alpha
        βk = beta

        for j=1:N
            ρkj = ρ0[j]
            if abs(ρkj) < eps
                continue
            end
            ρkmj =ρm[j]
            ρp[j] = ρkj*ρkmj*αm/(ρkmj*αm*(1.0+αk*vec_β[j])+αk*βm*(ρkmj-ρkj))
            αkj = (ρp[j]/ρkj)*αk

            add!(vec_x[j],αkj,vec_p[j])
            βkj = (ρp[j]/ρkj)^2*βk
            #add!(βkj,vec_p[j],ρp[j],r)
            add!(βkj,vec_p[j],ρp[j],s)

        end

        ρm[:] = ρ0[:]
        ρ0[:] = ρp[:]
        αm = αk
        βm = βk

        ρMAX = maximum(abs.(ρp))^2
        residual = abs(rr*ρMAX)
        println_verbose_level3(verbose,"$i-th eps: $residual")

        if abs(residual) < eps
            println_verbose_level3(verbose,"Converged at $i-th step. eps: $residual")
            println_verbose_level3(verbose,"--------------------------------------")
            return
        end

    end

    error("""
    The shifted CG is not converged! with maxsteps = $(maxsteps)
    residual is $residual
    maxsteps should be larger.""")

end

function reducedshiftedcg(leftvec,vec_β,x,A,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) #Ax=b
    println_verbose_level3(verbose,"--------------------------------------")
    println_verbose_level3(verbose,"shifted cg method")
    N = length(vec_β)
    temp1 = similar(b)
    r = deepcopy(b)
    p = deepcopy(b)
    q = similar(b)

    #=
    vec_r = Array{typeof(r),1}(undef,N)
    vec_p = Array{typeof(p),1}(undef,N)
    for j=1:N
        vec_r[j] = deepcopy(b)
        vec_p[j] = deepcopy(b)
    end
    =#

    Σ = leftvec ⋅ b

    θ = zeros(ComplexF64,N)
    Π = ones(ComplexF64,N) .* Σ

    αm = 1.0
    βm = 0.0

    ρm = ones(ComplexF64,N)
    ρ0 = ones(ComplexF64,N)
    ρp = ones(ComplexF64,N)
    residual = 0

    for i=1:maxsteps
        mul!(q,A,p)

        pAp = p ⋅ q
        rr = r*r
        αk = rr / pAp

        #! ...  x   = x   + alpha * p
        add!(x,αk,p)

        #...  r = r - alpha * q
        add!(r,-αk,q)

        βk = (r*r)/ rr
        add!(βk,p,1,r) #p = beta*p + r

        Σ = leftvec ⋅ r

        for j=1:N
            ρkj = ρ0[j]
            if abs(ρkj) < eps
                continue
            end
            ρkmj =ρm[j]
            ρp[j] = ρkj*ρkmj*αm/(ρkmj*αm*(1.0+αk*vec_β[j])+αk*βm*(ρkmj-ρkj))
            αkj = (ρp[j]/ρkj)*αk
            θ[j] += αkj*Π[j]

            #add!(vec_x[j],αkj,vec_p[j])
            βkj = (ρp[j]/ρkj)^2*βk
            Π[j] = βkj*Π[j] + ρp[j]*Σ

            #mul!(Π,ρp[j],Σ,1,βkj)
            #Π[:] =
            #add!(βkj,vec_p[j],ρp[j],r)

        end

        ρm[:] = ρ0[:]
        ρ0[:] = ρp[:]
        αm = αk
        βm = βk

        ρMAX = maximum(abs.(ρp))^2
        residual = abs(rr*ρMAX)
        println_verbose_level3(verbose,"$i-th eps: $residual")

        if abs(residual) < eps
            println_verbose_level3(verbose,"Converged at $i-th step. eps: $residual")
            println_verbose_level3(verbose,"--------------------------------------")
            return θ
        end

    end

    error("""
    The shifted CG is not converged! with maxsteps = $(maxsteps)
    residual is $residual
    maxsteps should be larger.""")

end

function shiftedbicg_2003(σ,A,b;maxsteps=3000,eps = 1e-15,verboselevel=2)
    N = length(σ)
    btype = typeof(b)
    q = zero(b)
    x = zero(b)
    r = deepcopy(b)
    rt = deepcopy(b)
    u = zero(b)
    ut = zero(b)
    vec_u = Vector{btype}(undef,N)
    vec_x = Vector{btype}(undef,N)
    for j=1:N
        vec_u[j] = zero(b)
        vec_x[j] = zero(b)
    end
    ρkold = ones(ComplexF64,N)
    ρk = ones(ComplexF64,N)
    πkold = ones(ComplexF64,N)
    πk = ones(ComplexF64,N)
    πknew = ones(ComplexF64,N)
    ρold = 1
    αold = 1
    residual = 0

    for i=1:maxsteps
        ρ = dot(rt,r)
        β = -ρ/ρold
        axpby!(1,r,-β,u)
        axpby!(1,rt,-β',ut)
        mul!(q,A,u)
        qrt = dot(rt,q)
        α = ρ/qrt
        axpy!(α,u,x)

        for k=1:N
            πknew[k] = (1+α*σ[k])*πk[k]+ (α*β/αold)*(πkold[k]- πk[k])
            βk = (πkold[k]/πk[k])^2*β
            αk =  (πk[k]/πknew[k])*α
            axpby!(1/πk[k],r,-βk,vec_u[k])
            axpy!(αk,vec_u[k],vec_x[k])
        end

        axpy!(-α,q,r)
        mul!(q,A',ut)
        axpy!(-α',q,rt)
        αold = α
        ρold  = ρ

        residual = dot(r,r)
        if verboselevel == 3
            println("$i-th step: ",residual)
        end

        if abs(residual) < eps
            if verboselevel >= 2
                println("Converged at $i-th step. eps: $residual")
                println("--------------------------------------")
            end
            return vec_x
        end

        for k=1:N
            πkold[k] = πk[k]
            πk[k] = πknew[k]
        end

    end

    error("bicg is not converged. The residual is $residual")

end


text:
module SmileExtra

using Smile

export adjacency_matrix, adjacency_matrix_to_net

function adjacency_matrix(net::Network)
	# returns a BitMatrix in which index [parent,child] is true if parent->child is true

	allnodes = get_all_nodes(net)
	n = length(allnodes)
	adj = falses(n,n)
	for p in allnodes
		for c in get_children(net, p)
			adj[p+1,c+1] = true
		end
	end
	adj::BitMatrix
end
function adjacency_matrix_to_net(adj::BitMatrix, dset::Dataset)

	n = size(adj,1)
	@assert(n == size(adj,2))
	@assert(n == get_number_of_variables(dset))

	pat = Pattern()
	set_size(pat, n)

	for p = 1 : n
		for c = 1 : n
			if adj[p,c]
				set_edge(pat, p, c, DSL_EDGETYPE_DIRECTED)
			end
		end
	end

	@assert(is_DAG(pat))

	net = Network()
	to_network(pat, dset, net)

	net
end

include("statistics.jl")
include("laplace_smoothing.jl")

end # end module


text:
# t3r3: fixed
Fixed(body1::AbstractBody{T},body2,p1,p2;offset = Quaternion{T}()) where T = Translational3{T}(body1, body2, p1 = p1, p2 = p2), Rotational3{T}(body1, body2, offset = offset)

# t2r3: prismatic
Prismatic(body1::AbstractBody{T},body2,p1,p2,axis;offset = Quaternion{T}()) where T = Translational2{T}(body1, body2, p1 = p1, p2 = p2, axis = axis), Rotational3{T}(body1, body2, offset = offset)

# t1r3: ?
# ?

# t0r3: fixed orientation (chicken's head)
FixedOrientation(body1::AbstractBody{T},body2;offset = Quaternion{T}()) where T = Translational0{T}(body1, body2), Rotational3{T}(body1, body2, offset = offset)

# t3r2: Revolute
Revolute(body1::AbstractBody{T},body2,p1,p2,axis;offset = Quaternion{T}()) where T = Translational3{T}(body1, body2, p1 = p1, p2 = p2), Rotational2{T}(body1, body2, axis = axis, offset = offset)

# t2r2: Cylindrical
Cylindrical(body1::AbstractBody{T},body2,p1,p2,axis;offset = Quaternion{T}()) where T = Translational2{T}(body1, body2, p1 = p1, p2 = p2, axis = axis), Rotational2{T}(body1, body2, axis = axis, offset = offset)

# t1r2: ?
# ?

# t0r2: ?
# ?

# What is r1?
# t3r1: universal?
# ?
# t2r1: ?
# ?
# t1r1: ?
# ?
# t0r1: ?
# ?

# t3r0: Spherical
Spherical(body1::AbstractBody{T},body2,p1,p2) where T = Translational3{T}(body1, body2, p1 = p1, p2 = p2), Rotational0{T}(body1, body2)

# t2r0: Point-on-Line
CylindricalFree(body1::AbstractBody{T},body2,p1,p2,axis) where T = Translational2{T}(body1, body2, p1 = p1, p2 = p2, axis = axis), Rotational0{T}(body1, body2)

# t1r0: Planar
Planar(body1::AbstractBody{T},body2,p1,p2,axis) where T = Translational1{T}(body1, body2, p1 = p1, p2 = p2, axis = axis), Rotational0{T}(body1, body2)

# t0r0: Free
OriginConnection(body1::AbstractBody{T},body2) where T = Translational0{T}(body1, body2), Rotational0{T}(body1, body2)


text:
#=
This shows basic data analysis using the Card data
on returns to education

See the github pages of the packages for documentation

=#

using CSV, DataFrames, DataFramesMeta, DelimitedFiles, StatsPlots, Econometrics
pyplot(reuse=false, show=true)
card = CSV.read("../Data/card.csv") # how to read CSV data into DataFrame
size(card)
names(card)
# see the first 6 lines
show(first(card,6))
# see some statistics using DataFrames functionality
show(describe(card))
# see some statistics using Econometrics functionality
println()
dstats(convert(Matrix{Float64},card))

# a simple way of looking at geographical and racial discrimination
# select north and white
nw = @where(card, :south .== 0, :black .== 0)
@df nw scatter(:educ, :wage,label="North-White")
# select south and black
sb = @where(card, :south .== 1, :black .== 1)
@df sb scatter!(:educ, :wage, label="South-Black",title="Wage vs. Educ",xlabel="Educ", ylabel="Wage")

# the following is similar, but maybe better, as it plots all points
@df card scatter(:educ, :wage, zcolor=:south+:black)

##
# density plots by groups
@df card density(:wage, group=(:black), fill=true, fillalpha=0.3,legend=:topright)

##
# violin and box plots
@df card boxplot(:wage) # puts is beside
@df card violin!(:wage)
# violin plot of wage, by education
@df card violin(:educ, :wage)
# add a boxplot on top
@df card boxplot!(:educ, :wage)
##
# convert the dataframe to an ordinary array,
# as this is usually better for more advanced
# modeling
data = convert(Matrix{Float64}, card)
writedlm("card.ascii", data) # write to plain ASCII file

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
module FoodCachingFitting

using Distributed, CMAEvolutionStrategy, SpecialFunctions, Serialization, LibGit2, Unitful
@everywhere using FoodCachingExperiments, FoodCachingModels, Distances,
                  Random, DataFrames
import FoodCachingModels: Population, setparameters!, beta, truncnorm
import FoodCachingExperiments: EXPERIMENTS, CLAYTON0103_EXPERIMENTS
import CMAEvolutionStrategy: NoiseHandling, population_mean, Optimizer

export simulate, logp_hat, fit, run_experiments, resume_fit

function simulate(population, experiments; N = 500, rng = Random.GLOBAL_RNG)
    seeds = rand(rng, 0:typemax(UInt), N)
    setup = shuffle!(rng, collect(Iterators.product(seeds, experiments)))
    @distributed (vcat) for (seed, e) in setup
        Random.seed!(seed)
        ms = rand(population, nbirds(e))
        data = run!(e, ms)
        res = results(e, data)
        DataFrame(seed = [seed], experiment = [e], results = [res])
    end
end

function logp_hat(population, experiments;
                  N = 500, k = 5, metric = Distances.Euclidean(),
                  rng = Random.GLOBAL_RNG)
    samples = simulate(population, experiments; N, rng)
    sum(combine(groupby(samples, :experiment),
                df -> DataFrame(ll = _logp_hat(target(df.experiment[1]),
                                               df.results, k; metric))).ll)
end
logvol(::Distances.Euclidean, n, ::Any) = log(π)*(n/2) - lgamma(n/2 + 1)
function logvol(m::Distances.WeightedEuclidean, n, x)
    logvol(Distances.Euclidean(), n, x) - sum(log.(sqrt.(m.weights)))
end
function _logp_hat(x, samples, k; metric = Distances.Euclidean())
    d = metric.(Ref(x), samples)
    sort!(d)
    _logp_hat(d, x, k, metric)
end
function _logp_hat(d, x, k, metric)
    n = length(x)
    M = length(d)
    log(k/M) - logvol(metric, n, x) - log(d[k]) * n
end

function noisehandler(n, N; minN = 20, maxN = 500, alphaN = 1.05)
    cb = minN == maxN ?
         s -> s > 0 :
         s -> begin
            if s > 0
                if N[] < maxN
                    N[] = min(maxN, round(Int, N[] * alphaN))
                    @show N[]
                end
                return true
            else
                if N[] > minN
                    N[] = max(minN, round(Int, N[] / alphaN^.25))
                    @show N[]
                end
                return false
            end
        end
    NoiseHandling(n, callback = cb)
end

function checkpoint_saver(population, simname, f, flog;
                          saveevery = 3600, sigma_threshold = 10^3,
                          data_dir = DATADIR, t = time())
    (o, _...) -> begin
        sigma = CMAEvolutionStrategy.sigma(o.p)
        if time() - t > saveevery || sigma > sigma_threshold
            setparameters!(population, population_mean(o))
            dict = Dict{Symbol, Any}(:rev => __REV__)
            save(joinpath(DATADIR, simname), population, dict)
            serialize(joinpath(DATADIR, "o_$simname.dat"),
                      Dict(:optimizer => o, :func => f, :rev => __REV__))
            flush(flog[])
            if sigma > sigma_threshold
                error("Sigma is larger than $sigma_threshold.")
            end
            t = time()
        end
    end
end

function simname(model, experiments, id, rev = __REV__)
    e = experiments == ALLEXPERIMENTS ? "all" : join(string.(experiments), "_")
    "$(model)_$(e)_$(id)_$rev"
end

default_lower(m::Population{<:Any, typeof(beta)}) =
    [fill(-10^3, length(m.m)); fill(.5414, length(m.s))]
default_lower(m::Population{<:Any, typeof(truncnorm)}) =
    fill(-20, length(m.m) + length(m.s))
default_upper(m::Population{<:Any, typeof(beta)}) =
    fill(10^3, length(m.m) + length(m.s))
default_upper(m::Population{<:Any, typeof(truncnorm)}) =
    fill(20, length(m.m) + length(m.s))

function optimizer(; model,
                   N0 = 20,
                   id = "0",
                   experiments = ALLEXPERIMENTS,
                   saveevery = 3600,
                   sigma_threshold = 10^3,
                   sigma0 = .1,
                   x0 = nothing,
                   lower = :default,
                   upper = :default,
                   seed = time_ns(),
                   kwargs...)
    name = simname(model, experiments, id)
    population = model(; experiments, kwargs...)
    x0 = x0 === nothing ? [population.m; population.s === nothing ? Float64[] : population.s] : x0
    if lower == :default
        lower = default_lower(population)
    elseif lower == -Inf
        lower = nothing
    end
    if upper == :default
        upper = default_upper(population)
    elseif upper == Inf
        upper = nothing
    end
    flog = open(joinpath(DATADIR, "log", "$name.log"), "a+")
    redirect_stdout(flog)
    println(join(ARGS, " "))
    flush(flog)
    N = Ref(N0)
    rng = MersenneTwister(seed)
    f = x -> -logp_hat(setparameters!(population, x), experiments;
                       N = N[], rng = rng)
    callback = checkpoint_saver(population, name, f, Ref(flog);
                                saveevery, sigma_threshold)
    noise_handling = noisehandler(length(x0), N)
    (o = Optimizer(x0, sigma0; noise_handling, callback, seed, lower, upper), f = f)
end

function fit(; kwargs...)
    o, f = optimizer(; kwargs...)
    CMAEvolutionStrategy.run!(o, f)
end

function resume_fit(simname)
    flog = open(joinpath(DATADIR, "log", "$simname.log"), "a+")
    redirect_stdout(flog)
    sim = deserialize(joinpath(DATADIR, "o_$simname.dat"))
    o = sim[:optimizer]
    f = sim[:func]
    f.population.constructor = FoodCachingModels.init(f.population.p)
    o.logger.callback.flog[] = flog
    CMAEvolutionStrategy.run!(o, f)
end

function run_experiments(simname, experiments, N; seed = time_ns())
    population = load(joinpath(DATADIR, simname))
    simulate(population, experiments, N = N, rng = MersenneTwister(seed))
end
function run_experiments(; models, experiments,
                           id, N, seed = time_ns(), rev = __REV__)
    vcat([begin
              res = run_experiments(simname(model, exs, id, rev), exs, N; seed)
              res.model = fill(Symbol(model), length(res.results))
              res
          end
          for (model, exs) in Iterators.product(models, experiments)]...)
end

function parse_args()
    kwargs = Dict{Symbol, Any}()
    for arg in ARGS
        eq = findfirst(isequal('='), arg)
        kwargs[Symbol(arg[1:eq-1])] = eval(Meta.parse(arg[eq+1:end]))
    end
    kwargs
end

__REV__ = ""
const DATADIR = joinpath(@__DIR__, "..", "data")
const ALLEXPERIMENTS = setdiff(keys(EXPERIMENTS), CLAYTON0103_EXPERIMENTS)

function __init__()
    dir = joinpath(@__DIR__, "..")
    LibGit2.isdirty(GitRepo(dir)) && @warn "git repo is dirty."
    global __REV__ = LibGit2.head(dir)[1:7]
end
end


text:
    """"
    hold miscellaneous functions with no downstream dependencies.
    """

    """
        avoid_activations(tensor_shaped_like_weights)
    randn! every weight as long as it's not an activation (defined as a 1-tensor, since normal weights are arranged as matrices between layers)
    """
    function avoid_activations(noiseTemp)
        [randn!(t) for t in noiseTemp if length(size(t)) > 1]
    end


text:
using Documenter, TableFrames

makedocs(
    modules=[TableFrames],
    format=:html,
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/rofinn/TableFrames.jl/blob/{commit}{path}#L{line}",
    sitename="TableFrames.jl",
    authors="Rory Finnegan",
    assets=[],
)

deploydocs(
    repo = "github.com/rofinn/TableFrames.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
export Bnum, Bbound

immutable Bnum <: Real
    num::BigFloat
    open::Bool
end

immutable Bbound <: Real
    lo::Bnum
    hi::Bnum
end

Bnum(x::Real,open::Bool) = Bnum(BigFloat(x),open)
Bnum(x::Real) = Bnum(x,false)

convert(::Type{Bnum},x::Bnum) = x
convert(::Type{Bnum},x::Real) = Bnum(x)

import Base.MPFR: to_mpfr, ClongMax

-(x::Bnum) = Bnum(-x.num,x.open)

function +(x::Bnum, y::Bnum, r::RoundingMode)
    znum = BigFloat()
    inex = ccall((:mpfr_add,:libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &znum, &x.num, &y.num, to_mpfr(r)) != 0

    zopen = (x.open | y.open | inex) & (isfinite(x.num) | x.open) & (isfinite(y.num) | y.open)
    Bnum(znum, zopen)
end

function -(x::Bnum, y::Bnum, r::RoundingMode)
    znum = BigFloat()
    inex = ccall((:mpfr_sub,:libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &znum, &x.num, &y.num, to_mpfr(r)) != 0

    zopen = (x.open | y.open | inex) & (isfinite(x.num) | x.open) & (isfinite(y.num) | y.open)
    Bnum(znum, zopen)
end

function *(x::Bnum, y::Bnum, r::RoundingMode)
    znum = BigFloat()
    inex = ccall((:mpfr_mul,:libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &znum, &x.num, &y.num, to_mpfr(r)) != 0

    zopen = (x.open | y.open | inex) & (isfinite(x.num) & (x.num!=0) | x.open) & (isfinite(y.num) & (y.num!=0) | y.open)
    Bnum(znum, zopen)
end

function /(x::Bnum, y::Bnum, r::RoundingMode)
    znum = BigFloat()
    inex = ccall((:mpfr_div,:libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &znum, &x.num, &y.num, to_mpfr(r)) != 0

    zopen = (x.open | y.open | inex) & (isfinite(x.num) & (x.num!=0) | x.open) & (!isfinite(y.num) | y.open)
    Bnum(znum, zopen)
end

function sqrt(x::Bnum, r::RoundingMode)
    znum = BigFloat()
    inex = ccall((:mpfr_sqrt,:libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &znum, &x.num, to_mpfr(r)) != 0

    zopen = (x.open | inex)
    Bnum(znum, zopen)
end

function ^(x::Bnum, y::ClongMax, r::RoundingMode)
    znum = BigFloat()
    inex = ccall((:mpfr_pow_si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &znum, &x.num, y, to_mpfr(r)) != 0

    zopen = y != 0 && (x.open | inex)
    Bnum(znum, zopen)
end

import Base: min, max
# Technically not correct as we need to know the orientation, but assumes
# min is being used for lower bound, and max for upper bound
function min(x::Bnum,y::Bnum)
    if x.num == y.num
        Bnum(x.num,x.open & y.open)
    elseif x.num < y.num
        x
    else
        y
    end
end
function max(x::Bnum,y::Bnum)
    if x.num == y.num
        Bnum(x.num,x.open & y.open)
    elseif x.num > y.num
        x
    else
        y
    end
end

import Base: isnan
isnan(x::Bnum) = isnan(x.num)

function convert(::Type{Bbound},x::Real)
    c = convert(Bnum,x)
    Bbound(c,c)
end

function print(io::IO, b::Bbound)
    nlo = b.lo.num
    flo = Float64(nlo)
    if flo == 0
        # don't print signed zero
        flo = abs(flo)
    end

    nhi = b.hi.num
    fhi = Float64(nhi)
    if fhi == 0
        fhi = abs(fhi)
    end

    if nlo == nhi
        print(io,nlo==flo?flo:nlo)
    else
        print(io,b.lo.open?'(':'[',nlo==flo?flo:nlo,',',nhi==fhi?fhi:nhi,b.hi.open?')':']')
    end
end

show(io::IO, b::Bbound) = print(io,typeof(b),'\n',b)

-(x::Bbound) = Bbound(-x.hi,-x.lo)

function +(x::Bbound, y::Bbound)
    Bbound(+(x.lo,y.lo,RoundDown), +(x.hi,y.hi,RoundUp))
end
function -(x::Bbound, y::Bbound)
    Bbound(-(x.lo,y.hi,RoundDown), -(x.hi,y.lo,RoundUp))
end

function *(x::Bbound, y::Bbound)
    if isposz(x) # x >= 0
        if isposz(y)
            Bbound(*(x.lo,y.lo,RoundDown),*(x.hi,y.hi,RoundUp))
        elseif isnegz(y)
            Bbound(*(x.hi,y.lo,RoundDown),*(x.lo,y.hi,RoundUp))
        else
            Bbound(*(x.hi,y.lo,RoundDown),*(x.hi,y.hi,RoundUp))
        end
    elseif isnegz(x) # x <= 0
        if isposz(y)
            Bbound(*(x.lo,y.hi,RoundDown),*(x.hi,y.lo,RoundUp))
        elseif isnegz(y)
            Bbound(*(x.hi,y.hi,RoundDown),*(x.lo,y.lo,RoundUp))
        else
            Bbound(*(x.lo,y.hi,RoundDown),*(x.lo,y.lo,RoundUp))
        end
    else
        if isposz(y)
            Bbound(*(x.lo,y.hi,RoundDown),*(x.hi,y.hi,RoundUp))
        elseif isnegz(y)
            Bbound(*(x.hi,y.lo,RoundDown),*(x.lo,y.lo,RoundUp))
        else
            Bbound(min(*(x.lo,y.hi,RoundDown),*(x.hi,y.lo,RoundDown)),
                     max(*(x.lo,y.lo,RoundUp),*(x.hi,y.hi,RoundUp)))
        end
    end
end

function /(x::Bbound, y::Bbound)
    if ispos(y) # b strictly positive
        if isposz(x)
            Bbound(/(x.lo,y.hi,RoundDown), /(x.hi,y.lo,RoundUp))
        elseif isnegz(x)
            Bbound(/(x.lo,y.lo,RoundDown), /(x.hi,y.hi,RoundUp))
        else
            Bbound(/(x.lo,y.lo,RoundDown), /(x.hi,y.lo,RoundUp))
        end
    elseif isneg(y)
        if isposz(x)
            Bbound(/(x.hi,y.hi,RoundDown), /(x.lo,y.lo,RoundUp))
        elseif isnegz(x)
            Bbound(/(x.hi,y.lo,RoundDown), /(x.lo,y.hi,RoundUp))
        else
            Bbound(/(x.hi,y.hi,RoundDown), /(x.lo,y.hi,RoundUp))
        end
    else
        Bbound(NaN)
    end
end

function abs(x::Bbound)
    if isposz(x)
        x
    elseif isnegz(x)
        -x
    else
        Bbound(zero(BigFloat),max(-x.lo,x.hi))
    end
end

function sqrt(x::Bbound)
    if isposz(x)
        Bbound(sqrt(x.lo,RoundDown),sqrt(x.hi,RoundUp))
    else
        Bbound(NaN)
    end
end

function ^(x::Bbound,y::Integer)
    if y == 0
        Bbound(1)
    elseif y > 0
        if isodd(y) || isposz(x)
            Bbound(^(x.lo,y,RoundDown),^(x.hi,y,RoundUp))
        elseif isnegz(x)
            Bbound(^(x.hi,y,RoundDown),^(x.lo,y,RoundUp))
        else
            Bbound(zero(BigFloat),^(max(-x.lo,x.hi),y,RoundUp))
        end
    else
        if isodd(y) || isposz(x)
            Bbound(^(x.hi,y,RoundDown),^(x.lo,y,RoundUp))
        elseif isnegz(x)
            Bbound(^(x.lo,y,RoundDown),^(x.hi,y,RoundUp))
        else
            Bbound(NaN)
        end
    end
end

isposz(x::Bbound) = x.lo.num >= 0
isnegz(x::Bbound) = x.hi.num <= 0
ispos(x::Bbound) = x.lo.num > 0 || x.lo.num == 0 && x.lo.open
isneg(x::Bbound) = x.hi.num < 0 || x.hi.num == 0 && x.hi.open

==(x::Bnum,y::Bnum) = x.num == y.num && x.open == y.open
==(x::Bbound,y::Bbound) = x.lo == y.lo && x.hi == y.hi

<(x::Bbound,y::Bbound) = x.hi.num < y.lo.num || x.hi.num == y.lo.num && (x.hi.open | y.lo.open)
<=(x::Bbound,y::Bbound) = x.hi.num <= y.lo.num

function in(x::Bbound,y::Bbound)
    (y.lo.num < x.lo.num || (y.lo.num == x.lo.num && (!y.lo.open || x.lo.open))) &&
    (x.hi.num < y.hi.num || (x.hi.num == y.hi.num && (!y.hi.open || x.hi.open)))
end


text:
using ChainedArrays
using Test

@testset "ChainedArrays.jl" begin
    v = ChainedVector([1, 2, 3], [4], [5, 6, 7, 8, 9, 10])

    @test all(v .== 1:10)
    @test collect(x for x in v) == 1:10

    @test v[7] == 7

    @test v[1:5] == 1:5

    @test (v[3] = 10; v[3] == 10)

    @test (v[4:10] .= 1:7; v[4:10] == 1:7)
end


text:
using RadeonProRender
using Documenter

DocMeta.setdocmeta!(RadeonProRender, :DocTestSetup, :(using RadeonProRender); recursive=true)

makedocs(; modules=[RadeonProRender], authors="Simon Danisch",
         repo="https://github.com/SimonDanisch/RadeonProRender.jl/blob/{commit}{path}#{line}",
         sitename="RadeonProRender.jl",
         format=Documenter.HTML(; prettyurls=get(ENV, "CI", "false") == "true",
                                canonical="https://SimonDanisch.github.io/RadeonProRender.jl", assets=String[]),
         pages=["Home" => "index.md"])

deploydocs(; repo="github.com/SimonDanisch/RadeonProRender.jl")


text:
using GeometryTypes
using LinearAlgebra

"""
    subdivide(msh::HomogenousMesh,f::Function)

Returns a subdived triangular mesh from passed mesh `msh` and interpolator function `f`. The interpolator `f` is expected to accept two vertex indicies of the edge and to return a tuple of coordinates of the middle point (as one wishes to define it).
"""
function subdivide(msh::HomogenousMesh,f::Function)
    edges = filter(x->x[1]<x[2],decompose(Face{2,Int},msh))
    epoint(v1,v2) = findfirst(x->x==(v2>v1 ? Face(v1,v2) : Face(v2,v1)), edges) + length(msh.vertices)

    newfaces = Face{3,Int}[]

    newvertices = copy(msh.vertices)
    resize!(newvertices,length(msh.vertices) + length(edges))

    for (v1,v2,v3) in msh.faces

        ev3 = epoint(v1,v2)
        ev1 = epoint(v2,v3)
        ev2 = epoint(v3,v1)

        ### Usually, does assignment twice but is important if the surface is not tightly connected
        newvertices[ev3] = Point(f(v1,v2)...)
        newvertices[ev1] = Point(f(v2,v3)...)
        newvertices[ev2] = Point(f(v3,v1)...)

        push!(newfaces,Face(v1,ev3,ev2))
        push!(newfaces,Face(v2,ev1,ev3))
        push!(newfaces,Face(v3,ev2,ev1))
        push!(newfaces,Face(ev1,ev2,ev3))
    end

    return HomogenousMesh(newvertices,newfaces)
end

"""
    unitsphere(n)

Returns a sphere mesh made out of icosahedron subdivided `n` times.
"""
function unitsphere(subdivissions::Int64)
    t = ( 1 + sqrt( 5 ) ) / 2;

    vertices = Point{3,Float64}[
        [ -1,  t,  0 ], [  1, t, 0 ], [ -1, -t,  0 ], [  1, -t,  0 ],
        [  0, -1,  t ], [  0, 1, t ], [  0, -1, -t ], [  0,  1, -t ],
        [  t,  0, -1 ], [  t, 0, 1 ], [ -t,  0, -1 ], [ -t,  0,  1 ]
    ]

    faces = Face{3,Int64}[
        [1, 12, 6], [1, 6, 2], [1, 2, 8], [1, 8, 11], [1, 11, 12], [2, 6, 10], [6, 12, 5],
        [12, 11, 3], [11, 8, 7], [8, 2, 9], [4, 10, 5], [4, 5, 3], [4, 3, 7], [4, 7, 9],
        [4, 9, 10], [5, 10, 6], [3, 5, 12], [7, 3, 11], [9, 7, 8], [10, 9, 2]
    ]

    msh = HomogenousMesh(vertices ./ sqrt(1+t^2),faces)

    function sf(msh,v1,v2)
        p1 = msh.vertices[v1]
        p2 = msh.vertices[v2]
        return (p1 + p2)/norm(p1+p2)
    end

    for i in 1:subdivissions
        msh = subdivide(msh,(v1,v2)->sf(msh,v1,v2))
    end

    return msh
end


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

# using Markdown
# using InteractiveUtils

# ╔═╡ ea96e5d7-7bc0-45cf-8bf1-4acfaf5507c9
#=╠═╡ notebook_exclusive
include("basics.jl")
  ╠═╡ notebook_exclusive =#

# ╔═╡ 0ba0aae0-5c29-449a-81db-a389baddf0dc
#=╠═╡ notebook_exclusive
#using PlutoDevMacros
  ╠═╡ notebook_exclusive =#

# ╔═╡ 2307b8e6-6308-4902-9921-135b85273f65
#=╠═╡ notebook_exclusive
import WhereTraits
  ╠═╡ notebook_exclusive =#

# ╔═╡ 9b6c7d3b-b70a-4cd4-9a77-f9b8d8501b9e
import MacroTools

# ╔═╡ 4495c328-5bff-41ec-97f3-2422ee8c7339
function _plutotraits(ex::Expr)
	if ex.head === :block && length(ex.args) == 2
		return _plutotraits(ex.args[2])
	elseif ex.head === :macrocall && ex.args[1] isa GlobalRef && ex.args[1].name === Symbol("@doc")
		ex.args[end], fname = _plutotraits(ex.args[end])
		return ex, fname
	else
		# We try to see if this expression is a function definition
		defdict = MacroTools.splitdef(ex)
		# If we reach this point, the expression is a funcdef or splitdef would error
		# Extract the linenumbernode from the definition
		body = defdict[:body]
		ln = body.args[1]
		# If the second arg is also a linenumbernode, remove the first on
		if body.args[2] isa LineNumberNode
			popfirst!(body.args)
		end
		# Add the call to traits
		# ex = Expr(:macrocall, Expr(:(.),:WhereTraits,QuoteNode(Symbol("@traits"))), ln, ex)
		ex = :(WhereTraits.@traits $(ex))
		# Change the line
		ex.args[2] = ln
		# push!(ex.args, ln, ex)
		return ex, defdict[:name]
	end
	error("The provided expression was not of the valid type")
end

# ╔═╡ b810a742-dda2-4bd0-b9d3-8b0b5ee7356c
"""
	@plutotraits f(args...;kwargs...) where {wargs..} = body...
Macro that allows to use the `WhereTraits.@traits` macro inside Pluto notebooks.

When ran inside a notebook, this reconstructs the expression `expr` one would get by calling
`WhereTraits.@traits f(args...;kwargs...) where {wargs..} = body...` but evaluates it in the top-level of the current workspace/module (using `Core.eval(__module__,Expr(:toplevel, expr))`).\\
This is needed to avoid Pluto's internal expression explorer throwing an error if multiple functions with the same name but different conditions are defined with `@traits` in separate cells.

The usual Pluto workaround of wrapping multiple definitions in a begin ... end does not work with @traits as every call to the @traits macro has to be executed at the top level.\\
This has the side-effect that function calls within the notebook are not re-computed automatically if the definition with `@plutotraits` is changed within the notebook.

When called outside of Pluto (or included in Pluto from another notebook using `@plutoinclude`), the macro simply resorts to calling directly the `Where.@traits` expression without resorting to the `Core.eval` call within the macro body.

See [`WhereTraits.@traits`](@ref) or [`https://github.com/schlichtanders/WhereTraits.jl`](https://github.com/schlichtanders/WhereTraits.jl) for understanding how to use the underlying `@traits` macro.
"""
macro plutotraits(expr)
	# Meta.dump(expr)
	ex, fname = try
		_plutotraits(expr)
	catch e
		# error("The macro only supports a single call to the `@traits` macro with optional docstrings (wrapped in a begin...end block))")
		error(e)
	end
	if is_notebook_local(String(__source__.file))
		Core.eval(__module__, Expr(:toplevel, ex))
		esc(fname)
	else
		esc(ex)
	end
end

# ╔═╡ 8f354575-6ca2-4f0f-9880-f50d785de8f9
export @plutotraits

# ╔═╡ d6e472a7-bbd3-4f18-aa27-7289e6c23eab
#=╠═╡ notebook_exclusive
@plutotraits begin
	"Test Documentation"
	g(a) where {iseven(a)} = "EVEN"
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ e3c59dc8-6894-487a-ae98-62594fd149ef
#=╠═╡ notebook_exclusive
@plutotraits g(a) where {isodd(a)} = "ODD"
  ╠═╡ notebook_exclusive =#

# ╔═╡ 51bdae6a-7749-4bee-a6d0-b1055fab5ece
#=╠═╡ notebook_exclusive
g(2)
  ╠═╡ notebook_exclusive =#

# ╔═╡ c826aa1b-fd70-452c-9716-3bfc6066e590
#=╠═╡ notebook_exclusive
g(1)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
WhereTraits = "c9d4e05b-6318-49cb-9b56-e0e2b0ceadd8"

[compat]
MacroTools = "~0.5.9"
WhereTraits = "~1.0.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0-rc2"
manifest_format = "2.0"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dce3e3fea680869eaa0b774b2e8343e9ff442313"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.40.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.3.0"

[[deps.DataTypesBasic]]
deps = ["Compat"]
git-tree-sha1 = "4cff12742dcd7a8639da323abaf6bd4722abc312"
uuid = "83eed652-29e8-11e9-12da-a7c29d64ffc9"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.ExprParsers]]
deps = ["Compat", "ProxyInterfaces", "SimpleMatch", "StructEquality"]
git-tree-sha1 = "03d3f97dad4bd2b10ca0febca5db6bef5e0b320a"
uuid = "c5caad1f-83bd-4ce8-ac8e-4b29921e994e"
version = "1.1.0"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.ProxyInterfaces]]
deps = ["Compat"]
git-tree-sha1 = "d85fef4db37288bd756c71d223e34406055b7414"
uuid = "9b3bf0c4-f070-48bc-ae01-f2584e9c23bc"
version = "1.0.0"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Setfield]]
deps = ["ConstructionBase", "Future", "MacroTools", "Requires"]
git-tree-sha1 = "fca29e68c5062722b5b4435594c3d1ba557072a3"
uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
version = "0.7.1"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.SimpleMatch]]
deps = ["Compat"]
git-tree-sha1 = "c1cc22bbe259ea4a159e30ff2cf5f01d378262e8"
uuid = "a3ae8450-d22f-11e9-3fe0-77240e25996f"
version = "1.0.0"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StructEquality]]
deps = ["Compat"]
git-tree-sha1 = "6e951cd0585cbe8f4ceb1cb09d69332b9484fec9"
uuid = "6ec83bb0-ed9f-11e9-3b4c-2b04cb4e219c"
version = "1.1.0"

[[deps.Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.WhereTraits]]
deps = ["Compat", "DataTypesBasic", "ExprParsers", "Markdown", "ProxyInterfaces", "Setfield", "SimpleMatch", "StructEquality", "Suppressor"]
git-tree-sha1 = "cb13380b76dbe6e68c433dcd097aca3325a26470"
uuid = "c9d4e05b-6318-49cb-9b56-e0e2b0ceadd8"
version = "1.0.0"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═0ba0aae0-5c29-449a-81db-a389baddf0dc
# ╠═2307b8e6-6308-4902-9921-135b85273f65
# ╠═9b6c7d3b-b70a-4cd4-9a77-f9b8d8501b9e
# ╠═ea96e5d7-7bc0-45cf-8bf1-4acfaf5507c9
# ╠═8f354575-6ca2-4f0f-9880-f50d785de8f9
# ╠═b810a742-dda2-4bd0-b9d3-8b0b5ee7356c
# ╠═4495c328-5bff-41ec-97f3-2422ee8c7339
# ╠═d6e472a7-bbd3-4f18-aa27-7289e6c23eab
# ╠═e3c59dc8-6894-487a-ae98-62594fd149ef
# ╠═51bdae6a-7749-4bee-a6d0-b1055fab5ece
# ╠═c826aa1b-fd70-452c-9716-3bfc6066e590
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
Network = SimpleNet

netparams = SimpleNetHP(
  width=200,
  depth_common=8,
  use_batch_norm=true,
  batch_norm_momentum=1.)

self_play = SelfPlayParams(
  num_games=10_000,
  num_workers=100,
  reset_mcts_every=4,
  mcts = MctsParams(
    num_iters_per_turn=400,
    cpuct=1.0,
    temperature=ConstSchedule(1.0),
    dirichlet_noise_ϵ=0.2,
    dirichlet_noise_α=1.0))

arena = ArenaParams(
  num_games=100,
  num_workers=100,
  use_gpu=false,
  reset_mcts_every=1,
  update_threshold=0.00,
  flip_probability=0.5,
  mcts = MctsParams(
    self_play.mcts,
    temperature=ConstSchedule(0.3),
    dirichlet_noise_ϵ=0.1))

learning = LearningParams(
  use_gpu=false,
  samples_weighing_policy=LOG_WEIGHT,
  l2_regularization=1e-4,
  optimiser=CyclicNesterov(
    lr_base=1e-3,
    lr_high=1e-2,
    lr_low=1e-3,
    momentum_high=0.9,
    momentum_low=0.8),
  batch_size=32,
  loss_computation_batch_size=2048,
  nonvalidity_penalty=1.,
  min_checkpoints_per_epoch=0,
  max_batches_per_checkpoint=5_000,
  num_checkpoints=1)

params = Params(
  arena=arena,
  self_play=self_play,
  learning=learning,
  num_iters=2,
  memory_analysis=MemAnalysisParams(
    num_game_stages=4),
  ternary_rewards=true,
  use_symmetries=true,
  mem_buffer_size=PLSchedule(80_000))

benchmark = [
  Benchmark.Duel(
    Benchmark.Full(self_play.mcts),
    Benchmark.MctsRollouts(self_play.mcts),
    num_games=400,
    num_workers=100,
    use_gpu=true,
    flip_probability=0.5),
  Benchmark.Duel(
    Benchmark.NetworkOnly(),
    Benchmark.MinMaxTS(depth=6, amplify_rewards=true, τ=1.),
    num_games=400,
    num_workers=100,
    use_gpu=true,
    flip_probability=0.5)]


text:
@testset "exnode" begin

    @testset "basic ops" begin
        nd = ExNode{:call}(:z, :(x + y))

        @test getvar(nd) == :z
        @test varname(nd) == :z
        @test varidxs(nd) == []
        @test getexpr(nd) == :(x + y)
        @test getvalue(nd) == nothing

        setvar!(nd, :(z[i]))
        @test getvar(nd) == :(z[i])
        @test varname(nd) == :z
        @test varidxs(nd) == [:i]

        setexpr!(nd, :(x[i] + y[i]))
        @test getexpr(nd) == :(x[i] + y[i])

        setvalue!(nd, [1, 2, 3.])
        @test getvalue(nd) == [1, 2, 3.]
    end

    @testset "simple node" begin
        nd = ExNode{:call}(:z, :(x + y))
        @test to_expr(nd) == :(z = x + y)
    end

    @testset "indexed node" begin
        nd = ExNode{:call}(:(z[i]), :(x[i] + y[i]))
        @test to_expr(nd) == :(z[i] = x[i] + y[i])
    end

    @testset "dependencies" begin
        @test dependencies(ExNode{:constant}(:z, 42)) == []

        @test dependencies(ExNode{:input}(:z, 42)) == []

        @test dependencies(ExNode{:(=)}(:z, :x)) == [:x]
        @test dependencies(ExNode{:(=)}(:(z[i]), :(x[i]))) == [:x]

        @test dependencies(ExNode{:call}(:z, :(x + y))) == [:x, :y]
        @test dependencies(ExNode{:call}(:(z[i]), :(x[i] + y[i]))) == [:x, :y]

        @test dependencies(ExNode{:bcast}(:z, :(f.(x)))) == [:x]
    end

    @testset "broadcasting" begin
        nd = ExNode{:bcast}(:z, :(f.(x)))
        @test isindexed(nd) == false

        nd = ExNode{:bcast}(:z, :(f.(x[i])))
        @test isindexed(nd) == true
    end

end


text:
isdefined(Base, :__precompile__) && __precompile__()

"""
Module for constructing self-contained power system objects.
"""
module PowerSystems

#################################################################################
# Exports

export System
export ConcreteSystem
export Bus
export LoadZones

export PowerSystemType
export Component
export Device
export Branch
export Injection
export Line
export MonitoredLine
export DCLine
export HVDCLine
export VSCDCLine
export Transformer2W
export TapTransformer
export PhaseShiftingTransformer

export Forecast
export Deterministic
export Scenarios
export Probabilistic

export Generator
export GenClasses

export HydroGen
export HydroFix
export HydroCurtailment
export HydroStorage
export TechHydro
export EconHydro

export RenewableGen
export TechRenewable
export EconRenewable
export RenewableFix
export RenewableCurtailment
export RenewableFullDispatch

export ThermalGen
export TechThermal
export EconThermal
export ThermalDispatch
export ThermalGenSeason

export ElectricLoad
export StaticLoad
export PowerLoad
export PowerLoadPF
export FixedAdmittance
export ControllableLoad
export InterruptibleLoad

export Storage
export GenericBattery

export Service
export Reserve
export ProportionalReserve
export StaticReserve
export Transfer

export parsestandardfiles
export parse_file
export ps_dict2ps_struct
export assign_ts_data
export read_data_files
export validate
export get_components
export get_mixed_components
export get_component_counts
export show_component_counts

#################################################################################
# Imports

import SparseArrays
import AxisArrays
import LinearAlgebra: LAPACK.getri!
import LinearAlgebra: LAPACK.getrf!
import LinearAlgebra: BLAS.gemm
import LinearAlgebra
import Dates
import TimeSeries
import DataFrames
import JSON
import CSV
import YAML

#################################################################################
# Includes

# supertype for all PowerSystems types
abstract type PowerSystemType end
abstract type Component <: PowerSystemType end
# supertype for "devices" (bus, line, etc.)
abstract type Device <: Component end
abstract type Injection <: Device end
# supertype for generation technologies (thermal, renewable, etc.)
abstract type TechnicalParams <: Component end

include("common.jl")

# Include utilities
include("utils/utils.jl")
include("utils/logging.jl")
include("utils/IO/base_checks.jl")

# PowerSystems models
include("models/topological_elements.jl")
include("models/forecasts.jl")
include("models/branches.jl")
#include("models/network.jl")

# Static types
include("models/generation.jl")
include("models/storage.jl")
include("models/loads.jl")
include("models/services.jl")

# Include Parsing files
include("parsers/pm_io.jl")
include("parsers/im_io.jl")
include("parsers/dict_to_struct.jl")
include("parsers/standardfiles_parser.jl")
include("parsers/cdm_parser.jl")
include("parsers/forecast_parser.jl")
include("parsers/pm2ps_parser.jl")

#Data Checks
include("utils/IO/system_checks.jl")
include("utils/IO/branchdata_checks.jl")

# Definitions of System
include("base.jl")
include("validation/powersystem.jl")

# Better printing
include("utils/print.jl")
include("utils/lodf_calculations.jl")

# Download test data
include("utils/data.jl")
import .UtilsData: TestData

end # module


text:
# """
#     gen_possible_hallucination(params::Video_Params, cat::Int64)
#
# This function takes a category and params and it returns the possible
# objects (as 2D Detections) of that category that could be detected
# """
# #hallucinate objects in 2D image
# @gen (static) function gen_possible_hallucination(params::Video_Params, cat::Int64)
#     x = @trace(uniform(0,params.image_dim_x), :x)
#     y = @trace(uniform(0,params.image_dim_y), :y)
#     return (x, y, cat)
# end
# possible_hallucination_map = Gen.Map(gen_possible_hallucination)

function within_frame(p::Detection2D)
    p[1] >= 0 && p[1] <= 256 && p[2] >= 0 && p[2] <= 256 #hard-codded frame size because I can't figure out how to use two arguments and filter
end

#check if p1 and p2 are withing radius r of each other. Euclidean space
function within_radius(p1::Detection2D, p2::Detection2D, r::Float64)
    sqrt((p1[1]-p2[1])^2 + (p1[2]-p2[2])^2) <= r
end

#first approximation
function update_alpha_beta(lesioned::Bool, alphas_old::Matrix{Int64}, betas_old::Matrix{Int64}, observations_2D, real_detections::Array{Detection2D})
    alphas = deepcopy(alphas_old)
    betas = deepcopy(betas_old)

    if lesioned
        return alphas, betas
    end

    observations_2D = filter!(within_frame, observations_2D)
    real_detections = filter!(within_frame, real_detections)

    #println(real_detections)

    #lets only do this by category
    #real_detections_cats = last.(real_detections)
    #observations_2D_cats = last.(observations_2D)

    #see if actually detected. update miss rate
    observations_2D_edited = deepcopy(observations_2D)
    #observations_2D_cats_edited = copy(observations_2D_cats)
    for i = 1:length(real_detections)
        real_detection = real_detections[i]
        cat = real_detection[3] #category
        alphas[cat, 2] = alphas[cat, 2] + 1 #increase alpha for detection/miss rate
        j = 1
        while j <= length(observations_2D_edited) #basically a for loop over observations_2D_edited while it changes sizes
            obs = observations_2D_edited[j]
            if obs[3] == cat && within_radius(real_detection, obs, 40.)#if same category and within a distance of each other. 40 matches std on multinormal distribution for detection location
                observations_2D_edited = deleteat!(observations_2D_edited, j)
                betas[cat, 2] = betas[cat, 2] + 1 #increase beta for detection/miss rate
                #keep j the same because something was deleted at j
            else
                j = j+1
            end
        end
    end

    #everything left in observations_2D_edited must have been hallucinated
    for i in 1:length(observations_2D_edited)
        alphas[observations_2D_edited[i][3], 1] = alphas[observations_2D_edited[i][3], 1] + 1
    end

    #update beta for hallucinations
    betas[:, 1] = betas[:, 1] .+ 1 #add one for each frame
    return alphas, betas
end

"""given a 3D detection, return BernoulliElement over a 2D detection"""
function render(params::Video_Params, camera_params::Camera_Params, object_3D::Object3D)
    cat = object_3D[4]
    object = Coordinate(object_3D[1], object_3D[2], object_3D[3])
    x, y = get_image_xy(camera_params, params, object)
    return (x, y, cat)
end

"""
    gen_camera(params::Video_Params)

Independently samples a camera location and camera focus from a
uniform distribution
"""
@gen (static) function gen_camera(params::Video_Params)
    #camera location
    camera_location_x = @trace(uniform(params.x_min,params.x_max), :camera_location_x)
    camera_location_y = @trace(uniform(params.y_min,params.y_max), :camera_location_y)
    camera_location_z = @trace(uniform(params.z_min,params.z_max), :camera_location_z)

    #camera focus focus
    camera_focus_x = @trace(uniform(params.x_min,params.x_max), :camera_focus_x)
    camera_focus_y = @trace(uniform(params.y_min,params.y_max), :camera_focus_y)
    camera_focus_z = @trace(uniform(params.z_min,params.z_max), :camera_focus_z)

    camera_params = Camera_Params(Coordinate(camera_location_x,camera_location_y,camera_location_z), Coordinate(camera_focus_x,camera_focus_y,camera_focus_z))
    return camera_params
end

"""
Generates the next frame given the current frame.

state is Tuple{Array{Any,1}, Matrix{Int64}, Matrix{Int64}}
"""
@gen (static) function frame_kernel(current_frame::Int64, state::Any, lesioned::Bool, params::Video_Params, v::Matrix{Real}, receptive_fields::Vector{Receptive_Field})
    ####Update 2D real objects

    ####Update camera location and pointing
    camera_params = @trace(gen_camera(params), :camera)

    #get locations of the objects in the image. basically want to input the list
    #of observations_3D [(x,y,z,cat), (x,y,z,cat)] and get out the [(x_image,y_image,cat)]
    scene = state[1]
    n_real_objects = length(scene)
    paramses = fill(params, n_real_objects)
    #vs = fill(v, n_real_objects)
    camera_paramses = fill(camera_params, n_real_objects)
    real_detections = map(render, paramses, camera_paramses, scene)
    real_detections = Array{Detection2D}(real_detections)
    #observations_2D will be what we condition on

    rfs_vec = get_rfs_vec(real_detections, params, v)

    #for loop over receptive fields
    #@show maximum(map(length, rfs_vec))
    #could re-write with map
    #@trace(Gen.Map(rfs)(rfs_vec), :observations_2D) #gets no method matching error
    observations_2D = @trace(rfs(rfs_vec), :observations_2D) #dirty shortcut because we only have one receptive field atm
    alphas, betas = update_alpha_beta(lesioned, state[2], state[3], observations_2D, real_detections)
    state = (scene, alphas, betas) #just keep sending the scene in.
    return state
end

frame_chain = Gen.Unfold(frame_kernel)

"""
Print
"""
function print_helper_alphas_betas(alphas::Matrix{Int64}, betas::Matrix{Int64})
    @show alphas
    @show betas
end

"""
Samples new values for lambda_fa.
"""
@gen (static) function update_lambda_fa(alpha::Int64, beta::Int64)
    #println("in update_lambda_fa")
    #println("alpha ", alpha)
    #println("beta ", beta)
    fa = @trace(gamma(alpha, 1/beta), :fa)
    #println("fa ", fa)
    return fa
end

"""
Samples new values for miss_rate.
"""
@gen (static) function update_miss_rate(alpha::Int64, beta::Int64)
    #println("in update_miss_rate")
    #println("alpha ", alpha)
    #println("beta ", beta)
    miss = @trace(beta(alpha, beta), :miss)
    #println("miss ", miss)
    return miss
end

"""
Samples a new v based on the previous v.
"""
@gen (static) function update_v_matrix(alphas::Matrix{Int64}, betas::Matrix{Int64})
    #v = Matrix{Real}(undef, dim(previous_v_matrix))
    #a = print_helper_alphas_betas(alphas, betas)
    fa = @trace(Map(update_lambda_fa)(alphas[:,1], betas[:,1]), :lambda_fa)
    miss = @trace(Map(update_miss_rate)(alphas[:,2], betas[:,2]), :miss_rate)
    #v[:, 1] = fa
    #v[:, 2] = miss
    v = hcat(fa, miss)
    return convert(Matrix{Real}, v)
    #return v
end

"""
Samples a new scene and a new v_matrix.
"""
@gen (static) function video_kernel(current_video::Int64, v_matrix_state::Any, lesioned::Bool, num_frames::Int64, params::Video_Params, receptive_fields::Array{Receptive_Field, 1})
    #for the scene. scenes are completely independent of each other
    #println("current video ", current_video)

    #rfs_element = GeometricElement{Object3D}(params.p_objects, object_distribution, (params,))
    #rfs_element = RFSElements{Object3D}([rfs_element]) #need brackets because rfs has to take an array
    init_scene = @trace(object_distribution_gaussian(params), :init_scene)
    #make the observations
    previous_v_matrix = v_matrix_state[1]
    previous_alphas = v_matrix_state[2]
    previous_betas = v_matrix_state[3]
    init_state = (init_scene, previous_alphas, previous_betas)

    state = @trace(frame_chain(num_frames, init_state, lesioned, params, previous_v_matrix, receptive_fields), :frame_chain)
    alphas = state[end][2]#not sure num_frames or end is better for index
    betas = state[end][3]
    #for the metacognition.
    # println("alphas ", alphas)
    # println("betas ", betas)
    v_matrix = @trace(update_v_matrix(alphas, betas), :v_matrix)
    v_matrix_state = (v_matrix, alphas, betas)
    return v_matrix_state
end

#video_chain = Gen.Unfold(video_kernel)
"""Creates scene chain"""
video_chain = Gen.Unfold(video_kernel)

"""Creates frame chain"""
frame_chain = Gen.Unfold(frame_kernel)

export video_chain
export render
#export within_frame
export update_alpha_beta

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
using CmdStan, DiffEqBayes, OrdinaryDiffEq, ParameterizedFunctions,
      RecursiveArrayTools, Distributions, Test

println("One parameter case")
f1 = @ode_def begin
  dx = a*x - x*y
  dy = -3y + x*y
end a
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5]
prob1 = ODEProblem(f1,u0,tspan,p)
sol = solve(prob1,Tsit5())
t = collect(range(1,stop=10,length=10))
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
priors = [Truncated(Normal(1.5,0.1),0,2)]

bayesian_result = stan_inference(prob1,t,data,priors;num_samples=300,
                                 num_warmup=500,likelihood=Normal)

@test mean(get(bayesian_result.chains,:theta1).theta1') ≈ 1.5 atol=3e-1

println("Four parameter case")
f1 = @ode_def begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob1 = ODEProblem(f1,u0,tspan,p)
sol = solve(prob1,Tsit5())
t = collect(range(1,stop=10,length=10))
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
priors = [Truncated(Normal(1.5,0.01),0,2),Truncated(Normal(1.0,0.01),0,1.5),
          Truncated(Normal(3.0,0.01),0,4),Truncated(Normal(1.0,0.01),0,2)]

bayesian_result = stan_inference(prob1,t,data,priors;num_samples=100,num_warmup=500,vars =(DiffEqBayes.StanODEData(),InverseGamma(4,1)))
@test mean(get(bayesian_result.chains,:theta1).theta1') ≈ 1.5 atol=1e-1
@test mean(get(bayesian_result.chains,:theta2).theta2') ≈ 1.0 atol=1e-1
@test mean(get(bayesian_result.chains,:theta3).theta3') ≈ 3.0 atol=1e-1
@test mean(get(bayesian_result.chains,:theta4).theta4') ≈ 1.0 atol=1e-1


text:
include("setup.jl")

@testset "factorizing permutations" begin
    @testset "symmetric graphs" begin
        for _ = 1:100
            n = rand(3:10)
            G = rand(0:1, n, n)
            G .= G .⊻ G'
            @test G == G'
            p = symgraph_factorizing_permutation(G)
            @test is_modular_permutation(G, p)
            T = sort!(StrongModuleTree(G, p))
            @test sort!(strong_modules(T)) == sort!(all_strong_modules(G))
            for _ = 1:10
                p′ = symgraph_factorizing_permutation(G, shuffle(1:n))
                @test is_modular_permutation(G, p′)
                T′ = sort!(StrongModuleTree(G, p′))
                @test T == T′
            end
        end
    end

    @testset "tournament graphs" begin
        for _ = 1:1000
            n = rand(3:10)
            T = Int[i != j && rand(Bool) for i=1:n, j=1:n]
            T .= T .⊻ T' .⊻ tril(ones(Int,n,n),-1)
            @test is_tournament(T)
            p = tournament_factorizing_permutation(T)
            modules = all_modules(T)
            @test is_modular_permutation(T, p, modules=modules)
        end
    end

    @testset "directed graphs" begin
        for _ = 1:1000
            n = rand(3:10)
            G = Int[i != j && rand(Bool) for i=1:n, j=1:n]
            p = digraph_factorizing_permutation(G)
            @test is_modular_permutation(G, p)
            p = graph_factorizing_permutation(G)
            @test is_modular_permutation(G, p)
        end
    end
end

function test_permutations(G::AbstractMatrix, T::StrongModuleTree, N::Integer=100)
    n = LinearAlgebra.checksquare(G)
    for _ = 1:N
        p = shuffle(1:n)
        G′ = G[p,p]
        @test sort!(p[StrongModuleTree(G′)]) == T
    end
end

@testset "symmetric graph example" begin
    # from Habib & Paul 2009 (and Wikipedia)
    G = sparse(
        [1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5,
         5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7,
         8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 11],
        [2, 3, 4, 1, 4, 5, 6, 7, 1, 4, 5, 6, 7, 1, 2, 3, 5, 6, 7, 2,
         3, 4, 6, 7, 2, 3, 4, 5, 8, 9, 10, 11, 2, 3, 4, 5, 8, 9, 10,
         11, 6, 7, 9, 10, 11, 6, 7, 8, 10, 11, 6, 7, 8, 9, 6, 7, 8, 9],
        1
    )
    T = sort!(StrongModuleTree(G))
    @test repr(T) == "{1 ((2 3) 4) 5 (6 7) (8 9 (10 11))}"
    test_permutations(G, T)
end

@testset "directed graph example" begin
    # from Capelle, Habib & Montgolfier 2002
    G = sparse(
        [1, 1, 1, 2, 3, 3, 3, 4, 5, 5, 6, 6, 6, 6, 7, 7, 7, 8, 8, 10, 10,
         10, 10, 11, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14, 14, 14],
        [3, 4, 5, 1, 2, 4, 5, 2, 4, 2, 7, 8, 9, 10, 8, 9, 10, 9, 10, 11,
         12, 13, 14, 9, 10, 14, 9, 10, 11, 13, 9, 10, 11, 12, 9, 10, 11],
        1
    )
    T = sort!(StrongModuleTree(G))
    @test repr(T) == "({1 2 [3 4 5]} {[6 7 8] 9 10 {11 (12 13) 14}})"
    test_permutations(G, T)
end