text:
import Base.showerror
struct CycleError <: Exception end
Base.showerror(io::IO, e::CycleError) = print(io, "Cycles are not allowed in this function.")

struct DFS <: TraversalAlgorithm end

function traverse_graph!(
    g::AbstractGraph{U},
    ss::AbstractVector,
    alg::DFS,
    state::AbstractTraversalState,
    neighborfn::Function=outneighbors
    ) where U<:Integer

    n = nv(g)
    visited = falses(n)
    S = Vector{U}()
    sizehint!(S, length(ss))
    @inbounds for s in ss
        us = U(s)
        visited[us] = true
        push!(S, us)
        initfn!(state, us) || return false
    end
    while !isempty(S)
        v = S[end]
        previsitfn!(state, v) || return false
        new_neighbor_found = false
        @inbounds for i in neighborfn(g, v)
            visitfn!(state, v, i) || return false
            if !visited[i]  # find the first unvisited neighbor
                newvisitfn!(state, v, i) || return false
                new_neighbor_found = true
                visited[i] = true
                push!(S, i)
                break
            end
        end
        postvisitfn!(state, v) || return false
        if !new_neighbor_found  # no more new neighbors. Let's go to the next source vertex.
            postlevelfn!(state) || return false
            pop!(S)
        end
    end
    return true
end

mutable struct TopoSortState{T<:Integer} <: AbstractTraversalState
    vcolor::Vector{UInt8}
    verts :: Vector{T}
    w::T
end

# 1 = visited
# 2 = vcolor 2
# @inline initfn!(s::TopoSortState{T}, u) where T = s.vcolor[u] = one(T)
@inline function previsitfn!(s::TopoSortState{T}, u) where T
    s.w = 0
    return true
end
@inline function visitfn!(s::TopoSortState{T}, u, v) where T
    return s.vcolor[v] != one(T)
end
@inline function newvisitfn!(s::TopoSortState{T}, u, v) where T
    s.w = v
    return true
end
@inline function postvisitfn!(s::TopoSortState{T}, u) where T
    if s.w != 0
        s.vcolor[s.w] = one(T)
    else
        s.vcolor[u] = T(2)
        push!(s.verts, u)
    end
    return true
end

@traitfn function topological_sort(g::AG::IsDirected) where {T, AG<:AbstractGraph{T}}
    vcolor = zeros(UInt8, nv(g))
    verts = Vector{T}()
    state = TopoSortState(vcolor, verts, zero(T))
    for v in vertices(g)
        state.vcolor[v] != 0 && continue
        state.vcolor[v] = 1
        if !traverse_graph!(g, v, DFS(), state)
            throw(CycleError())
        end
    end
    return reverse(state.verts)
end

mutable struct CycleState{T<:Integer} <: AbstractTraversalState
    vcolor::Vector{UInt8}
    w::T
end

@inline function previsitfn!(s::CycleState{T}, u) where T
    s.w = 0
    return true
end
@inline function visitfn!(s::CycleState{T}, u, v) where T
    return s.vcolor[v] != one(T)
end
@inline function newvisitfn!(s::CycleState{T}, u, v) where T
    s.w = v
    return true
end
@inline function postvisitfn!(s::CycleState{T}, u) where T
    if s.w != 0
        s.vcolor[s.w] = one(T)
    else
        s.vcolor[u] = T(2)
    end
    return true
end

@traitfn function is_cyclic(g::AG::IsDirected) where {T, AG<:AbstractGraph{T}}
    vcolor = zeros(UInt8, nv(g))
    state = CycleState(vcolor, zero(T))
    @inbounds for v in vertices(g)
        state.vcolor[v] != 0 && continue
        state.vcolor[v] = 1
        !traverse_graph!(g, v, DFS(), state) && return true
    end
    return false
end


text:
"""
    Classic(;method=median)
    Classic(method)

Classic PSF subtraction using the median of entire data cube. If another statistic is desired, like the `mean`, it can be passed as an argument as long as it supports slicing multi-dimensional arrays.

# References
1. [Marois et al. 2006](http://adsabs.harvard.edu/abs/2006ApJ...641..556M) Angular Differential Imaging: A Powerful High-Contrast Imaging Technique
"""
struct Classic <: ADIAlgorithm
    method
end
Classic(;method=median) = Classic(method)

"""
    ADI.ClassicDesign(n, frame)

Output for the [`Classic`](@ref) algorithm which contains the static frame unrolled into a vector (with size `(1, Npx)`). [`reconstruct`](@ref) will tile this vector in a non-allocating way `n` times to appear like a flattened cube. [`ADI.design`](@ref) will return the static frame .
"""
struct ClassicDesign{FT} <: ADIDesign
    n::Int # number of frames in original data
    frame::FT
end

design(des::ClassicDesign) = des.frame

reconstruct(des::ClassicDesign) = TiledArray(des.frame, des.n)

function fit(alg::Classic, data::AbstractMatrix; ref=data, kwargs...)
    n = size(data, 1)
    return ClassicDesign(n, alg.method(ref, dims=1))
end

#######################################

struct TiledArray{T,N,AT<:AbstractArray{T,N},RT} <: AbstractArray{T,N}
    parent::AT
    tilerange::RT
end

Base.size(t::TiledArray) = map(length, axes(t))
Base.axes(t::TiledArray) = (t.tilerange, Base.tail(axes(t.parent))...)
Base.parent(t::TiledArray) = t.parent

TiledArray(parent::AbstractArray, n::Int) = TiledArray(parent, Base.OneTo(n))

Base.@propagate_inbounds function Base.getindex(t::TiledArray{T,N}, idxs::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(t, idxs...)
    return t.parent[firstindex(t.parent, 1), Base.tail(idxs)...]
end

ext: jl


text:
module Utils
import Dates
import FFTW
import Logging
import LibGit2
import Pidfile: mkpidlock
import HDF5
import Luna: @hlock, settings

subzero = '\u2080'
subscript(digit::Char) = string(Char(codepoint(subzero)+parse(Int, digit)))
subscript(num::AbstractString) = prod([subscript(chi) for chi in num])
subscript(num::Int) = subscript(string(num))

function git_commit()
    try
        repo = LibGit2.GitRepo(lunadir())
        commit = string(LibGit2.GitHash(LibGit2.head(repo)))
        LibGit2.isdirty(repo) && (commit *= " (dirty)")
        return commit
    catch
        "unavailable (Luna is not checkout out for development)"
    end
end

function git_branch()
    try
        repo = LibGit2.GitRepo(lunadir())
        n = string(LibGit2.name(LibGit2.head(repo)))
        branch = split(n, "/")[end]
        return branch
    catch
        "unavailable (Luna is not checkout out for development)"
    end
end

srcdir() = dirname(@__FILE__)

lunadir() = dirname(srcdir())

datadir() = joinpath(srcdir(), "data")

cachedir() = joinpath(homedir(), ".luna")

function sourcecode()
    src = dirname(@__FILE__)
    luna = dirname(src)
    out = "#= Date: $(Dates.now())\n"
    out *= "git branch: $(git_branch())\n"
    out *= "git commit: $(git_commit())\n"
    out *= "hostname: $(gethostname())\n"
    out *= "=#"
    for folder in (src, luna)
        for obj in readdir(folder)
            if isfile(joinpath(folder, obj))
                if split(obj, ".")[end] in ("md", "jl", "txt", "toml") # avoid binary files
                    out *= "\n" * "#" * "="^8 * obj * "="^8 * "#" * "\n"^2
                    open(joinpath(folder, obj), "r") do file
                        out *= read(file, String)
                    end
                end
            end
        end
    end
    return out
end

function FFTWthreads()
    if Threads.nthreads() == 1
        1
    else
        settings["fftw_threads"] == 0 ? 4*Threads.nthreads() : settings["fftw_threads"]
    end
end

function loadFFTwisdom()
    FFTW.set_num_threads(FFTWthreads())
    fpath = joinpath(cachedir(), "FFTWcache_$(FFTWthreads())threads")
    lockpath = joinpath(cachedir(), "FFTWlock")
    isdir(cachedir()) || mkpath(cachedir())
    if isfile(fpath)
        Logging.@info("Found FFTW wisdom at $fpath")
        pidlock = mkpidlock(lockpath)
        ret = FFTW.import_wisdom(fpath)
        close(pidlock)
    else
        Logging.@info("No FFTW wisdom found")
    end
end

function saveFFTwisdom()
    fpath = joinpath(cachedir(), "FFTWcache_$(FFTWthreads())threads")
    lockpath = joinpath(cachedir(), "FFTWlock")
    pidlock = mkpidlock(lockpath)
    isfile(fpath) && rm(fpath)
    isdir(cachedir()) || mkpath(cachedir())
    FFTW.export_wisdom(fpath)
    close(pidlock)
    Logging.@info("FFTW wisdom saved to $fpath")
end

function save_dict_h5(fpath, d; force=false, rmold=false)
    if isfile(fpath) && rmold
        rm(fpath)
    end

    function dict2h5(k::AbstractString, v, parent)
        if HDF5.haskey(parent, k) && !force
            error("Dataset $k exists in $fpath. Set force=true to overwrite.")
        end
        parent[k] = v
    end

    function dict2h5(k::AbstractString, v::BitArray, parent)
        if HDF5.haskey(parent, k) && !force
            error("Dataset $k exists in $fpath. Set force=true to overwrite.")
        end
        parent[k] = Array{Bool, 1}(v)
    end

    function dict2h5(k::AbstractString, v::Nothing, parent)
        if HDF5.haskey(parent, k) && !force
            error("Dataset $k exists in $fpath. Set force=true to overwrite.")
        end
        parent[k] = Float64[]
    end

    function dict2h5(k::AbstractString, v::AbstractDict, parent)
        if !HDF5.haskey(parent, k)
            subparent = HDF5.create_group(parent, k)
        else
            subparent = parent[k]
        end
        for (kk, vv) in pairs(v)
            dict2h5(kk, vv, subparent)
        end
    end

    @hlock HDF5.h5open(fpath, "cw") do file
        for (k, v) in pairs(d)
            dict2h5(k, v, file)
        end
    end
end

function save_dict_h5(fpath, t::NamedTuple; kwargs...)
    d = Dict{String, Any}()
    for (k, v) in pairs(t)
        d[string(k)] = v
    end
    save_dict_h5(fpath, d; kwargs...)
end

function load_dict_h5(fpath)
    isfile(fpath) || error("Error loading file $fpath: file does not exist")

    function h52dict(x::HDF5.Dataset)
        return read(x)
    end

    function h52dict(x::Union{HDF5.Group, HDF5.File})
        dd = Dict{String, Any}()
        for n in keys(x)
            dd[n] = h52dict(x[n])
        end
        return dd
    end

    d = @hlock HDF5.h5open(fpath) do file
        h52dict(file)
    end
end

end


text:
#Adding Julia Packages to Jupyter

Pkg.add("DataFrames")
Pkg.add("RDatasets")
quit();


text:
# ---
# title: 101. Symmetric Tree
# id: problem101
# author: zhwang
# date: 2022-02-17
# difficulty: Easy
# categories: Tree, Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/symmetric-tree/description/>
# hidden: true
# ---
#
# Given a binary tree, check whether it is a mirror of itself (ie, symmetric
# around its center).
#
# For example, this binary tree `[1,2,2,3,4,4,3]` is symmetric:
#
#
#
#         1
#        / \
#       2   2
#      / \ / \
#     3  4 4  3
#
#
#
#
# But the following `[1,2,2,null,3,null,3]` is not:
#
#
#
#         1
#        / \
#       2   2
#        \   \
#        3    3
#
#
#
#
# **Follow up:** Solve it both recursively and iteratively.
#
#
## @lc code=start
using LeetCode

is_symmetric_tree(root::TreeNode{Int})::Bool = is_mirror_tree(root.left, root.right)
function is_mirror_tree(
    t1::Union{TreeNode{Int},Nothing}, t2::Union{TreeNode{Int},Nothing}
)::Bool
    isnothing(t1) && isnothing(t2) && return true
    (isnothing(t1) || isnothing(t2)) && return false
    t1.val != t2.val && return false
    return is_mirror_tree(t1.left, t2.right) && is_mirror_tree(t1.right, t2.left)
end
## @lc code=end


text:
#==============================================================================#
# KinesisVideo.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/kinesisvideo-2017-09-30.normal.json
#==============================================================================#

__precompile__()

module KinesisVideo

using AWSCore

"""
    using AWSSDK.KinesisVideo.create_stream
    create_stream([::AWSConfig], arguments::Dict)
    create_stream([::AWSConfig]; StreamName=, <keyword arguments>)

    using AWSCore.Services.kinesisvideo
    kinesisvideo([::AWSConfig], "POST", "/createStream", arguments::Dict)
    kinesisvideo([::AWSConfig], "POST", "/createStream", StreamName=, <keyword arguments>)

# CreateStream Operation

Creates a new Kinesis video stream.

When you create a new stream, Kinesis Video Streams assigns it a version number. When you change the stream's metadata, Kinesis Video Streams updates the version.

`CreateStream` is an asynchronous operation.

For information about how the service works, see [How it Works](http://docs.aws.amazon.com/kinesisvideostreams/latest/dg/how-it-works.html).

You must have permissions for the `KinesisVideo:CreateStream` action.

# Arguments

## `DeviceName = ::String`
The name of the device that is writing to the stream.

**Note**
> In the current implementation, Kinesis Video Streams does not use this name.

## `StreamName = ::String` -- *Required*
A name for the stream that you are creating.

The stream name is an identifier for the stream, and must be unique for each account and region.

## `MediaType = ::String`
The media type of the stream. Consumers of the stream can use this information when processing the stream. For more information about media types, see [Media Types](http://www.iana.org/assignments/media-types/media-types.xhtml). If you choose to specify the `MediaType`, see [Naming Requirements](https://tools.ietf.org/html/rfc6838#section-4.2) for guidelines.

To play video on the console, the media must be H.264 encoded, and you need to specify this video type in this parameter as `video/h264`.

This parameter is optional; the default value is `null` (or empty in JSON).

## `KmsKeyId = ::String`
The ID of the AWS Key Management Service (AWS KMS) key that you want Kinesis Video Streams to use to encrypt stream data.

If no key ID is specified, the default, Kinesis Video-managed key (`aws/kinesisvideo`) is used.

For more information, see [DescribeKey](http://docs.aws.amazon.com/kms/latest/APIReference/API_DescribeKey.html#API_DescribeKey_RequestParameters).

## `DataRetentionInHours = ::Int`
The number of hours that you want to retain the data in the stream. Kinesis Video Streams retains the data in a data store that is associated with the stream.

The default value is 0, indicating that the stream does not persist data.

When the `DataRetentionInHours` value is 0, consumers can still consume the fragments that remain in the service host buffer, which has a retention time limit of 5 minutes and a retention memory limit of 200 MB. Fragments are removed from the buffer when either limit is reached.

# Returns

`CreateStreamOutput`

# Exceptions

`AccountStreamLimitExceededException`, `DeviceStreamLimitExceededException`, `ResourceInUseException`, `InvalidDeviceException`, `InvalidArgumentException` or `ClientLimitExceededException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/kinesisvideo-2017-09-30/CreateStream)
"""
@inline create_stream(aws::AWSConfig=default_aws_config(); args...) = create_stream(aws, args)

@inline create_stream(aws::AWSConfig, args) = AWSCore.Services.kinesisvideo(aws, "POST", "/createStream", args)

@inline create_stream(args) = create_stream(default_aws_config(), args)

"""
    using AWSSDK.KinesisVideo.delete_stream
    delete_stream([::AWSConfig], arguments::Dict)
    delete_stream([::AWSConfig]; StreamARN=, <keyword arguments>)

    using AWSCore.Services.kinesisvideo
    kinesisvideo([::AWSConfig], "POST", "/deleteStream", arguments::Dict)
    kinesisvideo([::AWSConfig], "POST", "/deleteStream", StreamARN=, <keyword arguments>)

# DeleteStream Operation

Deletes a Kinesis video stream and the data contained in the stream.

This method marks the stream for deletion, and makes the data in the stream inaccessible immediately.

To ensure that you have the latest version of the stream before deleting it, you can specify the stream version. Kinesis Video Streams assigns a version to each stream. When you update a stream, Kinesis Video Streams assigns a new version number. To get the latest stream version, use the `DescribeStream` API.

This operation requires permission for the `KinesisVideo:DeleteStream` action.

# Arguments

## `StreamARN = ::String` -- *Required*
The Amazon Resource Name (ARN) of the stream that you want to delete.

## `CurrentVersion = ::String`
Optional: The version of the stream that you want to delete.

Specify the version as a safeguard to ensure that your are deleting the correct stream. To get the stream version, use the `DescribeStream` API.

If not specified, only the `CreationTime` is checked before deleting the stream.

# Returns

`DeleteStreamOutput`

# Exceptions

`ClientLimitExceededException`, `InvalidArgumentException`, `ResourceNotFoundException` or `NotAuthorizedException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/kinesisvideo-2017-09-30/DeleteStream)
"""
@inline delete_stream(aws::AWSConfig=default_aws_config(); args...) = delete_stream(aws, args)

@inline delete_stream(aws::AWSConfig, args) = AWSCore.Services.kinesisvideo(aws, "POST", "/deleteStream", args)

@inline delete_stream(args) = delete_stream(default_aws_config(), args)

"""
    using AWSSDK.KinesisVideo.describe_stream
    describe_stream([::AWSConfig], arguments::Dict)
    describe_stream([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.kinesisvideo
    kinesisvideo([::AWSConfig], "POST", "/describeStream", arguments::Dict)
    kinesisvideo([::AWSConfig], "POST", "/describeStream", <keyword arguments>)

# DescribeStream Operation

Returns the most current information about the specified stream. You must specify either the `StreamName` or the `StreamARN`.

# Arguments

## `StreamName = ::String`
The name of the stream.

## `StreamARN = ::String`
The Amazon Resource Name (ARN) of the stream.

# Returns

`DescribeStreamOutput`

# Exceptions

`InvalidArgumentException`, `ResourceNotFoundException`, `ClientLimitExceededException` or `NotAuthorizedException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/kinesisvideo-2017-09-30/DescribeStream)
"""
@inline describe_stream(aws::AWSConfig=default_aws_config(); args...) = describe_stream(aws, args)

@inline describe_stream(aws::AWSConfig, args) = AWSCore.Services.kinesisvideo(aws, "POST", "/describeStream", args)

@inline describe_stream(args) = describe_stream(default_aws_config(), args)

"""
    using AWSSDK.KinesisVideo.get_data_endpoint
    get_data_endpoint([::AWSConfig], arguments::Dict)
    get_data_endpoint([::AWSConfig]; APIName=, <keyword arguments>)

    using AWSCore.Services.kinesisvideo
    kinesisvideo([::AWSConfig], "POST", "/getDataEndpoint", arguments::Dict)
    kinesisvideo([::AWSConfig], "POST", "/getDataEndpoint", APIName=, <keyword arguments>)

# GetDataEndpoint Operation

Gets an endpoint for a specified stream for either reading or writing. Use this endpoint in your application to read from the specified stream (using the `GetMedia` or `GetMediaForFragmentList` operations) or write to it (using the `PutMedia` operation).

**Note**
> The returned endpoint does not have the API name appended. The client needs to add the API name to the returned endpoint.

In the request, specify the stream either by `StreamName` or `StreamARN`.

# Arguments

## `StreamName = ::String`
The name of the stream that you want to get the endpoint for. You must specify either this parameter or a `StreamARN` in the request.

## `StreamARN = ::String`
The Amazon Resource Name (ARN) of the stream that you want to get the endpoint for. You must specify either this parameter or a `StreamName` in the request.

## `APIName = "PUT_MEDIA", "GET_MEDIA", "LIST_FRAGMENTS", "GET_MEDIA_FOR_FRAGMENT_LIST" or "GET_HLS_STREAMING_SESSION_URL"` -- *Required*
The name of the API action for which to get an endpoint.

# Returns

`GetDataEndpointOutput`

# Exceptions

`InvalidArgumentException`, `ResourceNotFoundException`, `ClientLimitExceededException` or `NotAuthorizedException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/kinesisvideo-2017-09-30/GetDataEndpoint)
"""
@inline get_data_endpoint(aws::AWSConfig=default_aws_config(); args...) = get_data_endpoint(aws, args)

@inline get_data_endpoint(aws::AWSConfig, args) = AWSCore.Services.kinesisvideo(aws, "POST", "/getDataEndpoint", args)

@inline get_data_endpoint(args) = get_data_endpoint(default_aws_config(), args)

"""
    using AWSSDK.KinesisVideo.list_streams
    list_streams([::AWSConfig], arguments::Dict)
    list_streams([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.kinesisvideo
    kinesisvideo([::AWSConfig], "POST", "/listStreams", arguments::Dict)
    kinesisvideo([::AWSConfig], "POST", "/listStreams", <keyword arguments>)

# ListStreams Operation

Returns an array of `StreamInfo` objects. Each object describes a stream. To retrieve only streams that satisfy a specific condition, you can specify a `StreamNameCondition`.

# Arguments

## `MaxResults = ::Int`
The maximum number of streams to return in the response. The default is 10,000.

## `NextToken = ::String`
If you specify this parameter, when the result of a `ListStreams` operation is truncated, the call returns the `NextToken` in the response. To get another batch of streams, provide this token in your next request.

## `StreamNameCondition = [ ... ]`
Optional: Returns only streams that satisfy a specific condition. Currently, you can specify only the prefix of a stream name as a condition.
```
 StreamNameCondition = [
        "ComparisonOperator" =>  "BEGINS_WITH",
        "ComparisonValue" =>  ::String
    ]
```

# Returns

`ListStreamsOutput`

# Exceptions

`ClientLimitExceededException` or `InvalidArgumentException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/kinesisvideo-2017-09-30/ListStreams)
"""
@inline list_streams(aws::AWSConfig=default_aws_config(); args...) = list_streams(aws, args)

@inline list_streams(aws::AWSConfig, args) = AWSCore.Services.kinesisvideo(aws, "POST", "/listStreams", args)

@inline list_streams(args) = list_streams(default_aws_config(), args)

"""
    using AWSSDK.KinesisVideo.list_tags_for_stream
    list_tags_for_stream([::AWSConfig], arguments::Dict)
    list_tags_for_stream([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.kinesisvideo
    kinesisvideo([::AWSConfig], "POST", "/listTagsForStream", arguments::Dict)
    kinesisvideo([::AWSConfig], "POST", "/listTagsForStream", <keyword arguments>)

# ListTagsForStream Operation

Returns a list of tags associated with the specified stream.

In the request, you must specify either the `StreamName` or the `StreamARN`.

# Arguments

## `NextToken = ::String`
If you specify this parameter and the result of a `ListTagsForStream` call is truncated, the response includes a token that you can use in the next request to fetch the next batch of tags.

## `StreamARN = ::String`
The Amazon Resource Name (ARN) of the stream that you want to list tags for.

## `StreamName = ::String`
The name of the stream that you want to list tags for.

# Returns

`ListTagsForStreamOutput`

# Exceptions

`ClientLimitExceededException`, `InvalidArgumentException`, `ResourceNotFoundException`, `NotAuthorizedException` or `InvalidResourceFormatException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/kinesisvideo-2017-09-30/ListTagsForStream)
"""
@inline list_tags_for_stream(aws::AWSConfig=default_aws_config(); args...) = list_tags_for_stream(aws, args)

@inline list_tags_for_stream(aws::AWSConfig, args) = AWSCore.Services.kinesisvideo(aws, "POST", "/listTagsForStream", args)

@inline list_tags_for_stream(args) = list_tags_for_stream(default_aws_config(), args)

"""
    using AWSSDK.KinesisVideo.tag_stream
    tag_stream([::AWSConfig], arguments::Dict)
    tag_stream([::AWSConfig]; Tags=, <keyword arguments>)

    using AWSCore.Services.kinesisvideo
    kinesisvideo([::AWSConfig], "POST", "/tagStream", arguments::Dict)
    kinesisvideo([::AWSConfig], "POST", "/tagStream", Tags=, <keyword arguments>)

# TagStream Operation

Adds one or more tags to a stream. A *tag* is a key-value pair (the value is optional) that you can define and assign to AWS resources. If you specify a tag that already exists, the tag value is replaced with the value that you specify in the request. For more information, see [Using Cost Allocation Tags](http://docs.aws.amazon.com/awsaccountbilling/latest/aboutv2/cost-alloc-tags.html) in the *AWS Billing and Cost Management User Guide*.

You must provide either the `StreamName` or the `StreamARN`.

This operation requires permission for the `KinesisVideo:TagStream` action.

Kinesis video streams support up to 50 tags.

# Arguments

## `StreamARN = ::String`
The Amazon Resource Name (ARN) of the resource that you want to add the tag or tags to.

## `StreamName = ::String`
The name of the stream that you want to add the tag or tags to.

## `Tags = ::Dict{String,String}` -- *Required*
A list of tags to associate with the specified stream. Each tag is a key-value pair (the value is optional).

# Returns

`TagStreamOutput`

# Exceptions

`ClientLimitExceededException`, `InvalidArgumentException`, `ResourceNotFoundException`, `NotAuthorizedException`, `InvalidResourceFormatException` or `TagsPerResourceExceededLimitException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/kinesisvideo-2017-09-30/TagStream)
"""
@inline tag_stream(aws::AWSConfig=default_aws_config(); args...) = tag_stream(aws, args)

@inline tag_stream(aws::AWSConfig, args) = AWSCore.Services.kinesisvideo(aws, "POST", "/tagStream", args)

@inline tag_stream(args) = tag_stream(default_aws_config(), args)

"""
    using AWSSDK.KinesisVideo.untag_stream
    untag_stream([::AWSConfig], arguments::Dict)
    untag_stream([::AWSConfig]; TagKeyList=, <keyword arguments>)

    using AWSCore.Services.kinesisvideo
    kinesisvideo([::AWSConfig], "POST", "/untagStream", arguments::Dict)
    kinesisvideo([::AWSConfig], "POST", "/untagStream", TagKeyList=, <keyword arguments>)

# UntagStream Operation

Removes one or more tags from a stream. In the request, specify only a tag key or keys; don't specify the value. If you specify a tag key that does not exist, it's ignored.

In the request, you must provide the `StreamName` or `StreamARN`.

# Arguments

## `StreamARN = ::String`
The Amazon Resource Name (ARN) of the stream that you want to remove tags from.

## `StreamName = ::String`
The name of the stream that you want to remove tags from.

## `TagKeyList = [::String, ...]` -- *Required*
A list of the keys of the tags that you want to remove.

# Returns

`UntagStreamOutput`

# Exceptions

`ClientLimitExceededException`, `InvalidArgumentException`, `ResourceNotFoundException`, `NotAuthorizedException` or `InvalidResourceFormatException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/kinesisvideo-2017-09-30/UntagStream)
"""
@inline untag_stream(aws::AWSConfig=default_aws_config(); args...) = untag_stream(aws, args)

@inline untag_stream(aws::AWSConfig, args) = AWSCore.Services.kinesisvideo(aws, "POST", "/untagStream", args)

@inline untag_stream(args) = untag_stream(default_aws_config(), args)

"""
    using AWSSDK.KinesisVideo.update_data_retention
    update_data_retention([::AWSConfig], arguments::Dict)
    update_data_retention([::AWSConfig]; CurrentVersion=, Operation=, DataRetentionChangeInHours=, <keyword arguments>)

    using AWSCore.Services.kinesisvideo
    kinesisvideo([::AWSConfig], "POST", "/updateDataRetention", arguments::Dict)
    kinesisvideo([::AWSConfig], "POST", "/updateDataRetention", CurrentVersion=, Operation=, DataRetentionChangeInHours=, <keyword arguments>)

# UpdateDataRetention Operation

Increases or decreases the stream's data retention period by the value that you specify. To indicate whether you want to increase or decrease the data retention period, specify the `Operation` parameter in the request body. In the request, you must specify either the `StreamName` or the `StreamARN`.

**Note**
> The retention period that you specify replaces the current value.

This operation requires permission for the `KinesisVideo:UpdateDataRetention` action.

Changing the data retention period affects the data in the stream as follows:

*   If the data retention period is increased, existing data is retained for the new retention period. For example, if the data retention period is increased from one hour to seven hours, all existing data is retained for seven hours.

*   If the data retention period is decreased, existing data is retained for the new retention period. For example, if the data retention period is decreased from seven hours to one hour, all existing data is retained for one hour, and any data older than one hour is deleted immediately.

# Arguments

## `StreamName = ::String`
The name of the stream whose retention period you want to change.

## `StreamARN = ::String`
The Amazon Resource Name (ARN) of the stream whose retention period you want to change.

## `CurrentVersion = ::String` -- *Required*
The version of the stream whose retention period you want to change. To get the version, call either the `DescribeStream` or the `ListStreams` API.

## `Operation = "INCREASE_DATA_RETENTION" or "DECREASE_DATA_RETENTION"` -- *Required*
Indicates whether you want to increase or decrease the retention period.

## `DataRetentionChangeInHours = ::Int` -- *Required*
The retention period, in hours. The value you specify replaces the current value.

# Returns

`UpdateDataRetentionOutput`

# Exceptions

`ClientLimitExceededException`, `InvalidArgumentException`, `ResourceNotFoundException`, `ResourceInUseException`, `NotAuthorizedException` or `VersionMismatchException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/kinesisvideo-2017-09-30/UpdateDataRetention)
"""
@inline update_data_retention(aws::AWSConfig=default_aws_config(); args...) = update_data_retention(aws, args)

@inline update_data_retention(aws::AWSConfig, args) = AWSCore.Services.kinesisvideo(aws, "POST", "/updateDataRetention", args)

@inline update_data_retention(args) = update_data_retention(default_aws_config(), args)

"""
    using AWSSDK.KinesisVideo.update_stream
    update_stream([::AWSConfig], arguments::Dict)
    update_stream([::AWSConfig]; CurrentVersion=, <keyword arguments>)

    using AWSCore.Services.kinesisvideo
    kinesisvideo([::AWSConfig], "POST", "/updateStream", arguments::Dict)
    kinesisvideo([::AWSConfig], "POST", "/updateStream", CurrentVersion=, <keyword arguments>)

# UpdateStream Operation

Updates stream metadata, such as the device name and media type.

You must provide the stream name or the Amazon Resource Name (ARN) of the stream.

To make sure that you have the latest version of the stream before updating it, you can specify the stream version. Kinesis Video Streams assigns a version to each stream. When you update a stream, Kinesis Video Streams assigns a new version number. To get the latest stream version, use the `DescribeStream` API.

`UpdateStream` is an asynchronous operation, and takes time to complete.

# Arguments

## `StreamName = ::String`
The name of the stream whose metadata you want to update.

The stream name is an identifier for the stream, and must be unique for each account and region.

## `StreamARN = ::String`
The ARN of the stream whose metadata you want to update.

## `CurrentVersion = ::String` -- *Required*
The version of the stream whose metadata you want to update.

## `DeviceName = ::String`
The name of the device that is writing to the stream.

**Note**
> In the current implementation, Kinesis Video Streams does not use this name.

## `MediaType = ::String`
The stream's media type. Use `MediaType` to specify the type of content that the stream contains to the consumers of the stream. For more information about media types, see [Media Types](http://www.iana.org/assignments/media-types/media-types.xhtml). If you choose to specify the `MediaType`, see [Naming Requirements](https://tools.ietf.org/html/rfc6838#section-4.2).

To play video on the console, you must specify the correct video type. For example, if the video in the stream is H.264, specify `video/h264` as the `MediaType`.

# Returns

`UpdateStreamOutput`

# Exceptions

`ClientLimitExceededException`, `InvalidArgumentException`, `ResourceNotFoundException`, `ResourceInUseException`, `NotAuthorizedException` or `VersionMismatchException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/kinesisvideo-2017-09-30/UpdateStream)
"""
@inline update_stream(aws::AWSConfig=default_aws_config(); args...) = update_stream(aws, args)

@inline update_stream(aws::AWSConfig, args) = AWSCore.Services.kinesisvideo(aws, "POST", "/updateStream", args)

@inline update_stream(args) = update_stream(default_aws_config(), args)

end # module KinesisVideo

#==============================================================================#
# End of file
#==============================================================================#


text:
using StanSample, Random, Statistics

stan_chris = "
data {
    int n_rows;
    int<lower=1> n_cols;
    matrix<lower=0>[n_rows,n_cols] x;
}

parameters {
   real mu;
}

model {
    mu ~ normal(0, 1);
}";

seed = 65445
Random.seed!(seed)

n_rows = 10
n_cols = 2
data = [(x = rand(n_rows, n_cols), n_rows = n_rows, n_cols = n_cols) for i in 1:4]

tmpdir = joinpath(@__DIR__, "tmp")
sm = SampleModel("chris", stan_chris, tmpdir);
rc = stan_sample(sm; data, seed);

if success(rc)
	df = read_samples(sm , :dataframe)
    #df |> display
    [mean(df.mu), std(df.mu)] |> display
end

tmpdir = joinpath(@__DIR__, "tmp")
sm = SampleModel("chris", stan_chris, tmpdir);
rc = stan_sample(sm; data=data[1:3], seed);

if success(rc)
    df = read_samples(sm , :dataframe)
    #df |> display
    [mean(df.mu), std(df.mu)] |> display
end


text:
Field1(inner,::Vector{LagrangianEntry{SimpleNL.HessianNull}}) = inner

Hessian11a(::HessianNull,::MyRef{Float64}) = HESSIAN_NULL
Hessian22a(::HessianNull,::HessianNull,::MyRef{Float64},::MyRef{Float64}) = HESSIAN_NULL
Hessian22a(::HessianNull,h2,::MyRef{Float64},ref2::MyRef{Float64}) = Hessian11a(h2,ref2)
Hessian22a(h1,::HessianNull,ref1::MyRef{Float64},::MyRef{Float64}) = Hessian11a(h1,ref1)

# Hessian11(::FZero1,h1,h11,::MyRef{Float64},ref::MyRef{Float64}) = Hessian11a(h1,ref)
# Hessian111(ddf22::FZero2,e,d,indexer) = Hessian11a(Hessian(e.e2,d.d1,indexer),ref(d))
# Hessian112(ddf11::FZero2,e,d,indexer) = Hessian11a(Hessian(e.e1,d.d1,indexer),ref(d))


text:
export
    takewhile, dropwhile,
    tmap,
    NestedIterator, DefaultNestedIterator, nested

#
# BEGIN
# http://slendermeans.org/julia-iterators.html
#

abstract type WhileIterator end

struct TakeWhile{I, F} <: WhileIterator
    xs::I
    cond::F
end

struct DropWhile{I, F} <: WhileIterator
    xs::I
    cond::F
end

takewhile(cond, xs) = TakeWhile(xs, cond)
dropwhile(cond, xs) = DropWhile(xs, cond)

Base.start(it::TakeWhile) = begin
    current_state = start(it.xs)
    done(it.xs, current_state) && return true, nothing, nothing

    current_value, next_state = next(it.xs, current_state)
    false, current_value, next_state
end
Base.start(it::DropWhile) = begin
    current_state = start(it.xs)
    while !done(it.xs, current_state)
        current_value, next_state = next(it.xs, current_state)
        !it.cond(current_value) && return false, current_value, next_state
        current_state = next_state
    end
    true, nothing, nothing
end

Base.next(it::WhileIterator, state) = begin
    _, current_value, next_state = state
    done(it.xs, next_state) && return current_value, (true, nothing, nothing)

    next_value, nextnext_state = next(it.xs, next_state)
    current_value, (false, next_value, nextnext_state)
end

Base.done(it::TakeWhile, state) = begin
    current_done, current_value, _ = state
    current_done || !it.cond(current_value)
end
Base.done(::DropWhile, state) = state[1]

Base.eltype(it::TakeWhile) = Base.eltype(typeof(it))
Base.eltype(it::DropWhile) = Base.eltype(typeof(it))
Base.eltype(::Type{TakeWhile{I,F}}) where {I,F} = Base.eltype(I)
Base.eltype(::Type{DropWhile{I,F}}) where {I,F} = Base.eltype(I)

Base.iteratorsize(::TakeWhile) = Base.SizeUnknown()
Base.iteratorsize(::DropWhile) = Base.SizeUnknown()

#
# END
# http://slendermeans.org/julia-iterators.html
#

struct TMap{T}
    mapfunc::Base.Callable
    xs::Vector{Any}
end

tmap(mapfunc, resulttype::Type, it1, its...) =
    TMap{resulttype}(mapfunc, Any[it1, its...])

Base.start(it::TMap) = map(start, it.xs)
Base.next(it::TMap, state) = begin
    next_result = map(next, it.xs, state)
    (
        it.mapfunc(map(x -> x[1], next_result)...),
        map(x -> x[2], next_result)
    )
end
Base.done(it::TMap, state) = any(map(done, it.xs, state))

Base.eltype(it::TMap) = Base.eltype(typeof(it))
Base.eltype(::Type{TMap{T}}) where {T} = T

Base.length(it::TMap) = minimum(map(length, it.xs))

abstract type NestedIterator end

struct DefaultNestedIterator{T} <: NestedIterator
    source::T
    level::Int
end

nested(s::T, l::Int = 1) where {T} = DefaultNestedIterator{T}(s, l)

create(::NestedIterator, ::Any) = []
combine(::NestedIterator, outer, inner) = "$outer,$inner"

Base.start(it::NestedIterator) = nextstate(it, (Void, Void, Void))

Base.next(it::NestedIterator, state) = begin
    outerit = it.source
    outerstate, innerit, innerstate = state

    outervalue, _ = next(outerit, outerstate)
    innervalue, _ = next(innerit, innerstate)

    value = combine(it, outervalue, innervalue)
    state = nextstate(it, (outerstate, innerit, innerstate))

    value, state
end

Base.done(it::NestedIterator, state) = begin
    outerstate, _, _ = state
    done(it.source, outerstate)
end

nextstate(it::NestedIterator, state) = begin
    outerit = it.source
    outerstate, innerit, innerstate = state

    if outerstate == Void
        outerstate = start(outerit)
    end

    done(outerit, outerstate) && return (outerstate, innerit, innerstate)

    if innerit == Void
        innerit, innerstate = createinner(it, outerstate)
    else
        _, innerstate = next(innerit, innerstate)
    end

    while done(innerit, innerstate)
        _, outerstate = next(outerit, outerstate)
        done(outerit, outerstate) && break
        innerit, innerstate = createinner(it, outerstate)
    end

    (outerstate, innerit, innerstate)
end

createinner(it::NestedIterator, outerstate) = begin
    outervalue, _ = next(it.source, outerstate)
    innerit = create(it, outervalue)
    innerstate = start(innerit)

    innerit, innerstate
end

Base.iteratorsize(::NestedIterator) = Base.SizeUnknown()


text:
using Healpix, Libsharp
using PyCall
using Random, Statistics
using Printf, Tables, CSV
# %% Simulation parameters
nrz = 1500 # How many realisations?
Random.seed!(5147) # Initial random number seed. Useful if you need reproducible sequence
rsim = 0 # Tensor-to-scalar ratio used for the simulation
Alens = 1 # Lensing power spectrum amplitude (Alens = 1 for the fiducial)

# %% Specification of the experiments
# Reference: Simons Observatory Collaboration, JCAP, 02, 056 (2019), Table 1.
# Reference: CCAT-prime Collaboration, J.Low Temp.Phys., 199, 1089 (2020), Table 1.
ν = [27, 39, 93, 145, 225, 280, 350, 410, 850] # in GHz
nν, nνSO = length(ν), 6
FWHM = [91, 63, 30, 17, 11, 9, 0.58, 0.5, 0.23] # in arcmin
σ = FWHM * π / 10800 / √(8 * log(2)) # in radians
uKarcmin = [35, 21, 2.6, 3.3, 6.3, 16, 105, 372, 5.7e5] # in μK arcmin (for temperature; x√2 for pol)
lknee = [30, 30, 50, 50, 70, 100, 700, 700, 700]
αknee = [-2.4, -2.4, -2.5, -3, -3, -3, -1.4, -1.4, -1.4]
# temperature beam
bTl(ℓ, σb) = exp(-ℓ * (ℓ + 1) * σb^2 / 2)
# polarisation beam, Eq.(5.8) of Ng & Liu, Int.J.Mod.Phys.D, 8, 61 (1999)
bPl(ℓ, σb) = ifelse(ℓ >= 2, exp(-(ℓ * (ℓ + 1) - 4) * σb^2 / 2), 0)

# %% Read in the hits map and calculate weights for inhomogeneous noise
nhitsfile = "nhits_SAT_r7.FITS"
nhits = readMapFromFITS(nhitsfile, 1, Float64)
nside = nhits.resolution.nside
nhits /= maximum(nhits)
fsky = mean(nhits)^2 / mean(nhits .^ 2)
K = √(mean(nhits) / fsky)
weight = zeros(12 * nside^2)
for ip = 1:12*nside^2
    if nhits[ip] > 0
        weight[ip] = K / √nhits[ip]
    end
end

# %% Compute scalar and tensor power spectra using CLASS
lmax, mmax = 3 * nside - 1, 3 * nside - 1
include("compute_cl_class.jl")

# %% Basic parameters for spherical harmonics transform
geom_info = make_healpix_geom_info(nside, 1)
alm_info = make_triangular_alm_info(lmax, mmax, 1)
elm, blm = Alm{ComplexF64}(lmax, mmax), Alm{ComplexF64}(lmax, mmax)

# %% Read in and smooth foreground maps generated by "d1s1" model of Python Sky Model (PySM), in equatorial coordinates
# Reference: Thorne et al., MNRAS, 469, 2821 (2017), https://github.com/healpy/pysm
f_q, f_u = [], []
for iν = 1:nν
    filename = @sprintf("map_d1s1_equ_%03dghz_r7_uKcmb_qu.fits", ν[iν])
    q, u = (readMapFromFITS(filename, ic, Float64) for ic = 1:2)
    sharp_execute!(
        SHARP_MAP2ALM,
        2,
        [elm.alm, blm.alm],
        [q.pixels, u.pixels],
        geom_info,
        alm_info,
        SHARP_DP,
    )
    for l = 0:lmax
        for m = 0:l
            ilm = alm_index(alm_info, l, m) + 1
            elm.alm[ilm] *= bPl(l, σ[iν])
            blm.alm[ilm] *= bPl(l, σ[iν])
        end
    end
    sharp_execute!(
        SHARP_ALM2MAP,
        2,
        [elm.alm, blm.alm],
        [q.pixels, u.pixels],
        geom_info,
        alm_info,
        SHARP_DP,
    )
    push!(f_q, q)
    push!(f_u, u)
end

# %% Setup NaMaster for the power spectrum analysis on a partial sky
Δℓ = 10 # multipole binning size
maskfile = "mask_apodized_r7.fits"
include("setup_namaster.jl")
# The theory power spectrum must be binned into bandpowers in the same manner the data has.
# Generate an NmtWorkspace object that we use to compute and store the mode coupling matrix.
# Note that this matrix depends only on the masks of the two fields to correlate, but not on the maps themselves.
w = nmt.NmtWorkspace()
f = nmt.NmtField(mask, [f_q[1], f_u[1]], purify_b = true)
w.compute_coupling_matrix(f, f, b)
cls_th = [cls["ee"], zero(cls["ee"]), zero(cls["bb"]), cls["bb"]] * Tcmb^2
cls_th_binned = w.decouple_cell(w.couple_cell(cls_th))
clt_th = [clt["ee"], zero(clt["ee"]), zero(clt["bb"]), clt["bb"]] * Tcmb^2
clt_th_binned = w.decouple_cell(w.couple_cell(clt_th))
CSV.write(
    "scalar_eebb_binned.csv",
    Tables.table([ell_eff cls_th_binned[1, :] cls_th_binned[4, :]]),
    header = ["leff", "ee", "bb"],
)
CSV.write(
    "tensor_eebb_binned.csv",
    Tables.table([ell_eff clt_th_binned[1, :] clt_th_binned[4, :]]),
    header = ["leff", "ee", "bb"],
)

# %% Compute covariance matrices of the foreground EE and BB
# Used to compute the power spectra of residual foreground for checking the results,
# but this information cannot be used for the real data analysis.
f3 = []
for iν = 1:nν
    push!(
        f3,
        nmt.NmtField(
            mask,
            [f_q[iν], f_u[iν]],
            purify_b = true,
            beam = bPl.(0:lmax, σ[iν]),
        ),
    )
end
ce3, cb3 = zeros(nν, nν, nbands), zeros(nν, nν, nbands) # foreground
for iν = 1:nν, jν = iν:nν
    w.compute_coupling_matrix(f3[iν], f3[jν], b)
    cl3 = compute_master(f3[iν], f3[jν], w)
    ce3[iν, jν, :], cb3[iν, jν, :] = cl3[1, :], cl3[4, :]
end
# Write out to binary files
open("cov_bb_foreg.dat", "w") do io
    write(io, cb3)
end
open("cov_ee_foreg.dat", "w") do io
    write(io, ce3)
end

# %% Loop over realisations
for irz = 1:nrz
    @show irz
    # Create spherical harmonics coefficients of the CMB
    for l = 2:lmax
        ee = cls_th[1][l+1] + (rsim / rclass) * clt_th[1][l+1]
        bb = Alens * cls_th[4][l+1] + (rsim / rclass) * clt_th[4][l+1]
        ilm = alm_index(alm_info, l, 0) + 1
        elm.alm[ilm] = √ee * randn(Float64)
        blm.alm[ilm] = √bb * randn(Float64)
        for m = 1:l
            ilm = alm_index(alm_info, l, m) + 1
            elm.alm[ilm] = √ee * randn(ComplexF64)
            blm.alm[ilm] = √bb * randn(ComplexF64)
        end
    end
    # Smooth CMB and create spherical harmonics coefficients of the noise
    celm, cblm = Alm{ComplexF64}(lmax, mmax), Alm{ComplexF64}(lmax, mmax)
    nelm, nblm = Alm{ComplexF64}(lmax, mmax), Alm{ComplexF64}(lmax, mmax)
    c_q, c_u = HealpixMap{Float64,RingOrder}(nside), HealpixMap{Float64,RingOrder}(nside)
    n_q, n_u = HealpixMap{Float64,RingOrder}(nside), HealpixMap{Float64,RingOrder}(nside)
    f1, f2 = [], [] # List of NaMaster fields
    for iν = 1:nν
        for l = 0:lmax
            ee =
                2 *
                (uKarcmin[iν] * π / 10800)^2 *
                (1 + (l / lknee[iν])^αknee[iν])
            bb =
                2 *
                (uKarcmin[iν] * π / 10800)^2 *
                (1 + (l / lknee[iν])^αknee[iν])
            if l >= 2
                ilm = alm_index(alm_info, l, 0) + 1
                celm.alm[ilm] = elm.alm[ilm] * bPl(l, σ[iν])
                cblm.alm[ilm] = blm.alm[ilm] * bPl(l, σ[iν])
                nelm.alm[ilm] = √ee * randn(Float64)
                nblm.alm[ilm] = √bb * randn(Float64)
                for m = 1:l
                    ilm = alm_index(alm_info, l, m) + 1
                    celm.alm[ilm] = elm.alm[ilm] * bPl(l, σ[iν])
                    cblm.alm[ilm] = blm.alm[ilm] * bPl(l, σ[iν])
                    nelm.alm[ilm] = √ee * randn(ComplexF64)
                    nblm.alm[ilm] = √bb * randn(ComplexF64)
                end
            else
                for m = 0:l
                    ilm = alm_index(alm_info, l, m) + 1
                    celm.alm[ilm], cblm.alm[ilm] = 0, 0
                    nelm.alm[ilm], nblm.alm[ilm] = 0, 0
                end
            end
        end
        sharp_execute!(
            SHARP_ALM2MAP,
            2,
            [celm.alm, cblm.alm],
            [c_q.pixels, c_u.pixels],
            geom_info,
            alm_info,
            SHARP_DP,
        )
        sharp_execute!(
            SHARP_ALM2MAP,
            2,
            [nelm.alm, nblm.alm],
            [n_q.pixels, n_u.pixels],
            geom_info,
            alm_info,
            SHARP_DP,
        )
        # Inhomogenous noise for SO-SAT and homogeneous noise for CCAT-prime
        if iν <= nνSO
            n_q .*= weight
            n_u .*= weight
        end
        ## Create NaMaster fields for the pseudo-Cℓ
        # f1: foreground + CMB + noise
        push!(
            f1,
            nmt.NmtField(
                mask,
                [f_q[iν] + c_q + n_q, f_u[iν] + c_u + n_u],
                purify_b = true,
                beam = bPl.(0:lmax, σ[iν]),
            ),
        )
        # f2: noise
        push!(
            f2,
            nmt.NmtField(
                mask,
                [n_q, n_u],
                purify_b = true,
                beam = bPl.(0:lmax, σ[iν]),
            ),
        )
    end
    ## Compute covariance matrices of EE and BB
    ce1, cb1 = zeros(nν, nν, nbands), zeros(nν, nν, nbands) # total
    ce2, cb2 = zeros(nν, nν, nbands), zeros(nν, nν, nbands) # noise
    for iν = 1:nν, jν = iν:nν
        w.compute_coupling_matrix(f1[iν], f1[jν], b)
        cl1 = compute_master(f1[iν], f1[jν], w)
        cl2 = compute_master(f2[iν], f2[jν], w)
        ce1[iν, jν, :], cb1[iν, jν, :] = cl1[1, :], cl1[4, :]
        ce2[iν, jν, :], cb2[iν, jν, :] = cl2[1, :], cl2[4, :]
    end
    # Write out to binary files
    open(@sprintf("cov_bb_total_irz%04d.dat", irz), "w") do io
        write(io, cb1)
    end
    open(@sprintf("cov_bb_noise_irz%04d.dat", irz), "w") do io
        write(io, cb2)
    end
    open(@sprintf("cov_ee_total_irz%04d.dat", irz), "w") do io
        write(io, ce1)
    end
    open(@sprintf("cov_ee_noise_irz%04d.dat", irz), "w") do io
        write(io, ce2)
    end
end


text:
@testset "Testing robustlq" begin
    rough_kwargs = Dict(:atol => 1e-4, :rtol => 1e-4)

    # set up
    a_0     = 100
    a_1     = 0.5
    ρ     = 0.9
    sigma_d = 0.05
    β    = 0.95
    c       = 2
    γ   = 50.0
    θ   = 0.002
    ac      = (a_0 - c) / 2.0

    R = [0   ac    0
         ac  -a_1  0.5
         0.  0.5   0]

    R = -R
    Q = γ / 2

    A = [1. 0. 0.
         0. 1. 0.
         0. 0. ρ]
    B = [0.0 1.0 0.0]'
    C = [0.0 0.0 sigma_d]'

    rblq = RBLQ(Q, R, A, B, C, β, θ)
    lq = LQ(Q, R, A, B, C, β)

    Fr, Kr, Pr = robust_rule(rblq)

    # test stuff
    @testset "test robust vs simple" begin
        Fs, Ks, Ps = robust_rule_simple(rblq, Pr; tol=1e-12)

        @test isapprox(Fr, Fs; rough_kwargs...)
        @test isapprox(Kr, Ks; rough_kwargs...)
        @test isapprox(Pr, Ps; rough_kwargs...)
    end

    @testset "test f2k and k2f" begin
        K_f2k, P_f2k = F_to_K(rblq, Fr)
        F_k2f, P_k2f = K_to_F(rblq, Kr)

        @test isapprox(K_f2k, Kr; rough_kwargs...)
        @test isapprox(F_k2f, Fr; rough_kwargs...)
        @test isapprox(P_f2k, P_k2f; rough_kwargs...)
    end

    @testset "test evaluate F" begin
        Kf, Pf, df, Of, of =  evaluate_F(rblq, Fr)

        @test isapprox(Pf, Pr; rough_kwargs...)
        @test isapprox(Kf, Kr; rough_kwargs...)
    end

    @testset "test no run-time error in robust_rule_simple" begin
        # this will just print out a warning
        robust_rule_simple(rblq, ones(Pr); tol=eps(), max_iter=1)
    end
end  # @testset


text:
"""
    c1(img, HΦ)

Calculate the average value inside the function

# Arguments
- img: input image
- HΦ: Heaviside of the level set function `Φ`
"""
function c1(img, HΦ)
    ∫H𝚽 = sum(HΦ)
    ∫u₀H𝚽 = sum(img .* HΦ)
    return c₁ = ∫u₀H𝚽 / ∫H𝚽
end

"""
    c2(img, HΦ)

Calculate the average value outside the function

# Arguments
- img: input image
- HΦ: Heaviside of the level set function `Φ`
"""
function c2(img, HΦ)
    H𝚽ⁱ = 1 .- HΦ
    ∫H𝚽ⁱ = sum(H𝚽ⁱ)
    ∫u₀H𝚽ⁱ = sum(img .* H𝚽ⁱ)
    return c₂ = ∫u₀H𝚽ⁱ / ∫H𝚽ⁱ
end

"""
    δₕ(x, h=1)

Modified Dirac Delta function

# Arguments
- x: input
"""
δₕ(x, h=1) = h ./ (h .^ 2 .+ x .^ 2)

# Equations for `calculate_variation`
raw"""
```math
\begin{aligned}
	\phi_{i+} 	&= \phi_{i, j+1}^n - \phi_{i, j}^n \\
	\phi_{i-} 	&= \phi_{i, j}^n - \phi_{i, j - 1}^n \\
	\phi_{i} 	&= \frac{\phi_{i, j+1}^n - \phi_{i, j-1}^n}{2} \\
	\phi_{j+} 	&= \phi_{i+1, j}^n - \phi_{i, j}^n \\
	\phi_{j-} 	&= \phi_{i, j}^n - \phi_{i-1, j}^n \\
	\phi_{j} 	&= \frac{\phi_{i+1, j}^n - \phi_{i-1, j}^n}{2} \\
\end{aligned}
\tag{9a}
```
```math
\begin{aligned}
	C_1 &= \frac{1}{\sqrt{\phi_{i+}^2 - \phi_{j}^2}} \\
	C_2 &= \frac{1}{\sqrt{\phi_{i-}^2 - \phi_{j}^2}} \\
	C_3 &= \frac{1}{\sqrt{\phi_{i}^2 - \phi_{j+}}} \\
	C_4 &= \frac{1}{\sqrt{\phi_{i}^2 - \phi_{j+}}}
\end{aligned}
\tag{9b}
```
```math
	C = C_1 * \phi_{i, j+1}^n + C_2 * \phi_{i, j-1}^n + C_3 * \phi_{i+1, j}^n + C_4 * \phi_{i-1, j}^n
\tag{9c}
```
```math
\begin{aligned}
D = -\lambda_{1} (u_{0_{i,j}} - c_1*(\phi_{i,j}^n))^2 + \lambda_{2} (u_{0_{i,j}} - c_2 * (\phi_{i,j}^n))^2
\end{aligned}
\tag{9d}
```
```math
\begin{aligned}
	\phi_{i,j}^{n+1} = \frac{\phi_{i,j}^n + \Delta t * \delta(\phi_{i,j}^n) * (\mu * C + D)}{ 1 + \mu * \Delta t * \delta(\phi_{i,j}^n) * (C_1 + C_2 + C_3 + C_4)}
\end{aligned}
\tag{9e}
```
"""
"""
	calculate_variation(img, 𝚽ⁿ, μ, λ₁, λ₂, Δt)

Calculates the variation of level set `𝚽ⁿ` based on algorithm parameters.
Solves the set of equations from the given paper below:

# Arguments
- img: input image
- 𝚽ⁿ: level set function
- μ: "edge length" weight parameter. Higher `μ` values will
	produce a 'round' edge, while values closer to zero will
	detect smaller objects.
- λ₁: "difference from average" weight parameter for the output
	region with value 'true'. If it is lower than `λ₂`, this
	region will have a larger range of values than the other.
- λ₂: "difference from average" weight parameter for the output
	region with value 'false'. If it is lower than `λ₁`, this
	region will have a larger range of values than the other.
- Δt: A multiplication factor applied at calculations for each step,
	serves to accelerate the algorithm. While higher values may
	speed up the algorithm, they may also lead to convergence
	problems.

# Citation
'The Chan-Vese Algorithm' [Rami Cohen](https://arxiv.org/pdf/1107.2782.pdf)
"""
function calculate_variation(img, 𝚽ⁿ, μ, λ₁, λ₂, Δt)
    ϵ = 1e-16
    𝚽⁺ = padarray(𝚽ⁿ, Pad(1, 1))

    # 9a
    𝚽ᵢ₊ = 𝚽⁺[1:(end - 1), 2:end] .- 𝚽⁺[1:(end - 1), 1:(end - 1)]
    𝚽ᵢ₋ = 𝚽⁺[1:(end - 1), 1:(end - 1)] .- 𝚽⁺[1:(end - 1), 0:(end - 2)]
    𝚽ᵢ = (𝚽⁺[1:(end - 1), 2:end] .- 𝚽⁺[1:(end - 1), 0:(end - 2)]) / 2

    𝚽ⱼ₊ = 𝚽⁺[2:end, 1:(end - 1)] - 𝚽⁺[1:(end - 1), 1:(end - 1)]
    𝚽ⱼ₋ = 𝚽⁺[1:(end - 1), 1:(end - 1)] - 𝚽⁺[0:(end - 2), 1:(end - 1)]
    𝚽ⱼ = (𝚽⁺[2:end, 1:(end - 1)] - 𝚽⁺[0:(end - 2), 1:(end - 1)]) / 2

    # 9b
    C₁ = 1 ./ sqrt.(ϵ .+ 𝚽ᵢ₊ .^ 2 .+ 𝚽ⱼ .^ 2)
    C₂ = 1 ./ sqrt.(ϵ .+ 𝚽ᵢ₋ .^ 2 .+ 𝚽ⱼ .^ 2)
    C₃ = 1 ./ sqrt.(ϵ .+ 𝚽ᵢ .^ 2 .+ 𝚽ⱼ₊ .^ 2)
    C₄ = 1 ./ sqrt.(ϵ .+ 𝚽ᵢ .^ 2 .+ 𝚽ⱼ₋ .^ 2)

    # 9c
    C =
        𝚽⁺[1:(end - 1), 2:end] .* C₁ .+ 𝚽⁺[1:(end - 1), 0:(end - 2)] .* C₂ .+
        𝚽⁺[2:end, 1:(end - 1)] .* C₃ .+ 𝚽⁺[0:(end - 2), 1:(end - 1)] .* C₄

    H𝚽 = 1 .* (𝚽ⁿ .> 0)
    c₁, c₂ = c1(img, H𝚽), c2(img, H𝚽)

    # 9d
    D = -λ₁ .* (img .- c₁) .^ 2 .+ λ₂ .* (img .- c₂) .^ 2

    # 9
    return 𝚽ⁿ⁺¹ =
        (𝚽ⁿ .+ Δt .* δₕ(𝚽ⁿ) .* (μ * C .+ D)) ./
        (1 .+ μ .* Δt .* δₕ(𝚽ⁿ) .* (C₁ .+ C₂ .+ C₃ .+ C₄))
end

# Equations for `calculate_reinitial`
raw"""
```math
\begin{aligned}
	a &= (\Psi_{i,j} - \Psi_{i-1,j}) \\
	b &= (\Psi_{i,j} - \Psi_{i+1,j}) \\
	c &= (\Psi_{i,j} - \Psi_{i,j-1}) \\
	d &= (\Psi_{i,j} - \Psi_{i,j+1})
\end{aligned}
\tag{12a}
```
```math
\begin{aligned}
	G(\Psi_{i,j}^{n}) &=
    \begin{cases}
		\sqrt{max((a^+)^2, (b^-)^2) + max((c^+)^2, (d^-)^2)} - 1,
			\ &\Phi(x_i, j_i, t) > 0 \\
		\sqrt{max((a^-)^2, (b^+)^2) + max((c^-)^2, (d^+)^2)} - 1,
			\ &\Phi(x_i, j_i, t) < 0 \\
      0 \ &\text{otherwise}
    \end{cases}
\end{aligned}
\tag{12b}
```
```math
\begin{aligned}
	sign(\Phi(x,y,t)) &= \frac{\Phi(x,y,t)}{\sqrt{\Phi(x,y,t)^2}}
\end{aligned}
\tag{11a}
```
```math
\begin{aligned}

	\Psi_{i,j}^{n+1} &= \Psi_{i,j}^{n} - \Delta \tau * sign(\Phi(x,y,t) * G(\Psi_{i,j}^{n})
\end{aligned}
\tag{11b}
```
"""
"""
	calculate_reinitial(𝚽, Δt)

Calculates the way in which to reinitialize the level set `𝚽` based on
algorithm parameters. Solves the set of equations from the given paper below:

# Arguments
- 𝚽: level set function
- Δt: A multiplication factor applied at calculations for each step,
	serves to accelerate the algorithm. While higher values may
	speed up the algorithm, they may also lead to convergence
	problems.

# Citation
'The Chan-Vese Algorithm' [Rami Cohen](https://arxiv.org/pdf/1107.2782.pdf)
"""
function calculate_reinitial(𝚽, Δt)
    ϵ = 1e-8
    𝚽⁺ = padarray(𝚽, Pad(1, 1))

    # 12a
    a = 𝚽⁺[1:(end - 1), 1:(end - 1)] - 𝚽⁺[1:(end - 1), 0:(end - 2)]
    b = 𝚽⁺[1:(end - 1), 2:end] - 𝚽⁺[1:(end - 1), 1:(end - 1)]
    c = 𝚽⁺[1:(end - 1), 1:(end - 1)] - 𝚽⁺[0:(end - 2), 1:(end - 1)]
    d = 𝚽⁺[2:end, 1:(end - 1)] - 𝚽⁺[1:(end - 1), 1:(end - 1)]

    a⁺ = max(a, 0)
    a⁻ = min(a, 0)
    b⁺ = max(b, 0)
    b⁻ = min(b, 0)
    c⁺ = max(c, 0)
    c⁻ = min(c, 0)
    d⁺ = max(d, 0)
    d⁻ = min(d, 0)

    G = zeros(size(𝚽))
    index⁺ = 𝚽 .> 0
    index⁻ = 𝚽 .< 0

    # 12b
    G =
        (sqrt.(max.(a⁺ .^ 2, b⁻ .^ 2) + max.(c⁺ .^ 2, d⁻ .^ 2)) .- 1) .* index⁺ .+
        (sqrt.(max.(a⁻ .^ 2, b⁺ .^ 2) .+ max.(c⁻ .^ 2, d⁺ .^ 2)) .- 1) .* index⁻

    # 11a
    sign𝚽 = 𝚽 ./ sqrt.(𝚽 .^ 2 .+ ϵ)

    # 11b
    return 𝚿 = 𝚽 .- Δt .* sign𝚽 .* G
end

function reinitialize(𝚽, Δt, max_reiter=5)
    iter = 0
    while iter < max_reiter
        𝚽 = calculate_reinitial(𝚽, Δt)
        iter += 1
    end

    return 𝚽
end

"""
	classical_chan_vese(
		img; μ=0.25, λ₁=1.0, λ₂=1.0, tol=1e-3, max_iter=500, Δt=0.5, reinitial_flag=false
	)

Chan-Vese segmentation algorithm. Active contour model by evolving a level set.
Can be used to segment objects without clearly defined boundaries.

# Arguments
- img: input image
- μ: "edge length" weight parameter. Higher `μ` values will
	produce a 'round' edge, while values closer to zero will
	detect smaller objects.
- λ₁: "difference from average" weight parameter for the output
	region with value 'true'. If it is lower than `λ₂`, this
	region will have a larger range of values than the other.
- λ₂: "difference from average" weight parameter for the output
	region with value 'false'. If it is lower than `λ₁`, this
	region will have a larger range of values than the other.
- tol: Level set variation tolerance between iterations. If the
	L2 norm difference between the level sets of successive
	iterations normalized by the area of the image is below this
	value, the algorithm will assume that the solution was
	reached.
- max_iter: Maximum number of iterations allowed before the algorithm
	interrupts itself.
- Δt: A multiplication factor applied at calculations for each step,
	serves to accelerate the algorithm. While higher values may
	speed up the algorithm, they may also lead to convergence
	problems.
- reinitial_flag: not sure what this does yet

# Citation
'Active contours without edges' [Chan, Vese](10.1109/83.902291)

'Chan–Vese Segmentation' [Pascal Getreuer](https://www.ipol.im/pub/art/2012/g-cv/article.pdf)

'The Chan-Vese Algorithm' [Rami Cohen](https://arxiv.org/pdf/1107.2782.pdf)
"""
function classical_chan_vese(
    img; μ=0.25, λ₁=1.0, λ₂=1.0, tol=1e-3, max_iter=500, Δt=0.5, reinitial_flag=false
)
    iter = 0
    𝚽ⁿ = checkerboard(size(img))
    δ = tol + 1
    img .= img .- minimum(img)

    if maximum(img) != 0
        img = img ./ maximum(img)
    end

    while (δ > tol) & (iter < max_iter)
        𝚽ⁿ⁺¹ = calculate_variation(img, 𝚽ⁿ, μ, λ₁, λ₂, Δt)
        δ = sqrt(meanfinite((𝚽ⁿ⁺¹ .- 𝚽ⁿ) .^ 2, (1, 2))[1])
        if reinitial_flag
            𝚽ⁿ .= reinitialize(𝚽ⁿ⁺¹, Δt)
        else
            𝚽ⁿ .= 𝚽ⁿ⁺¹
        end

        iter += 1
    end

    return 𝚽ⁿ, iter
end


text:
using SimpleDamage ; import SimpleDamage as SD
using NamedArrays
import DamagedShearBand as DSB
using Test

pp =   [8.4e9, 0.8, 0.55, 2e-3, 0.4, 20.0, 1.014e6, 0.001, 0.4]
pp_names = [:G,   :γ,   :μ,   :a,  :D₀, :n,   :K₁c,   :l̇₀,   :Dᵢ]
pn = NamedArray(pp,pp_names)
r = DSB.Rheology(μ=pn[:μ], ψ=45, a=pn[:a], D₀=pn[:D₀], n=pn[:n], K₁c=pn[:K₁c], l̇₀=pn[:l̇₀])
pc = 1e6
ϵ̇ = 1e-5
s = 1e8
C = 0.0
τy_const = 2e6

plas_coulomb_minthres = Plasticity(MinViscosityThreshold(nothing),
                                    CoulombYieldStress(;μ=pn[:μ],C)
                                    )
plas_const_minthres = Plasticity(MinViscosityThreshold(nothing),
                                    ConstantYieldStress(τy_const)
                                    )
damp = PrincipalKICharlesLaw(r)
dami = InvariantsKICharlesLaw(r)
sd = SimpleCharlesLaw()

lw = LinearWeakening(pn[:γ])

elas = IncompressibleElasticity(pn[:G])

cmp = ConstitutiveModel(;weakening=lw, damage=damp, elasticity=elas, plasticity=plas_coulomb_minthres)
cmi = ConstitutiveModel(;weakening=lw, damage=dami, elasticity=elas, plasticity=plas_coulomb_minthres)

csr = ConstantStrainRate(ϵ̇)
cs =ConstantStress(s)

loading2D = TriaxialSetup(; geom=Geom2D(), control=csr, pc)
loading3D = TriaxialSetup(; geom=Geom3D(), control=csr, pc)
creep2D = TriaxialSetup(; geom=Geom2D(), control=cs, pc)
creep3D = TriaxialSetup(; geom=Geom3D(), control=cs, pc)

mpl2D = Model(;setup=loading2D, cm=cmp)
mil2D = Model(;setup=loading2D, cm=cmi)
mil3D = Model(;setup=loading3D, cm=cmi)

@testset "helper functions" begin
    @test SD.confining_pressure(loading2D) == SD.confining_pressure(mpl2D) == pc

    @test SD.plastic_yield_stress(plas_coulomb_minthres,0) == plas_coulomb_minthres.σy.C * cos(atan(plas_coulomb_minthres.σy.μ))
    @test SD.plastic_yield_stress(plas_const_minthres,0) == plas_const_minthres.σy.val

    Di_less = r.D₀ - 0.1
    Di_more = r.D₀ + 0.1
    @test init_variables(damp,Di_less)[2] == damp.r.D₀
    @test init_variables(damp,Di_more)[2] == Di_more
end

@testset "stress invariants" begin
    @test SD.stress_invariants(0, loading2D) == (0, -loading2D.pc)
    @test SD.stress_invariants(0, loading3D) == (0, -loading3D.pc)
    @test SD.stress_invariants(-1, loading2D) == (1, -1 - loading2D.pc)
    @test SD.stress_invariants(-1, loading3D) == (sqrt(3)/2, -0.5 - loading3D.pc)
end

@testset "strain invariants" begin
    @test SD.ϵII_func(1, mil2D)  == 1
    @test SD.ϵII_func(1, mil3D) == sqrt(3)/2
end

@testset "weakening functions and derivatives" begin
    # simple damage
    let cm = ConstitutiveModel(;weakening=lw, damage=sd, elasticity=elas, plasticity=plas_coulomb_minthres)
        γ = cm.weakening.γ
        @test SD.weakening_func(0,cm) == 1
        @test SD.weakening_func(1,cm) == γ
        @test SD.weakening_derivative(0,cm) == SD.weakening_derivative(1,cm) == γ - 1
    end

    # micromechanical damage
    let cm = ConstitutiveModel(;weakening=lw, damage=damp, elasticity=elas, plasticity=plas_coulomb_minthres)
        γ = cm.weakening.γ
        D₀ = cm.damage.r.D₀
        @test SD.weakening_func(D₀,cm) == 1
        @test SD.weakening_func(1,cm) == γ
        @test SD.weakening_derivative(0,cm) == SD.weakening_derivative(1,cm) == (γ-1)/(1-D₀)
    end

end

@testset "axial stress proxies conversions" begin

    let geom = Geom2D()
        # pc = 0
        let sₐₓ = -1, pc = 0
            setup = TriaxialSetup(; geom, control=csr, pc)
            σₐₓtrue = -2
            Δσtrue = σₐₓtrue

            @test SD.sₐₓ2σₐₓ(setup,sₐₓ) == σₐₓtrue
            @test SD.σₐₓ2sₐₓ(setup,σₐₓtrue) == -1
            @test SD.sₐₓ2Δσ(setup,sₐₓ) == Δσtrue
            @test SD.Δσ2sₐₓ(setup,Δσtrue) == sₐₓ
        end
        # pc = 2
        let sₐₓ = -1, pc = 2
            setup = TriaxialSetup(; geom, control=csr, pc)
            σₐₓtrue = -4
            Δσtrue = -2

            @test SD.sₐₓ2σₐₓ(setup,sₐₓ) == σₐₓtrue
            @test SD.σₐₓ2sₐₓ(setup,σₐₓtrue) == sₐₓ
            @test SD.sₐₓ2Δσ(setup,sₐₓ) == Δσtrue
            @test SD.Δσ2sₐₓ(setup,Δσtrue) == sₐₓ
        end
    end

    let geom = Geom3D()
        # pc = 0
        let sₐₓ = -1, pc = 0
            setup = TriaxialSetup(; geom, control=csr, pc)
            σₐₓtrue = -3/2
            Δσtrue = σₐₓtrue
            @test SD.sₐₓ2σₐₓ(setup,sₐₓ) == σₐₓtrue
            @test SD.σₐₓ2sₐₓ(setup,σₐₓtrue) == -1
            @test SD.sₐₓ2Δσ(setup,sₐₓ) == Δσtrue
            @test SD.Δσ2sₐₓ(setup,Δσtrue) == sₐₓ
        end
        # pc = 1
        let sₐₓ = -1, pc = 2
            setup = TriaxialSetup(; geom, control=csr, pc)
            σₐₓtrue = -3/2 -2
            Δσtrue = -3/2
            @test SD.sₐₓ2σₐₓ(setup,sₐₓ) == σₐₓtrue
            @test SD.σₐₓ2sₐₓ(setup,σₐₓtrue) == sₐₓ
            @test SD.sₐₓ2Δσ(setup,sₐₓ) == Δσtrue
            @test SD.Δσ2sₐₓ(setup,Δσtrue) == sₐₓ
        end
    end

end

@testset "KI" begin
    # TEST KI WITH INVARIANTS
    let m = Model(;setup=loading2D, cm=cmi)
        a = m.cm.damage.r.a
        μ = m.cm.damage.r.μ
        D₀ = m.cm.damage.r.D₀
        pc = m.setup.pc
        ψ = m.cm.damage.r.ψ

        let s = -μ*pc/(1-μ)#- μ*pc/(0.5*(sqrt(3) - μ))
            τ, σₘ = SD.stress_invariants(s, m)
            # at D = D0 and τ/σ=μ -> KI = O
            @test DSB.compute_KI(r,σₘ,τ,D₀) ≈ 0
            @test SD.get_KI_from_s(s, D₀, m) ≈ 0
        end
    end

    let m = Model(;setup=loading3D, cm=cmi)
        a = m.cm.damage.r.a
        μ = m.cm.damage.r.μ
        D₀ = m.cm.damage.r.D₀
        pc = m.setup.pc
        ψ = m.cm.damage.r.ψ

        let s = - μ*pc/(0.5*(sqrt(3) - μ))
            τ, σₘ = SD.stress_invariants(s, m)
            # at D = D0 and τ/σ=μ -> KI = O
            @test isapprox(DSB.compute_KI(r,σₘ,τ,D₀), 0, atol=1e-10)
            @test isapprox(SD.get_KI_from_s(s, D₀, m), 0, atol=1e-10)
        end
    end

    # TEST KI WITH PRINCIPAL STRESSES
    let m = Model(;setup=loading2D, cm=cmp)
        a = m.cm.damage.r.a
        μ = m.cm.damage.r.μ
        D₀ = m.cm.damage.r.D₀
        pc = m.setup.pc
        ψ = m.cm.damage.r.ψ
        S = (sqrt(1+μ^2) + μ) / (sqrt(1+μ^2) - μ)

        # TEST KI INDEPENDENCE ON SIGN CONVENTION

        # σ₁ & σ₃ negative
        let σ₁ = -S*pc
            σ₃ = -pc
            σs = SD.flaw_shear_stress(σ₁,σ₃,ψ)
            σn = SD.flaw_normal_stress(σ₁,σ₃,ψ)
            @info "friction coef cancelling KI ($(abs(σs/σn))) is not equal to internal friction ($μ)"

            # at D = D0 and σ₁=Sσ₃ -> KI = O
            @test SD.get_KI_from_s(SD.σₐₓ2sₐₓ(m,σ₁), D₀, m) == 0
            # at D > D0 and σ₁=Sσ₃ -> KI = - sign(σ3)*c3*σ₃
            @test SD.get_KI_from_s(SD.σₐₓ2sₐₓ(m,σ₁), D₀ + 0.5(1-D₀), m)  ≈  -sign(σ₁) * sqrt(π*a) * σ₃ * SD.compute_c3(r,D₀ + 0.5(1-D₀))
            @test SD.get_KI_from_s(SD.σₐₓ2sₐₓ(m,σ₁), D₀ + 0.8(1-D₀), m)  ≈  -sign(σ₁) * sqrt(π*a) * σ₃ * SD.compute_c3(r,D₀ + 0.8(1-D₀))
        end
    end

end

# m = Model(;setup=loading2D, cm=cmp)
#         a = m.cm.damage.r.a
#         μ = m.cm.damage.r.μ
#         D₀ = m.cm.damage.r.D₀
#         pc = m.setup.pc
#         ψ = m.cm.damage.r.ψ
# S_func(μ) = (sqrt(1+μ^2) + μ) / (sqrt(1+μ^2) - μ)
#         σ₁ = -S*pc
#         σ₃ = -pc
#         σs = SD.flaw_shear_stress(σ₁,σ₃,ψ)
#         σn = SD.flaw_normal_stress(σ₁,σ₃,ψ)

# μ_vec = 0.3:0.01:0.8
# σ1_vec = [-S_func(μi)*pc for μi in μ_vec]

# σs_vec = [SD.flaw_shear_stress(σ₁i,σ₃,ψ) for σ₁i in σ1_vec]
# σn_vec = [SD.flaw_normal_stress(σ₁i,σ₃,ψ) for σ₁i in σ1_vec]

# lw = 4
# f,ax,pl = lines(μ_vec,abs.(σs_vec./σn_vec),linewidth=lw)
# lines!(ax,μ_vec[1]..μ_vec[end],x->x,linestyle=:dash,linewidth=lw)
# ax.xlabel = "μ"
# ax.ylabel = "τ/σₙ at KI=0"
# ax.xlabelsize = ax.ylabelsize = 30
# DataInspector();

max_stars_repo_name: Leooop/SimpleDamage

max_issues_repo_name: Leooop/SimpleDamage

max_forks_repo_name: Leooop/SimpleDamage


text:
export representatives, genus_representatives

# Genus representatives for quadratic lattices,
#
# With permission ported from Magma package of Markus Kirschmer:
# http://www.math.rwth-aachen.de/~Markus.Kirschmer/magma/lat.html

################################################################################
#
#  Neighbours methods
#
################################################################################

function _all_row_span(M)
  F = base_ring(M)
  rows = Vector{Vector{elem_type(F)}}(undef, nrows(M))
  for i in 1:nrows(M)
    rows[i] = elem_type(F)[M[i, j] for j in 1:ncols(M)]
  end
  n = nrows(M)
  it = Iterators.product([F for i in 1:n]...)
  res = Vector{Vector{elem_type(F)}}()
  for c in it
    z = Ref(c[1]) .* rows[1]
    for i in 2:n
      z = z .+ (Ref(c[i]) .* rows[i])
    end
    push!(res, z)
  end
  return res
end

@doc Markdown.doc"""
    smallest_neighbour_prime(L::HermLat) -> Bool, NfRelOrdIdl, Vector{NfOrdIdl}

Given a hermitian lattice `L`, return `def, P0, bad` such that:

- `def` is `true` if `L` is definite, else `false`;
- `P0` is a prime ideal in the base ring `O` of `L` which is not bad, such that
  `L` is isotropic at `minimum(P0)` and `P0` has smallest minimum among the primes
  satisfying these properties; if `L` is indefinite, `P0` is set to be the trivial ideal;
- `bad` is a vector of prime ideals `p` in the maximal order of the fixed field
  of `L` such that `L_p` is not modular or `p` is dyadic and is not coprime to
  the discriminant of `O`.
"""
function smallest_neighbour_prime(L)
  S = base_ring(L)
  R = base_ring(S)
  lp = bad_primes(L)
  bad = ideal_type(R)[p for p in lp if !ismodular(L, p)[1] ]
  for (p,_) in factor(discriminant(S))
    if isdyadic(p) && !(p in bad)
      push!(bad, p)
    end
  end

  if !isdefinite(L)
    return false, 1*S, bad
  end

  # TODO: This does not find the prime ideal with smallest norm,
  # but with smallest minimum ...

  m = rank(L)
  p = 1
  P = ideal_type(R)[]
  while length(P) == 0
    p = next_prime(p)
    lp = ideal_type(R)[ p[1] for p in prime_decomposition(R, p)]
    P = setdiff(lp, bad)
    if m == 2
      P = filter(p -> isisotropic(L, p), P)
    end
  end
  Q = prime_decomposition(S, P[1])[1][1]

  if isempty(bad)
    I = 1 * S
  else
    I = prod(bad) * S
  end
  n = absolute_norm(Q)
  if n >= 1000
    PP = prime_ideals_up_to(S, 1000)
    for QQ in PP
      if !iscoprime(QQ, I)
        continue
      end

      if isisotropic(L, QQ)
        return true, QQ, bad
      end
    end
  end
  PP = prime_ideals_up_to(S, n)
  for QQ in PP
    if !iscoprime(QQ, I)
      continue
    end
    if isisotropic(L, QQ)
      return true, QQ, bad
    end
  end
  throw(error("Impossible"))
end

function _neighbour(L, B, xG, x, h, P, CC, split)
  R = order(P)
  K = nf(R)
  n = nrows(B)
  a = involution(L)

  local C::ideal_type(R)

  if CC isa Int
    C = split ? a(P) : P
  else
    C = CC
  end

  I = Int[ i for i in 1:n if !iszero(xG[i])]
  i = I[1]
  M = zero_matrix(K, n - length(I), ncols(B))
  for (k, nk) in enumerate(setdiff(1:n, I))
    for j in 1:ncols(B)
      M[k, j] = B[nk, j]
    end
  end
  CI = fractional_ideal_type(R)[ K(1) * R for j in 1:(n - length(I)) ]
  for j in I
    if j != i
      M = vcat(M, matrix(K, 1, ncols(B), elem_type(K)[B[j, k] - a(h\(divexact(xG[j], xG[i]))) * B[i, k] for k in 1:ncols(B)]))
      push!(CI, K(1) * R)
    end
  end
  M = vcat(M, sub(B, i:i, 1:ncols(B)))
  push!(CI, fractional_ideal(order(P), P))
  M = vcat(M, matrix(K, 1, ncols(B), x))
  push!(CI, inv(C))
  pm = pseudo_matrix(M, CI)
  _M = _sum_modules(L, _module_scale_ideal((split ? P : P * C), pseudo_matrix(L)), pm)
  LL = lattice(ambient_space(L), _M)

  @assert islocally_isometric(L, LL, P)
  return LL
end

function stdcallback(list, L)
  keep = all(LL -> !isisometric(LL,L)[1], list)
  return keep, true
end

function eqcallback(list, L)
  keep = all(LL -> LL != L, list)
  return keep, true
end

function _neighbours(L, P, result, max, callback = eqcallback, use_auto = true)
  ok, scale = ismodular(L, P)
  @req ok "The lattice must be locally modular"
  R = base_ring(L)
  K = nf(R)
  a = involution(L)

  if !isdefinite(L) && use_auto
    use_auto = false
  end

  T = local_basis_matrix(L, P, type = :submodule)
  p = minimum(P)
  lp = prime_decomposition(R, p)
  if length(lp) == 1
    C = P
  else
    if lp[1][1] == P
      C = lp[2][1]
    elseif lp[2][1] == P
      C = lp[1][1]
    else
      error("This should not happen.")
    end
  end
  k, h = ResidueField(R, C)
  hext = extend(h, K)
  local form::dense_matrix_type(K)
  form = gram_matrix(ambient_space(L))
  special = false
  if scale != 0
    if isramified(R, minimum(P))
      special = isodd(scale)
      scale = div(scale + 1, 2)
    end
    form = K(elem_in_nf(uniformizer(minimum(P))))^(-scale) * form
  end
  n = rank(L)
  W = vector_space(k, n)

  if use_auto
    G = automorphism_group_generators(L)
    Tinv = inv(T)
    adjust_gens = eltype(G)[T * g * Tinv for g in G]
    adjust_gens_mod_p = dense_matrix_type(k)[map_entries(hext, g) for g in adjust_gens]
    adjust_gens_mod_p = dense_matrix_type(k)[x for x in adjust_gens_mod_p if !isdiagonal(x)]
    if length(adjust_gens_mod_p) > 0
      _LO = line_orbits(adjust_gens_mod_p)
      LO = Vector{eltype(k)}[x[1] for x in _LO]
    else
      LO = enumerate_lines(k, n)
    end
  else
    LO = enumerate_lines(k, n)
  end

  keep = true
  cont = true

  if P != C
    _G = T * form * _map(transpose(T), a)
    G = map_entries(hext, _G)
    pi = p_uniformizer(P)
    pih = h(pi)
    for w in LO
      x = elem_type(K)[ sum(T[i, j] * (hext\w[i]) for i in 1:n) for j in 1:ncols(T)]
      LL = _neighbour(L, T, pih * matrix(k, 1, length(w), w) * G, K(pi) .* x, hext, P, C, true)
      keep, cont = callback(result, LL)
      if keep
        push!(result, LL)
      end
      if !cont || (length(result) >= max)
        break
      end
    end
  elseif special
    pi = uniformizer(P)
    _G = elem_in_nf(pi) * T * form * _map(transpose(T), a)
    G = map_entries(hext, _G)
    for w::Vector{fq} in LO
      Gw = G * matrix(k, length(w), 1, w)
      ok = 0
      for d in 1:n
        if !iszero(Gw[d])
          ok = d
          break
        end
      end
      @assert ok != 0
      x = elem_type(K)[ sum(T[i, j] * (hext\w[i]) for i in 1:n) for j in 1:ncols(T)]
      nrm = _inner_product(form, x, x, a)
      b = hext\(-hext(nrm) // (2*Gw[ok, 1]))
      for j in 1:ncols(T)
        x[j] = x[j] + b * elem_in_nf(pi) * T[ok, j]
      end
      nrm = _inner_product(form, x, x, a)
      LL = _neighbour(L, T, matrix(k, 1, length(w), w) * G, x, hext, P, P, false)
      keep, cont = callback(result, LL)
      if keep
        push!(result, LL)
      end
      if !cont || length(result) >= max
        break
      end
    end
  else
    _G = T * form * _map(transpose(T), a)
    G = map_entries(hext, _G)
    ram = isramified(R, minimum(P))
    if ram
      pi = uniformizer(P)
      S = [ h\x for x in k ]
    else
      p = minimum(P)
      pi = uniformizer(p)
      kp, hp = ResidueField(order(p), p)
      alpha = h\(degree(k) == 1 ? one(k) : gen(k))
      Tram = matrix(kp, 2, 1, [2, hp(tr(alpha))])
    end
    for w::Vector{fq} in LO
      __w = [ (hext\w[i]) for i in 1:n]
      x = [ sum(T[i, j] * (__w[i]) for i in 1:n if !iszero(w[i])) for j in 1:ncols(T)]
      nrm = _inner_product(form, x, x, a)
      if !(nrm in P)
        continue
      end
      wG = matrix(k, 1, length(w), w) * G
      ok = 0
      for j in 1:n
        if !iszero(wG[j])
          ok = j
          break
        end
      end
      @assert ok != 0
      if !ram
        el = order(p)(base_field(K)(nrm)//pi)
        b, s, V = can_solve_with_kernel(Tram, matrix(kp, 1, 1, [hp(-el)]), side = :left)
        @assert b
        @assert s * Tram == matrix(kp, 1, 1, [hp(-el)])
        _kernel = [ matrix(kp, 1, 2, v) for v in _all_row_span(V)]
        l = a(hext\(inv(wG[ok])))
        S = elem_type(K)[ l * (hext\((s + v)[1]) + (hext\(s + v)[2])*alpha) for v in _kernel ]
      end
      for s in S
        LL = _neighbour(L, T, wG, elem_type(K)[x[o] + K(elem_in_nf(pi))*s*T[ok, o] for o in 1:ncols(T)], hext, P, P, false)
        keep, cont = callback(result, LL)
        if keep
          push!(result, LL)
        end
        if !cont || (length(result) >= max)
          break
        end
      end
    end
  end
  return result
end

@doc Markdown.doc"""
    neighbours(L::HermLat, P::NfRelOrdIdl, max = inf) -> Vector{HermLat}

Return the immediate `P`-neighbours of `L`. At most `max` neighbours are returned.

If `L` is definite, this function uses by default the automorphism group of `L`. If
`L` is indefinite, the use of the automorphism group is automatically disabled.
"""
function neighbours(L::HermLat, P, max = inf)
  @req order(P) == base_ring(L) "Arguments are incompatible"
  @req isprime(P) "Second argument must be prime"
  @req !isramified(order(P), minimum(P)) || !Hecke.isdyadic(minimum(P)) "Second argument cannot be a ramified prime over 2"
  @req ismodular(L, P)[1] "The lattice must be locally modular"
  @req rank(L) >= 2 "The rank of the lattice must be at least 2"
  @req Hecke.isisotropic(L, P) "The lattice must be locally isotropic"

  return _neighbours(L, P, [], max)
end

@doc Markdown.doc"""
    iterated_neighbours(L:HermLat, P::NfRelOrdIdl; use_auto = false, max = inf,
				                   callback = eqcallback,
						   missing_mass = Ref{fmpq}(zero(fmpq)))
                                                                            -> Vector{HermLat}

Return a set of representatives of `N(L,P)` (see [Kir16, Definition 5.2.6]). At most
`max` representatives are returned.

The use of the automorphism group of `L` is disabled by default. If `use_auto` is set on
`true`, the function uses the automorphism group in the definite case; in the indefinite
case, this keyword has no effect.
If `callback = false`, it uses `stdcallback` in the case where `L` is definite, `eqcallback`
otherwise. By defaut, the use of the mass is disabled.
"""
function iterated_neighbours(L::HermLat, P; use_auto = false, max = inf,
                                            callback = false,
                                            missing_mass = Ref{fmpq}(zero(fmpq)))
  @req order(P) == base_ring(L) "Arguments are incompatible"
  @req isprime(P) "Second argument must be prime"
  @req !isramified(order(P), minimum(P)) || !Hecke.isdyadic(minimum(P)) "Second argument cannot be a ramified prime over 2"
  @req ismodular(L, P)[1] "The lattice must be locally modular"
  @req rank(L) >= 2 "The rank of the lattice must be at least 2"
  @req Hecke.isisotropic(L, P) "The lattice must be locally isotropic"

  if callback == false && isdefinite(L)
    _callback = stdcallback
  elseif callback == false && !isdefinite(L)
    _callback = eqcallback
  else
    _callback = callback
  end

  result = typeof(L)[ L ]

  use_mass = !iszero(missing_mass[])

  if use_mass
    _mass = missing_mass[] - 1//automorphism_group_order(L)
  end

  i = 1
  oldlength = length(result)
  while length(result) < max && i <= length(result)
    result = _neighbours(result[i], P, result, max, _callback, use_auto)
    no_lattices = length(result) - oldlength
    oldlength = length(result)
    if use_mass && no_lattices > 0
      _mass = _mass - sum(fmpq[1//automorphism_group_order(result[i]) for i in (length(result) - no_lattices + 1):length(result)])
      if iszero(_mass)
        break
      end
    end
    if use_mass && _mass < 0
      throw(error("This should not happen"))
    end
    i = i + 1
  end
  if use_mass
    missing_mass[] = _mass
  end
  return result
end

@doc Markdown.doc"""
    neighbours_with_ppower(L::HermLat, P::NfRelOrdIdl, e::Integer, use_auto = true)
                                                                      -> Vector{HermLat}

Return a sequence of `P`-neighbours of length `e`, `L=L_1, L_2, \dots, L_e` such that
`L_{i-1} != L_{i+1}` for `i = 2, \dots, e-1` (see [Kir19, Algorithm 4.7.]).

If the lattice is definite, the use of the automorphism group is by default enabled.
In the indefinite case, the automorphism group is not used.
"""
function neighbours_with_ppower(L, P, e)
  result = typeof(L)[]
  for i = 1:e
    if i == 1
      L = neighbours(L, P, 1)[1]
    else
      N = neighbours(L, P,  2)
      L = N[1] == result[end] ? N[2] : N[1]
    end
  push!(result, L)
  end
  return result
end

################################################################################
#
#  Genus representatives for hermitian lattices
#
################################################################################

@doc Markdown.doc"""
    genus_generators(L::HermLat) -> Vector{Tuple{NfRelOrdIdl, fmpz}}, Bool,
                                    NfRelOrdIdl

Given a hermitian lattice `L`, return `gens, def, P0` such that:

- `gens` is a vector of tuples `(P,e)` consisting of a prime ideal `P` in the base ring of `L`
  and an integer `e \geq 2` which can be used to compute the ideal `\mathfrak A` in line 11
  of [Kir19, Algorithm 4.7.]);
- `def` is `true` if `L` is definite, else `false`;
- `P0` is a prime ideal in the base ring of `L` which is not bad, such that
  `L` is isotropic at `minimum(P0)` and `P0` has smallest minimum among the primes
  satisfying these properties.
"""
function genus_generators(L::HermLat)
  R = base_ring(L)
  E = nf(R)
  D = different(R)
  a = involution(L)
  def, P0, bad = smallest_neighbour_prime(L)

  local bad_prod::ideal_type(base_ring(R))

  if isempty(bad)
    bad_prod = 1 * base_ring(R)
  else
    bad_prod = prod(bad)
  end

  # First the ideals coming from the C/C0 quotient
  Eabs, EabstoE = absolute_simple_field(ambient_space(L))
  Rabs = maximal_order(Eabs)
  C, h = class_group(Rabs)
  RR = base_ring(R)
  C0 = support(D)
  CC, hh = class_group(RR)
  for p in find_gens(pseudo_inv(hh), PrimesSet(2, -1))[1]
    if !(p in C0)
      push!(C0, sum(R * R(EabstoE(elem_in_nf(b))) for b in basis(p)))
    end
  end
  Q0, q0 = quo(C, elem_type(C)[ h\ideal(Rabs, [Rabs(EabstoE\b) for b in absolute_basis(i)]) for i in C0])
  q00 = pseudo_inv(q0) * h
  PP = ideal_type(R)[]

  local F::GaloisField

  local W::Generic.QuotientModule{gfp_elem}

  if iseven(rank(L))
    for (P, e) in factor(D)
      p = minimum(P)
      G = genus(L, p)
      if any(i -> isodd(rank(G, i)), 1:length(G))
        continue
      elseif !isdyadic(p)
        if any(i -> iseven(scale(G, i)), 1:length(G))
          continue
        end
      else
        if any(i -> isodd(e  + scale(G, i)) || (e + scale(G, i) != G.ni[i]), 1:length(G))
          continue
        end
      end
      push!(PP, P)
    end

    if !isempty(PP)
      U, f = unit_group_fac_elem(Rabs)
      UU, ff = unit_group_fac_elem(RR)
      nnorm = hom(U, UU, GrpAbFinGenElem[ff\FacElem(nf(RR)(norm(f(U[i])))) for i in 1:ngens(U)])
      l = length(PP)
      VD = Int[ valuation(D, P) for P in PP ]
      K, k = kernel(nnorm)
      F = GF(2, cached = false)
      V = VectorSpace(F, length(PP))
      S = elem_type(V)[]
      for u in gens(K)
        z = elem_type(F)[]
        for i in 1:length(PP)
          zz = R(EabstoE(evaluate(f(k(u))))) - 1
          if iszero(zz) || (valuation(zz, PP[i]) >= VD[i])
            push!(z, F(0))
          else
            push!(z, F(1))
          end
        end
        push!(S, V(z)::elem_type(V))
      end
      _T, _ = sub(V, S)
      W, w = quo(V, _T)
      if dim(W) == 0
        PP = ideal_type(R)[]
      end
    end
  end

  Gens = Tuple{ideal_type(R), fmpz}[]

  if isempty(PP)
    S = GrpAbFinGenElem[]
    Q, q = quo(Q0, S)::Tuple{GrpAbFinGen, GrpAbFinGenMap}
    Work = def ? typeof(P0)[ P0 ] : typeof(P0)[]
    p = 2
    while order(Q) > 1
      while isempty(Work)
        p = next_prime(p)
        Work = ideal_type(R)[ QQ for QQ in support(p * R) if length(prime_decomposition(R,minimum(QQ))) == 2 && valuation(bad_prod, minimum(QQ)) == 0 ]
      end
      P = popfirst!(Work)
      c = (q00\(EabstoE\P))::GrpAbFinGenElem
      o = order(q(c))::fmpz
      if !isone(o)
        push!(S, c)
        Q, q = quo(Q0, S)::Tuple{GrpAbFinGen, GrpAbFinGenMap}
        push!(Gens, (P, o))
      end
    end
  else
    ll = Int(order(Q0))
    cocycle = Matrix{elem_type(W)}(undef, ll, ll)
    for i in 1:ll
      for j in 1:ll
        cocycle[i,j] = zero(W)
      end
    end
    C = collect(Q0)
    ideals = ideal_type(Rabs)[ q00(C[i]) for i in 1:length(C) ]
    for i in 1:ll
      for j in 1:i
        ij = findfirst(isequal(C[i] + C[j]), C)
        Iabs = ideals[i] * ideals[j] * inv(ideals[ij])
        I = EabstoE(Iabs)
        J = I * inv(a(I))
        Jabs = EabstoE\J
        ok, x = isprincipal(Jabs)
        u = f(nnorm\(-(ff\FacElem(nf(RR)(norm(x))))))
        x = x * u
        @assert norm(x) == 1
        if evaluate(x) == 1
          y = w(V(zeros(F,length(PP))))
        else
          y = w(V([ valuation(EabstoE(evaluate(x) - 1), PP[i]) >= VD[i] ? F(0) : F(1) for i in 1:length(PP)]))
        end
        cocycle[i, j] = y
        cocycle[j, i] = y
      end
    end

    S = Tuple{elem_type(Q0), Generic.QuotientModuleElem{elem_type(F)}}[(id(Q0), zero(W))]
    Work = def ? typeof(P0)[ P0 ] : typeof(P0)[]
    p = 2
    sizW = dim(W) * characteristic(W.base_ring)
    while length(S) != order(Q0) * sizW
      while isempty(Work)
        p = next_prime(p)
        Work = ideal_type(R)[ QQ for QQ in support(p * R) if length(prime_decomposition(R, minimum(QQ))) == 2 && valuation(bad_prod, minimum(QQ)) == 0 ]
      end
      P = popfirst!(Work)
      c = q00\(EabstoE\P)
      i = findfirst(isequal(c), C)
      Iabs = EabstoE\P * inv(ideals[i])
      I = EabstoE(Iabs)
      J = I * inv(a(I))
      Jabs = EabstoE\J
      ok, x = isprincipal(Jabs)
      u = f(nnorm\(-(ff\FacElem(nf(RR)(norm(x))))))
      x = x * u
      @assert norm(x) == 1
      if evaluate(x) == 1
        y = zeros(F,length(PP))
      else
        y = [ valuation(EabstoE(evaluate(x) - 1), PP[i]) >= VD[i] ? F(0) : F(1) for i in 1:length(PP)]
      end
      idx = findfirst(isequal(P), PP)
      if idx !== nothing
        y = elem_type(F)[i == idx ? y[i] + 1 : y[i]  for i in 1:dim(V)]
      end
      elt = (c, w(V(y)))
      elt1 = elt
      o = 1
      siz = length(S)
      while !(elt1 in S)
        j = findfirst(isequal(elt1[1]), C)
        @assert !isnothing(j)
        for l in 1:siz
          elt2 = S[l]
          k = findfirst(isequal(elt2[1]), C)
          @assert !isnothing(k)
          prod = (elt1[1] + elt2[1], elt1[2] + elt2[2] + cocycle[j, k])
          if !(prod in S)
            push!(S, prod)
          end
        end
        elt1 = (elt[1] + elt1[1], elt[2] + elt1[2] + cocycle[i, j])
        o = o + 1
      end
      @assert length(S) == siz * o
      if o != 1
        push!(Gens, (P, o))
      end
    end
  end
  return Gens, def, P0
end

@doc Markdown.doc"""
    genus_representatives(L::HermLat; max = inf, use_auto = true,
                                                 use_mass = false)
                                                          -> Vector{HermLat}

Compute representatives for the isometry classes in the genus of $L$. At most
`max` representatives are returned.

If `L` is definite, the use of the automorphism group of `L` is enabled by default.
It can be disabled by `use_auto = false`. In the case where `L` is indefinite, the entry
`use_auto` has no effect. The computation of the mass can be enabled by `use_mass = true`.
"""
function genus_representatives(L::HermLat; max = inf, use_auto::Bool = true,
                                                      use_mass::Bool = false)
  @req rank(L) >= 2 "Lattice must have rank >= 2"
  R = base_ring(L)
  gens, def, P0 = genus_generators(L)
  a = involution(L)
  LL = typeof(L)[ L ]
  for g in gens
    if def && g[1] == P0
      continue
    end
    I = g[1]^Int(g[2] - 1)
    J = inv(a(I))
    N = neighbours_with_ppower(L, g[1], g[2] - 1)
    inter = typeof(L)[]
    for i in 2:length(LL)
      M = pseudo_matrix(LL[i])
      IM = _module_scale_ideal(I,M)
      JM = _module_scale_ideal(J,M)
      inter = append!(inter, lattice(ambient_space(L), _intersect_modules(_sum_modules(IM, pseudo_matrix(x)), JM) for x in N))
    end
    LL = vcat(LL, N, inter)
  end
  @assert length(LL) == prod(Int[g[2] for g in gens if !def || g[1] != P0])
  @assert all(X -> genus(X) == genus(L), LL)

  if use_mass
    mass = Hecke.mass(L)
  else
    mass = zero(fmpq)
  end

  missing_mass = Ref(mass)

  local result::Vector{typeof(L)}

  if def
    result = typeof(L)[]
    for L in LL
      neig = iterated_neighbours(L, P0, max = max, use_auto = use_auto, missing_mass = missing_mass)
      append!(result, neig)
      max = max - length(result)
    end
  else
    result = LL
  end
  return result
end

@doc Markdown.doc"""
    representatives(G::GenusHerm) -> Vector{HermLat}

Return representatives for the isometry classes in $G$.
"""
function representatives(G::GenusHerm)
  return genus_representatives(representative(G))
end


text:
using Plots, Test

@testset "Shape Copy" begin
    square = Shape([(0,0),(1,0),(1,1),(0,1)])
    square2 = Shape(square)
    @test square2.x == square.x
    @test square2.y == square.y
end

@testset "Shape Center" begin
    square = Shape([(0,0),(1,0),(1,1),(0,1)])
    @test Plots.center(square) == (0.5,0.5)
end

@testset "Shape Translate" begin
    square = Shape([(0,0),(1,0),(1,1),(0,1)])
    squareUp = Shape([(0,1),(1,1),(1,2),(0,2)])
    squareUpRight = Shape([(1,1),(2,1),(2,2),(1,2)])

    @test Plots.translate(square,0,1).x == squareUp.x
    @test Plots.translate(square,0,1).y == squareUp.y

    @test Plots.center(translate!(square,1)) == (1.5,1.5)
end

@testset "Brush" begin
    @testset "Colors" begin
        baseline = brush(1, RGB(0,0,0))
        @test brush(:black) == baseline
        @test brush("black") == baseline
    end
    @testset "Weight" begin
        @test brush(10).size == 10
        @test brush(0.1).size == 1
    end
    @testset "Alpha" begin
        @test brush(0.4).alpha == 0.4
        @test brush(20).alpha == nothing
    end
    @testset "Bad Argument" begin
        # using test_logs because test_warn seems to not work anymore
        @test_logs (:warn,"Unused brush arg: nothing (Nothing)") begin
            brush(nothing)
        end
    end
end


text:
abstract type Manifold end
struct TangentBundle{M<:Manifold} <: Manifold end
const T = TangentBundle
struct ProductManifold{M1<:Manifold,M2<:Manifold} <: Manifold end
const PM = ProductManifold

abstract type ChartID end
# Usually for regular manifolds, I <: ChartID or type(I) = Int
# For product manifolds, I <: Tuple{Int,Int}
struct Chart{I, M<:Manifold} end
ChartList = Vector{Chart}

abstract type TaskMapX{A,B,S} end
abstract type TaskMap{M<:Manifold, N<:Manifold, S} <: TaskMapX{M,N,S} end
abstract type TaskMapT{M<:Manifold, N<:Manifold, S} <: TaskMapX{M,N,S} end
Base.eltype(::TaskMapX{A,B,S}) where {A,B,S} = S

abstract type TaskX{F} end
abstract type Task{F<:TaskMapX} <: TaskX{F} end
abstract type TaskGDS{F<:TaskMapT} <: TaskX{F} end
Base.eltype(::Task{F}) where F = eltype(F)

struct ProductTaskMap{T1<:TaskMap,T2<:TaskMap,S} <: TaskMapX{T1,T2,S}
    map1::T1
    map2::T2
end
const PTM = ProductTaskMap
PTM(map1::T1,map2::T2) where {T1,T2} = PTM{T1,T2,eltype(T1)}(map1,map2)
Base.eltype(::Type{PTM{T1,T2,S}}) where {T1,T2,S} = S
BaseTaskMap = Union{TaskMap, ProductTaskMap}

struct ProductTaskMapT{T1<:TaskMap,T2<:TaskMap,S} <: TaskMapX{T1,T2,S}
    base_map::ProductTaskMap{T1,T2,S}
end
ProductTaskMapT(T1,T2) = ProductTaskMapT(ProductTaskMap(T1,T2))
const PTMT = ProductTaskMapT

ProductTaskMapX = Union{PTM, PTMT}
PTMX = ProductTaskMapX

TaskList = Vector{Task{<:BaseTaskMap}}
TaskTList = Vector{Task{<:TaskMapT}}
TaskGDSList = Vector{TaskGDS{<:TaskMapT}}
TaskXList = Union{TaskList, TaskTList, TaskGDSList}

# Traits for specifying base definitions of task maps
# Also for specifying coordinate representation used for manifold chart selection,
# task trajectory logging
abstract type CoordinateRep end
struct EmbRep <: CoordinateRep end
struct ChartRep <: CoordinateRep end


text:
function walkup(f::Function, dir0::AbstractString;
        root = homedir(),
        onerr::Function = _default_onerr,
    )

    root = abspath(root)
    dir0 = abspath(dir0)

    content = try; readdir(dir0)
        catch err; (onerr(dir0, err) === true) && return
    end

    # walk content
    path::String = ""
    subi = firstindex(content)
    sub1 = lastindex(content)
    while subi <= sub1
        try
            for _ in subi:sub1
                name = content[subi]
                subi += 1

                path = joinpath(dir0, name)
                (f(path) === true) && return
            end
        catch err;
            (onerr(path, err) === true) && return
        end
    end

    # base
    (dir0 == root) && return nothing

    # recursive
    walkup(f, dirname(dir0); root, onerr)
end


text:
module ForwardCurveSmoother

using DataFrames
using CSV
using Dates
using JuMP
using OSQP
using GLM
using GLPK
using TimeSeriesInterface
using Lasso

include("structs.jl")
include("optimization.jl")
include("utils.jl")

function fit_structural_model(contracts::DataFrame, spot::DataFrame, structural_parameters::StructuralParameters)

    # Creating input and running the smoothing algorithm
    input_structural, coeficients = create_structural_model_input(contracts, spot,
                                                                  structural_parameters)
    output_structural = forward_structural_model(input_structural)

    # Handle results
    reconstructed_prices = handle_data_frame_of_reconstructed_prices(contracts, output_structural.reconstructed_forward_prices)
    elementary_prices    = output_structural.structural_dinamic .+ output_structural.structural_error
    elementary_errors    = output_structural.structural_error

    number_of_maturities = size(elementary_prices, 2)
    timestamps           = contracts[:, 1]
    names_columns        = [:timestamps; [Symbol("maturity_$i") for i in 1:number_of_maturities]]

    df_elementary_prices = DataFrame(hcat(timestamps, elementary_prices), names_columns)
    df_elementary_errors = DataFrame(hcat(timestamps, elementary_errors), names_columns)

    # Saving CSVs on the current path
    @info("Saving reconstructed_prices.csv...")
    CSV.write("Results/reconstructed_prices.csv", reconstructed_prices; delim = ';')
    @info("Saving elementary_prices.csv...")
    CSV.write("Results/elementary_prices.csv", df_elementary_prices; delim = ';')
    @info("Saving elementary_errors.csv...")
    CSV.write("Results/elementary_errors.csv", df_elementary_errors; delim = ';')

    return input_structural, output_structural
end

end

max_stars_repo_name: dietzemarina/ForwardCurveSmoother

max_issues_repo_name: dietzemarina/ForwardCurveSmoother

max_forks_repo_name: dietzemarina/ForwardCurveSmoother


text:
 using JuPOT

             # Collect data from ARGS array
             data = float(split(ARGS[1], ","))
             returns = round(float(split(ARGS[2], ",")),3)
             tickers = split(ARGS[3], ",")
             target_return = float(ARGS[4])

             # Calculate number of stocks and days from data
             stocks = length(tickers)
             days = Int(length(data)/stocks)

             # Form data into matrix
             matrix = reshape(data,days,stocks)

             # Produce covariance matrix
             covariance = cov(matrix)

             # Assets data structure containing, tickers, expected returns and covarariance
             assets = AssetsCollection(tickers, returns, covariance)

             # Create simple MVO object using assets and target_return
             mvo = SimpleMVO(assets, target_return; short_sale=false)

             # Optimise MVO object and print result to file
             outfile = "output.txt"
                          f = open(outfile, "w")
                          println(f, optimize(mvo))
                          close(f)

             # No. Of Assets
             # n = 3
             # No. of days data
             # m = 6
             # d = [0.5,0.7,0.3,0.2,0.56,0.26,0.2,0.4,0.6,0.2,0.15,0.58,0.8,0.3,0.9,0.2,0.51,0.15]
             # Returns a matrix of size(n) with entries between 0-1s
             # returns =  [0.5,0.4,0.6]
             # matrix = reshape(d,m,n)
             #covariance matrix of returns
             # covariance = cov(matrix)
             # List of asset names
             # tickers = [randstring(3) for i in 1:n]

             # Assets data structure containing, names, expected returns, covarariance
             # assets = AssetsCollection(tickers, returns, covariance)

             # target_return = 0.2
             # mvo = SimpleMVO(assets, target_return; short_sale=false)

             # println(optimize(mvo))


text:
mutable struct CUDA_Resize
    input::CuArray{Float32,4}

    xpass::CuArray{Float32,4}
    ypass::CuArray{Float32,4}
    gauss_x::CuArray{Float32,1}
    gauss_y::CuArray{Float32,1}

    resized::CuArray{Float32,4}
    flip_xy::Bool
    output::CuArray{Float32,4}
end

function CUDA_Resize(in_w,in_h,out_w,out_h,n,flip_xy=true)
    input = convert(CuArray,zeros(Float32,in_w,in_h,1,n))

    xpass = convert(CuArray,zeros(Float32,in_w,in_h,1,n))
    ypass = convert(CuArray,zeros(Float32,in_w,in_h,1,n))
    (gauss_x,gauss_y) = get_xy_gauss_cu((in_w,in_h),(out_w,out_h))

    resized = convert(CuArray,zeros(Float32,out_w,out_h,1,n))
    output = convert(CuArray,zeros(Float32,out_w,out_h,1,n))

    CUDA_Resize(input,xpass,ypass,gauss_x,gauss_y,resized,flip_xy,output)
end

function lowpass_resize(cr::CUDA_Resize,input::AbstractArray{T,4}) where T

    cr.input[:] = convert(CuArray,input)

    CUDA_blur_x(cr.input,cr.xpass,cr.gauss_x)
    CUDA_blur_y(cr.xpass,cr.ypass,cr.gauss_y)
    CUDA_resize4(cr.ypass,cr.resized)
    if cr.flip_xy
        CUDA_flip_xy(cr.resized,cr.output)
    else
        cr.output[:] = cr.resized
    end

    nothing
end


text:
export SearchClient, create_index, get_search_client, drop_index

"""
    Client(index_name::AbstractString; kwargs....) -> Jedis.Client

Set the RediSearch Client. Clients contain an index name and a Jedis Client Struct.
Kwargs are included as any Jedis client kwargs

# Fields
- `index_name::String`: Index name that will be created with the client.
- `client::Jedis.Client`: Client object made through the Jedis package.

# Examples
```julia-repl
julia> client = Client("myIdx");

julia> client = Client("myIdx"; host="localhost", port=6379);
```
"""
mutable struct SearchClient
    index_name::String
    client::Jedis.Client
end

function SearchClient(index_name::AbstractString; kwargs...)::SearchClient
    client = SearchClient(index_name, Jedis.Client(; kwargs...))
    SEARCHCLIENT[] = client
    return client
end

"""
    Global client of RediSearch. Used to reference the base client for RediSearch made
    through `Client()`.
"""
const SEARCHCLIENT = Ref{SearchClient}()

""" Return the client object for RediSearch"""
get_search_client() = SEARCHCLIENT[]

"""
    create_index(
        fields;
        client=get_search_client(),
        no_term_offsets=false,
        no_field_flags=false,
        stopwords=nothing,
        definition=nothing,
        max_text_fields=false,
        temporary=nothing,
        no_highlight=false,
        no_term_frequencies=false,
        skip_initial_scan=false,
    )
    create_index(field::AbstractField; kwargs...) = create_index([field]; kwargs...)

Creates an index with the given spec. If an index already exists, this will error.

# Arguments
- `fields`: List of Field objects

# Keywords
- `client`: Search client
- `no_term_offsets`: Bool for offsets.
- `no_field_flags`: Bool for field flags
- `stopwords`: List of stopwords. If nothing, default stopwords are used.
- `definition`: Index Definition object.
- `max_text_fields`: Forces RediSearch to encode indexes as if there were more than 32
text attributes.
- `temporary`: Create a lightweight temporary index which will expire after the specified
period of inactivity.
- `no_highlight`: Conserves storage space and memory by disabling highlighting support.
- `no_term_frequencies`: If true, we avoid saving the term frequencies in the index.
- `skip_initial_scan`: If true, we do not scan and index.
"""
function create_index(
    fields;
    client=get_search_client(),
    no_term_offsets=false,
    no_field_flags=false,
    stopwords=nothing,
    definition=nothing,
    max_text_fields=false,
    temporary=nothing,
    no_highlight=false,
    no_term_frequencies=false,
    skip_initial_scan=false,
)
    args = []
    append!(args, [SearchCommands.CREATE_CMD, client.index_name])
    redis_client = client.client

    if !isnothing(definition)
        append!(args, definition.args)
    end

    if max_text_fields
        push!(args, SearchCommands.MAXTEXTFIELDS)
    end

    if !isnothing(temporary) && temporary isa Integer
        append!(args, [SearchCommands.TEMPORARY, temporary])
    end

    if no_term_offsets
        push!(args, SearchCommands.NOOFFSETS)
    end

    if no_highlight
        push!(args, SearchCommands.NOHL)
    end

    if no_field_flags
        push!(args, SearchCommands.NOFIELDS)
    end

    if no_term_frequencies
        push!(args, SearchCommands.NOFREQS)
    end

    if skip_initial_scan
        push!(args, SearchCommands.SKIPINITIALSCAN)
    end

    if !isnothing(stopwords) && isa_list(stopwords)
        num_stopwords = length(stopwords)
        append!(args, [SearchCommands.STOPWORDS, num_stopwords])

        if num_stopwords > 0
            append!(args, [word for word in stopwords])
        end
    end

    push!(args, "SCHEMA")
    schema_args = []

    for field in fields
        append!(schema_args, redis_args(field))
    end

    append!(args, schema_args)
    return Jedis.execute(args, redis_client)
end

"""
    drop_index(index; client=get_search_client())

Drop an index against a Client.

# Arguments
- `index`: Index name to be dropped

# Keywords
- `client::SearchClient`: RediSearch Client
"""
function drop_index(index; client::SearchClient=get_search_client())
    return Jedis.execute([SearchCommands.DROP_CMD, index], client.client)
end
drop_index() = drop_index(get_search_client().index_name)