text:
##### Beginning of file

module GitHubSrcHost # Begin submodule Snapshots.Hosts.GitHubSrcHost

__precompile__(true)

import ..Types
import ..Utils

import Dates
import GitHub
import TimeZones

import ..delayederror

function new_github_session(
        ;
        github_organization::AbstractString,
        github_bot_username::AbstractString,
        github_bot_personal_access_token::AbstractString,
        )::Function

    _github_organization::String = strip(
        convert(String, github_organization)
        )
    _provided_github_bot_username::String = strip(
        convert(String, github_bot_username,)
        )
    _github_bot_personal_access_token::String = strip(
        convert(String, github_bot_personal_access_token)
        )
    function _get_github_username(auth::GitHub.Authorization)::String
        user_information::AbstractDict = GitHub.gh_get_json(
            GitHub.DEFAULT_API,
            "/user";
            auth = auth,
            )
        username::String = user_information["name"]
        username_stripped::String = strip(username)
        return username_stripped
    end

    @info("Attempting to authenticate to GitHub...")
    auth::GitHub.Authorization = GitHub.authenticate(
        _github_bot_personal_access_token
        )
    _github_username::String = _get_github_username(auth)
    if lowercase(strip(_github_username)) !=
            lowercase(strip(_provided_github_bot_username))
        delayederror(
            string(
                "Provided GitHub username ",
                "(\"$(_provided_github_bot_username)\") ",
                "does not match ",
                "actual GitHub username ",
                "(\"$(_github_username)\").",
                )
            )
    else
        @info(
            string(
                "Provided GitHub username matches ",
                "actual GitHub username.",
                ),
            _provided_github_bot_username,
            _github_username,
            )
    end
    @info("Successfully authenticated to GitHub :)")

    @info(
        string(
            "GitHub username: ",
            "$(_get_github_username(auth))",
            )
        )
    @info(
        string(
            "GitHub organization: ",
            "$(_github_organization)",
            )
        )

    repository_owner = GitHub.owner(
        _github_organization,
        true;
        auth = auth,
        )

    function _list_all_repos()::Vector{String}
        @info("Loading the list of all of repos in the organization")
        full_repo_list::Vector{GitHub.Repo} = GitHub.Repo[]
        need_to_continue::Bool = true
        current_page_number::Int = 1
        while need_to_continue
            repos, page_data = GitHub.repos(
                _github_organization,
                true;
                auth = auth,
                page_limit = 1,
                params = Dict(
                    "per_page" => 100,
                    "page" => current_page_number,
                    ),
                )
            if length(repos) == 0
                need_to_continue = false
            else
                append!(
                    full_repo_list,
                    repos,
                    )
                need_to_continue = true
                current_page_number += 1
            end
        end
        unique_repo_list::Vector{GitHub.Repo} = unique(full_repo_list)
        unique_name_list = [x.name for x in unique_repo_list]
        unique!(unique_name_list)
        sort!(unique_name_list)
        return unique_name_list
    end

    function _repo_name_with_org(
            ;
            repo::AbstractString,
            org::AbstractString,
            )::String
        repo_name_without_org::String = _repo_name_without_org(
            ;
            repo = repo,
            org = org,
            )
        org_stripped::String = strip(
            strip(strip(strip(strip(convert(String, org)), '/')), '/')
            )
        result::String = string(
            org_stripped,
            "/",
            repo_name_without_org,
            )
        return result
    end

    function _repo_name_without_org(
            ;
            repo::AbstractString,
            org::AbstractString,
            )::String
        repo_stripped::String = strip(
            strip(strip(strip(strip(convert(String, repo)), '/')), '/')
            )
        org_stripped::String = strip(
            strip(strip(strip(strip(convert(String, org)), '/')), '/')
            )
        if length(org_stripped) == 0
            result = repo_stripped
        else
            repo_stripped_lowercase::String = lowercase(repo_stripped)
            org_stripped_lowercase::String = lowercase(org_stripped)
            org_stripped_lowercase_withtrailingslash::String = string(
                org_stripped_lowercase,
                "/",
                )
            if startswith(repo_stripped_lowercase,
                    org_stripped_lowercase_withtrailingslash)
                index_start =
                    length(org_stripped_lowercase_withtrailingslash) + 1
                result = repo_stripped[index_start:end]
            else
                result = repo_stripped
            end
        end
        return result
    end

    function _get_src_url(
            ;
            repo_name::String,
            credentials::Symbol,
            )::String
        repo_name_without_org::String = _repo_name_without_org(
            ;
            repo = repo_name,
            org = _github_organization,
            )
        result::String = ""
        if credentials == :with_auth
            result = string(
                "https://",
                _github_username,
                ":",
                _github_bot_personal_access_token,
                "@",
                "github.com/",
                _github_organization,
                "/",
                repo_name_without_org,
                )
        elseif credentials == :with_redacted_auth
            result = string(
                "https://",
                _github_username,
                ":",
                "*****",
                "@",
                "github.com/",
                _github_organization,
                "/",
                repo_name_without_org,
                )
        elseif credentials == :without_auth
            result =string(
                "https://",
                "github.com/",
                _github_organization,
                "/",
                repo_name_without_org,
                )
        else
            delayederror("$(credentials) is not a supported value for credentials")
        end
        return result
    end

    function _github_provider(task::Symbol)::Function
        if task == :list_all_repos
            return _list_all_repos
        elseif task == :get_src_url
            return _get_src_url
        else
            delayederror("$(task) is not a valid task")
        end
    end

    return _github_provider
end

end # End submodule Snapshots.Hosts.GitHubSrcHost

##### End of file


text:
struct IntegrationPath{T,P}
  center::P #
  inner_radius::T
  outer_radius::T
  segments::Vector{Tuple{P,P}}
  arcs::Vector{Tuple{P,P,Int}}
  circles::Vector{Int}
  normal::P
  height::T #
  projected_center::P
  plane_inner_radius::T
  plane_outer_radius::T
end

const CLOCKWISE = -1
const COUNTERCLOCKWISE = +1

function contour(p1, p2, p3, center)

    Vertex = typeof(p1)
    segments = Tuple{Vertex,Vertex}[
        (p1,p2),
        (p2,p3),
        (p3,p1),
    ]
    arcs = Tuple{Vertex,Vertex,Int}[]
    circles = Int[]

    normal = (p1-p3) × (p2-p3)
    normal /= norm(normal)

    height = (center - p1) ⋅ normal
    plane_center = center - height * normal

    inner_radius = outer_radius = zero(eltype(p1))
    plane_inner_radius = plane_outer_radius = zero(eltype(p1))

    IntegrationPath(
      center,
      inner_radius,
      outer_radius,
      segments,
      arcs,
      circles,
      normal,
      height,
      plane_center,
      plane_inner_radius,
      plane_outer_radius,
    )
end

mutable struct WorkSpace{V}
  inner_in::Vector{V}
  inner_out::Vector{V}
  outer_in::Vector{V}
  outer_out::Vector{V}
  inner_in_edge::Vector{Int}
  inner_out_edge::Vector{Int}
  outer_in_edge::Vector{Int}
  outer_out_edge::Vector{Int}
  segments::Vector{Tuple{V,V}}
  arcs::Vector{Tuple{V,V,Int}}
  circles::Vector{Int}
end

function contour(p1, p2, p3, center, inner_radius, outer_radius)
  ws = workspace(typeof(p1))
  contour!(p1, p2, p3, center, inner_radius, outer_radius, ws)
end

function workspace(V::Type)
  ws = WorkSpace(
    Vector{V}(undef,3), Vector{V}(undef,3), Vector{V}(undef,3), Vector{V}(undef,3),
    Vector{Int}(undef,3), Vector{Int}(undef,3), Vector{Int}(undef,3), Vector{Int}(undef,3),
    Vector{Tuple{V,V}}(undef,0), Vector{Tuple{V,V,Int}}(undef,0), Vector{Int}(undef,0)
  )
  sizehint!(ws.segments,8)
  sizehint!(ws.arcs,8)
  sizehint!(ws.circles,2)
  return ws
end

function contour!(p1, p2, p3, center, inner_radius, outer_radius, ws)

  deci = eltype(p1)
  Vertex = typeof(p1)

  #t_inner = zeros(deci,2)
  #t_outer = zeros(deci,2)

  #inner_in  = Vector{Vertex}(3)
  #inner_out = Vector{Vertex}(3)
  #outer_in  = Vector{Vertex}(3)
  #outer_out = Vector{Vertex}(3)

  #inner_in_edge  = zeros(Int,3)
  #inner_out_edge = zeros(Int,3)
  #outer_in_edge  = zeros(Int,3)
  #outer_out_edge = zeros(Int,3)
  fill!(ws.inner_in_edge, 0)
  fill!(ws.inner_out_edge, 0)
  fill!(ws.outer_in_edge, 0)
  fill!(ws.outer_out_edge, 0)

  is_not_empty = number_of_segments = number_of_arcs = number_of_circles = 0
  inner = outer = inner_inc = inner_outc = outer_inc = outer_outc = 0

  # centra en stralen berekenen
  inner_radius = (inner_radius > 0.0) ? inner_radius : 0.0
  outer_radius = (outer_radius > 0.0) ? outer_radius : 0.0

  # Vertex normal = f.normal()
  normal = (p1-p3) × (p2-p3)
  normal /= norm(normal)

  d = (center - p1) ⋅ normal
  absd = abs(d)
  plane_center = center - dot(normal, center - p1) * normal

  discr = inner_radius*inner_radius - d * d
  plane_inner_radius = (discr > 0.0) ? sqrt(discr) : -1.0
  discr = outer_radius*outer_radius - d * d
  plane_outer_radius = (discr > 0.0) ? sqrt(discr) : -1.0

  #segments = Vector{Tuple{Vertex,Vertex}}()
  #arcs = Vector{Tuple{Vertex,Vertex,Int}}()
  #circles = Vector{Int}()
  resize!(ws.segments,0)
  resize!(ws.arcs,0)
  resize!(ws.circles,0)

  # loop over de edges om segmenten te vinden
  first = zero(Vertex)
  last  = zero(Vertex)
  F = (p1,p2,p3)
  L = (p2,p3,p1)
  for i in 1:3
    first, last = F[i], L[i]
    z_inner, t_inner = incidenceLineWithSphere(center, first, last, inner_radius)
    z_outer, t_outer = incidenceLineWithSphere(center, first, last, outer_radius)
    code = 10 * z_outer + z_inner

    if code == 00
      d1 = norm(first - center)
      d2 = norm(last - center)
      D = 0.5 * (d1 + d2)
      if !(D < inner_radius || D > outer_radius)
        push!(ws.segments, (first,last))
      end

    elseif code == 01
      t1 = t_inner[1]
      v1 = (1-t1) * first + t1 * last
      d1 = norm(first - center)
      d2 = norm(last - center)
      if d1 > d2
        push!(ws.segments, (first, v1))
        inner_inc += 1
        ws.inner_in[inner_inc] = v1
        ws.inner_in_edge[inner_inc] = i
      else
        push!(ws.segments, (v1, last))
        inner_outc += 1
        ws.inner_out[inner_outc] = v1
        ws.inner_out_edge[inner_outc] = i
      end

    elseif code == 02
      t1 = t_inner[1]
      t2 = t_inner[2]
      v1 = (1-t1) * first +  t1 * last
      inner_inc += 1
      ws.inner_in[inner_inc] = v1
      ws.inner_in_edge[inner_inc] = i
      v2 = (1-t2) * first + t2 * last
      inner_outc += 1
      ws.inner_out[inner_outc] = v2
      ws.inner_out_edge[inner_outc] = i
      push!(ws.segments, (first,v1))
      push!(ws.segments, (v2,last))

    elseif code == 10
      t1 = t_outer[1]
      v1 = (1-t1) * first +  t1 * last
      d1 = norm(first - center)
      d2 = norm(last - center)
      if d1 < d2
          push!(ws.segments, (first, v1))
          outer_inc += 1
          ws.outer_in[outer_inc] = v1
          ws.outer_in_edge[outer_inc] = i
      else
          push!(ws.segments, (v1, last))
          outer_outc += 1
          ws.outer_out[outer_outc] = v1
          ws.outer_out_edge[outer_outc] = i
      end

    elseif code == 11
      t1 = t_inner[1]
      t2 = t_outer[1]
      v1 = (1-t1) * first + t1 * last
      v2 = (1-t2) * first + t2 * last
      if t1 < t2
        push!(ws.segments, (v1, v2))
        inner_outc += 1
        ws.inner_out[inner_outc] = v1
        ws.inner_out_edge[inner_outc] = i
        outer_inc += 1
        ws.outer_in[outer_inc] = v2
        ws.outer_in_edge[outer_inc] = i
      else
        push!(ws.segments, (v2, v1))
        inner_inc += 1
        ws.inner_in[inner_inc] = v1
        ws.inner_in_edge[inner_inc] = i
        outer_outc += 1
        ws.outer_out[outer_outc] = v2
        ws.outer_out_edge[outer_outc] = i
      end

    elseif code == 12
      t1 = t_inner[1]
      t2 = t_inner[2]
      t3 = t_outer[1]
      v1 = (1-t1) * first + t1 * last
      v2 = (1-t2) * first + t2 * last
      v3 = (1-t3) * first + t3 * last
      if t1 < t3
        push!(ws.segments, (first, v1))
        push!(ws.segments, (v2, v3))
        inner_inc += 1
        ws.inner_in[inner_inc] = v1
        ws.inner_in_edge[inner_inc] = i
        inner_outc += 1
        ws.inner_out[inner_outc] = v2
        ws.inner_out_edge[inner_outc] = i
        outer_inc += 1
        ws.outer_in[outer_inc] = v3
        ws.outer_in_edge[outer_inc] = i
      else
        push!(ws.segments, (v3, v1))
        push!(ws.segments, (v2, last))
        inner_inc += 1
        ws.inner_in[inner_inc] = v1
        ws.inner_in_edge[inner_inc] = i
        inner_outc += 1
        ws.inner_out[inner_outc] = v2
        ws.inner_out_edge[inner_outc] = i
        outer_outc += 1
        ws.outer_out[outer_outc] = v3
        ws.outer_out_edge[outer_outc] = i
      end

    elseif code == 20
      t1 = t_outer[1]
      t2 = t_outer[2]
      v1 = (1-t1) * first + t1 * last
      outer_outc += 1
      ws.outer_out[outer_outc] = v1
      ws.outer_out_edge[outer_outc] = i
      v2 = (1-t2) * first + t2 * last
      outer_inc += 1
      ws.outer_in[outer_inc] = v2
      ws.outer_in_edge[outer_inc] = i
      push!(ws.segments, (v1, v2))

    elseif code == 22
      t1 = t_inner[1]
      t2 = t_inner[2]
      t3 = t_outer[1]
      t4 = t_outer[2]
      v1 = (1-t1) * first + t1 * last
      v2 = (1-t2) * first + t2 * last
      v3 = (1-t3) * first + t3 * last
      v4 = (1-t4) * first + t4 * last
      inner_inc += 1
      ws.inner_in[inner_inc] = v1
      ws.inner_in_edge[inner_inc] = i
      inner_outc += 1
      ws.inner_out[inner_outc] = v2
      ws.inner_out_edge[inner_outc] = i
      outer_outc += 1
      ws.outer_out[outer_outc] = v3
      ws.outer_out_edge[outer_outc] = i
      outer_inc += 1
      ws.outer_in[outer_inc] = v4
      ws.outer_in_edge[outer_inc] = i
      push!(ws.segments, (v3, v1))
      push!(ws.segments, (v2, v4))
    end

    #@show code

  end # end for i

  # bereken het aantal snijpunten met de twee cirkels
  inner = inner_inc + inner_outc
  outer = outer_inc + outer_outc

  #@show inner_inc inner_outc
  #@show outer_inc outer_outc

  try
  @assert inner_inc == inner_outc
  @assert outer_inc == outer_outc
  catch
      # @show inner_inc inner_outc
      # @show outer_inc outer_outc
      # @show p1, p2, p3
      # @show center, inner_radius, outer_radius
      error("inconsistent contour construction!")
  end

  # construct the inner arcs
  if inner == 0
  elseif inner == 2
    push!(ws.arcs, (ws.inner_in[1], ws.inner_out[1], CLOCKWISE))
  elseif inner == 4
    if ws.inner_in_edge[1] == ws.inner_out_edge[1]
      push!(ws.arcs, (ws.inner_in[1],ws.inner_out[2],CLOCKWISE))
      push!(ws.arcs, (ws.inner_in[2], ws.inner_out[1], CLOCKWISE))
    elseif ws.inner_in_edge[1] == ws.inner_out_edge[2]
      push!(ws.arcs, (ws.inner_in[1],ws.inner_out[1],CLOCKWISE))
      push!(ws.arcs, (ws.inner_in[2], ws.inner_out[2],CLOCKWISE))
    else
      ed1 = ws.inner_in_edge[1]
      ed2 = mod1(ed1-1, 3)
      if ws.inner_out_edge[1] == ed2
        # This branch -I think- cannot be reached
        push!(ws.arcs, (ws.inner_in[1], ws.inner_out[1], CLOCKWISE))
        push!(ws.arcs, (ws.inner_in[2], ws.inner_out[2], CLOCKWISE))
      else
        push!(ws.arcs, (ws.inner_in[1], ws.inner_out[2], CLOCKWISE))
        push!(ws.arcs, (ws.inner_in[2], ws.inner_out[1], CLOCKWISE))
      end
    end
  elseif inner == 6
    for i in 1:3
      ed1 = ws.inner_in_edge[i]
      ed2 = mod1(ed1-1,3)
      for j in 1:3
        if ed2 == ws.inner_out_edge[j]
          push!(ws.arcs, (ws.inner_in[i], ws.inner_out[j], CLOCKWISE))
        end
      end
    end
  end

  # construct the outer arcs
  if outer == 0
  elseif outer == 2
    push!(ws.arcs, (ws.outer_in[1], ws.outer_out[1], COUNTERCLOCKWISE))
  elseif outer == 4
    if ws.outer_in_edge[1] == ws.outer_out_edge[1]
      push!(ws.arcs, (ws.outer_in[1], ws.outer_out[2], COUNTERCLOCKWISE))
      push!(ws.arcs, (ws.outer_in[2], ws.outer_out[1], COUNTERCLOCKWISE))
    elseif ws.outer_in_edge[1] == ws.outer_out_edge[2]
      # This branch cannot be reached
      push!(ws.arcs, (ws.outer_in[1], ws.outer_out[1], COUNTERCLOCKWISE))
      push!(ws.arcs, (ws.outer_in[2], ws.outer_out[2], COUNTERCLOCKWISE))
    else
      ed1 = ws.outer_in_edge[1]
      ed2 = mod1(ed1+1,3)
      if ws.outer_out_edge[1] == ed2
        push!(ws.arcs, (ws.outer_in[1], ws.outer_out[1], COUNTERCLOCKWISE))
        push!(ws.arcs, (ws.outer_in[2], ws.outer_out[2], COUNTERCLOCKWISE))
      else
        push!(ws.arcs, (ws.outer_in[1], ws.outer_out[2], COUNTERCLOCKWISE))
        push!(ws.arcs, (ws.outer_in[2], ws.outer_out[1], COUNTERCLOCKWISE))
      end
    end
  elseif outer == 6
    for i in 1:3
      ed1 = ws.outer_in_edge[i]
      ed2 = mod1(ed1+1, 3)
      for j in 1:3
        if ed2==ws.outer_out_edge[j]
          push!(ws.arcs, (ws.outer_in[i], ws.outer_out[j], COUNTERCLOCKWISE))
        end
      end
    end
  end

  # construct the circle contributions
  if inside(plane_center,p1,p2,p3,normal)
    d1 = distancetoline(plane_center, p1, p2)
    d2 = distancetoline(plane_center, p2, p3)
    d3 = distancetoline(plane_center, p3, p1)
    if d1>plane_inner_radius && d2>plane_inner_radius && d3>plane_inner_radius && plane_inner_radius>0
      push!(ws.circles, CLOCKWISE)
    end
    if d1>plane_outer_radius && d2>plane_outer_radius && d3>plane_outer_radius && plane_outer_radius>0
      push!(ws.circles, COUNTERCLOCKWISE)
    end
  end

  IntegrationPath(
    center,
    inner_radius,
    outer_radius,
    ws.segments,
    ws.arcs,
    ws.circles,
    normal,
    d,
    plane_center,
    plane_inner_radius,
    plane_outer_radius
  )
end

function inside(v,p1,p2,p3,n)
  dot((p2-p1)×(v-p1), n) <= 0 && return false
  dot((p3-p2)×(v-p2), n) <= 0 && return false
  dot((p1-p3)×(v-p3), n) <= 0 && return false
  return true
end

function distancetoline(p,a,b)
  t = b-a
  t /= norm(t)
  u = p-a
  norm(u - dot(u,t)*t)
end

"""
  r, t = incidenceLineWithSphere(v, first, last, r)

Returns the number of crossings and their barycentric coordinates for a circle
with center v and radius r and a segment [a,b]. If the circle is tangent to the
segment, zero crossings are reported.
"""
function incidenceLineWithSphere(v, first, last, r)

    ϵ = eps(typeof(r)) * 1e3
    z = zero(typeof(r))

    r<=0 && return 0, (z,z)

    a = dot((last - first) , (last - first))
    b =  2 * dot( (first - v) , (last - first))
    c = dot((first - v) , (first - v)) - r*r
    d = b * b - 4 * a * c
    d < ϵ && return 0, (z,z)
    f = sqrt(d)

    t1 = (-b - f) / (2 * a)
    p1 = (0 < t1 <= 1)

    t2 = (-b + f) / (2 * a)
    p2 = (0 < t2 <= 1)

    if p1
        if p2
            return 2, (t1,t2)
        else
            return 1, (t1,z)
        end
    else
        if p2
            return 1, (t2,z)
        else
            return 0, (z,z)
        end
    end

    return 0, (z,z)
end

ext: jl


text:
struct Raw <: FileFormat end

# Raw Parser
# ==========

function _get_seqs(io::Union{IO, AbstractString})
    SEQS = String[]
    for line::String in lineiterator(io)
        push!(SEQS, line)
    end
    SEQS
end

function Base.parse(io::Union{IO, AbstractString},
                    format::Type{Raw},
                    output::Type{T};
                    deletefullgaps::Bool=true)::T where T <: MSAMatrix
    SEQS = _get_seqs(io)
    _strings_to_msa(T, SEQS, deletefullgaps)
end

function Base.parse(io::Union{IO, AbstractString},
                    format::Type{Raw},
                    output::Type{AnnotatedMultipleSequenceAlignment};
                    generatemapping::Bool=false,
                    deletefullgaps::Bool=true,
                    keepinserts::Bool=false)
    SEQS = _get_seqs(io)
    annot = Annotations()
    if keepinserts
        _keepinserts!(SEQS, annot)
    end
    if generatemapping
        MSA, MAP = _to_msa_mapping(SEQS) # It checks sequence lengths
        setannotfile!(annot, "NCol", string(size(MSA,2)))
        setannotfile!(annot, "ColMap", join(vcat(1:size(MSA,2)), ','))
        for i in 1:length(SEQS)
            setannotsequence!(annot, string(i), "SeqMap", MAP[i])
        end
    else
        MSA = NamedArray(convert(Matrix{Residue}, SEQS)) # It checks sequence lengths
    end
    msa = AnnotatedMultipleSequenceAlignment(MSA, annot)
    if deletefullgaps
        deletefullgapcolumns!(msa)
    end
    msa
end

function Base.parse(io::Union{IO, AbstractString},
                    format::Type{Raw},
                    output::Type{MultipleSequenceAlignment};
                    deletefullgaps::Bool=true)
    msa = parse(io, Raw, NamedResidueMatrix{Array{Residue,2}}; deletefullgaps=deletefullgaps)
    MultipleSequenceAlignment(msa)
end

function Base.parse(io::Union{IO, AbstractString}, format::Type{Raw};
                    deletefullgaps::Bool=true)
    parse(io, Raw, AnnotatedMultipleSequenceAlignment; deletefullgaps=deletefullgaps)
end

# Print Raw
# =========

function Base.print(io::IO, msa::AbstractMatrix{Residue}, format::Type{Raw})
    for i in 1:nsequences(msa)
        println(io, stringsequence(msa, i))
    end
end

Base.print(msa::AbstractMatrix{Residue}, format::Type{Raw}) = print(stdout, msa, Raw)


text:
module var"#Internal"
public(x::String) = false
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @interpret sum(rand(10))
    expr = quote
        public(x::Integer) = true
        module Private
            private(y::String) = false
        end
        const threshold = 0.1
    end
    for (mod, ex) in ExprSplitter(var"#Internal", expr)
        frame = Frame(mod, ex)
        debug_command(frame, :c, true)
    end
end


text:
###################
# DISPLAY METHODS #
###################
using LaTeXStrings
plot_seq(seq::Sequence) = begin
	idx = ["Gx" "Gy" "Gz"]
	M, N = size(seq.GR)
	O, _ = size(seq.RF)
	G = [j <= M ? seq.GR[j,floor(Int,i/2)+1].A : 0. for i=0:2*N-1, j=1:3]
	R = [seq.RF[j,floor(Int,i/2)+1].A for i=0:2*N-1, j=1:O]
	D = [seq.DAC[floor(Int,i/2)+1].N > 0 for i=0:2*N-1]
	T = [seq.GR[1,i].T for i=1:N]
	t = [sum(T[1:i]) for i=1:N]
	t = [t[floor(Int,i/2)+1] for i=0:2*N-1]
	t = [0; t[1:end-1]]

	l = PlotlyJS.Layout(;yaxis_title="G [mT/m]",#, hovermode="x unified", #title="Sequence",
			xaxis_title="t [ms]",height=300, modebar=attr(orientation="v"),
			legend=attr(orientation="h",yanchor="bottom",xanchor="right",y=1.02,x=1),
			xaxis=attr(
				range=[0.,10],
				rangeslider=attr(visible=true),
				rangeselector=attr(
					buttons=[
						attr(count=1,
						label="1m",
						step=10,
						stepmode="backward"),
						attr(step="all")
						]
					)
				)
			)
	p = [PlotlyJS.scatter() for j=1:(3+O+1)]
	#GR
	for j=1:3
		p[j] = PlotlyJS.scatter(x=t*1e3, y=G[:,j]*1e3,name=idx[j],line_shape="hv",hovertemplate="%{y:.1f} mT/m")
	end
	#RF
	for j=1:O
		p[j+3] = PlotlyJS.scatter(x=t*1e3, y=abs.(R[:,j])*1e6,name="RF_$j",line_shape="hv",hovertemplate="%{y:.1f} μT")
	end
	#DAC
	p[O+3+1] = PlotlyJS.scatter(x=t*1e3, y=D*1., name="DAC")
	PlotlyJS.plot(p, l)#, options=Dict(:displayModeBar => false))
end

plot_grads_moments(seq::Sequence, idx::Int=1; title="", mode="quick") = begin
	if mode == "quick"
		plotter = PlotlyJS.scattergl
	else
		plotter = PlotlyJS.scatter
	end
	names = ["Gx", "Gy", "Gz"]
	M, N = size(seq.GR)
	G = [seq.GR[j,floor(Int,i/2)+1].A for i=0:2*N-1, j=1:M]
	T = [seq.GR[1,i].T for i=1:N]
	t = [sum(T[1:i]) for i=1:N]
	t = [t[floor(Int,i/2)+1] for i=0:2*N-1]
	t = [0; t[1:end-1]]

	M0, M1, M2 = get_M0_M1_M2(seq)
	M0t = [M0(t)[idx] for t=t][:]; M0max = maximum(abs.(M0t))*1.1e6
	M1t = [M1(t)[idx] for t=t][:]; M1max = maximum(abs.(M1t))*1.1e9
	M2t = [M2(t)[idx] for t=t][:]; M2max = maximum(abs.(M2t))*1.1e12
	Gmax = maximum(abs.(G))*1.1e3

	p = [plotter() for k = 1:4]

	l = PlotlyJS.Layout(;
	title=title,
	xaxis=attr(domain = [0, .75]),
	yaxis=attr(title=attr(text="G [mT/m]", standoff=0),
		side="left", tickfont=attr(color="#636efa"), titlefont=attr(color="#636efa"),range=[-Gmax, Gmax]),
	yaxis2=attr(title=attr(text="M0 [mT/m⋅ms]", standoff=1, position=1), showgrid=false,
		anchor="free", overlaying="y",position=0.75, side="right",
		tickfont=attr(color="#ef553b"), titlefont=attr(color="#ef553b"), range=[-M0max, M0max]),
	yaxis3=attr(title=attr(text="M1 [mT/m⋅ms²]", standoff=1), showgrid=false,
		anchor="free",overlaying="y",position=0.75+.25/3,side="right",
		tickfont=attr(color="#45d9b2"), titlefont=attr(color="#45d9b2"), range=[-M1max, M1max]),
	yaxis4=attr(title=attr(text="M2 [mT/m⋅ms³]", standoff=1), showgrid=false,
		anchor="free",overlaying="y",position=0.75+.5/3,side="right",
		tickfont=attr(color="#b373fa"), titlefont=attr(color="#b373fa"), range=[-M2max, M2max]),
	xaxis_title="t [ms]", height=300)
	p[1] = plotter(x=t*1e3, y=G[:,idx]*1e3, name=names[idx] , line_shape="hv")
	p[2] = plotter(x=t*1e3, y=M0t*1e6, name="M0", yaxis="y2",line=attr(dash="dash"))
	p[3] = plotter(x=t*1e3, y=M1t*1e9, name="M1", yaxis="y3",line=attr(dash="dash"))
	p[4] = plotter(x=t*1e3, y=M2t*1e12, name="M2", yaxis="y4",line=attr(dash="dash"))
	PlotlyJS.plot(p, l)
end

# plot_Phantom(obj::Phantom,filename::String) = begin
# 	# Phantom
# 	p1 = heatmap(obj.x*1e2,obj.y*1e2,obj.ρ,aspect_ratio=:equal)#,clims=(0,maximum(T2[:])))
# 	xaxis!("x [cm]"); yaxis!("y [cm]")
# 	title!("Proton density")
# 	p2 = heatmap(obj.x*1e2,obj.y*1e2,obj.T2*1e3,aspect_ratio=:equal)#,clims=(0,maximum(T2[:])))
# 	xaxis!("x [cm]"); yaxis!("y [cm]")
# 	title!("T2 [ms]")
# 	p3 = heatmap(obj.x*1e2,obj.y*1e2,obj.Δw/2π,aspect_ratio=:equal)#,clims=(0,maximum(T2[:])))
# 	xaxis!("x [cm]"); yaxis!("y [cm]")
# 	title!("df [Hz]")
# 	P(i,j) = rotz(Dθ[i,j])[1:2,1:2]; D(i,j) = [obj.Dλ1[i,j] 0;0 obj.Dλ2[i,j]]
# 	nx = [1;0]; ny = [0;1]
# 	Dx = [nx'*P(i,j)'*D(i,j)*P(i,j)*nx for i=1:size(obj.Dλ1,1),j=1:size(obj.Dλ1,2)]
# 	Dy = [ny'*P(i,j)'*D(i,j)*P(i,j)*ny for i=1:size(obj.Dλ1,1),j=1:size(obj.Dλ1,2)]
# 	p4 = heatmap(obj.x*1e2,obj.y*1e2,Dx*1e12,aspect_ratio=:equal)#,clims=(0,maximum(T2[:])))
# 	xaxis!("x [cm]"); yaxis!("y [cm]")
# 	title!("Dx [um2/s]")
# 	p5 = heatmap(obj.x*1e2,obj.y*1e2,Dy*1e12,aspect_ratio=:equal)#,clims=(0,maximum(T2[:])))
# 	xaxis!("x [cm]"); yaxis!("y [cm]")
# 	title!("Dy [um2/s]")
# 	p = plot(p1,p2,p3,p4,p5,size=(1300,500),layout=@layout [a b c; d e])
# 	savefig(p,filename)
# end
# plot_sim_res(obj::Phantom,SEQ::Array{Sequence},S::Array{ComplexF64},
# 	t::Array{Float64,2},filename::String,t_k0::Float64,Ga::Float64) = begin
# 	T2m = sum(obj.T2.*obj.ρ)/sum(obj.ρ)
# 	b, n = get_bvalue(SEQ[1])
# 	println("Sequence with b-value: "*string(round(b*1e-6,digits=2))*" s/mm2") # s/mm2
# 	P(i) = rotz(obj.Dθ[i])[1:2,1:2]; D(i) = [obj.Dλ1[i] 0;0 obj.Dλ2[i]]
# 	Deq = [n'*P(i)'*D(i)*P(i)*n for i=1:prod(size(obj.Dλ1))]
# 	Dm = sum(Deq.*obj.ρ)/sum(obj.ρ) #MODIFYY

# 	p = plot_grads(SEQ,t,t_k0,Ga)
# 	p3 = plot([t_k0; t_k0]*1e3,[0; 1.1],linewidth=0.25,color=:black,label="k=0")
# 	plot!(t[:]*1e3,abs.(S),label="|S|")
# 	plot!(t[:]*1e3,exp.(-t[:]/T2m),linestyle=:dash,label="exp(-t/T2)")
# 	# exp(-b*D) <-> exp(-4*π^2*(Δ-δ/3)*q'*D*q)
# 	plot!(t[:]*1e3,exp.(-t[:]/T2m.-b*Dm).*(t[:].>=(Δ+δ)),linestyle=:dash,label="exp(-t/T2-bD)")

# 	xlabel!("Time [ms]")
# 	ylabel!("Signal [a.u]")
# 	xlims!((minimum(t),minimum(t).+dur(sum(SEQ))).*1e3)
# 	p = plot(p[1],p[2],p3,size=(800,600),layout=@layout [a ; b; c])
# 	savefig(p,filename)
# end
# plot_ksapce_trajectory(ACQ::Sequence,t::Array{Float64,2},filename::String) = begin
# 	k = get_designed_kspace(ACQ)
# 	p = plot(k[:,1],k[:,2],legend=:none,aspect_ratio=:equal)
# 	k = get_actual_kspace(ACQ,t)
# 	scatter!(k[:,1],k[:,2],legend=:none,markersize=1)
# 	xlabel!("kx [1/m]"); ylabel!("ky [1/m]")
# 	savefig(p,filename)
# end
# plot_recon(kdata::Array{ComplexF64},rec::Array{ComplexF64},
# 	Δx_pix::Float64,Δy_pix::Float64,filename::String,title::String) = begin
# 	Nx, Ny = size(rec)
# 	xr = -Δx_pix*(Nx-1)/2:Δx_pix:Δx_pix*(Nx-1)/2
# 	yr = -Δy_pix*(Ny-1)/2:Δy_pix:Δy_pix*(Ny-1)/2
# 	Wx, Wy = 1/Δx_pix, 1/Δy_pix
# 	kx = range(-Wx/2,stop=Wx/2,length=Nx)
# 	ky = range(-Wy/2,stop=Wy/2,length=Ny)
# 	p1 = heatmap(kx,ky,log.(abs.(kdata).+1),aspect_ratio=:equal,
# 				legend=:none,size=(400,400))
# 	xaxis!("kx [1/m]"); yaxis!("ky [1/m]")
# 	title!("k-space")
# 	p2 = heatmap(xr*1e2,yr*1e2,abs.(rec),aspect_ratio=:equal,
# 				legend=:none,size=(400,400))
# 	xaxis!("x [cm]"); yaxis!("y [cm]")
# 	title!(title)
# 	#savefig(p1,filename*"_ksp.pdf")
# 	savefig(p2,filename*".pdf")
# end
# plot_Eq(vx,vy,S,S0,SEQ) = begin
# 	Nq, Nθ = size(SEQ)
# 	rs = range(0,stop=1,length=Nq)
# 	θs = range(0,stop=π,length=Nθ); θs = [θs[1:end-1]; π.+θs[1:end-1]];
#     Eq(vx,vy) = begin
#         Eq = [S[i,j][vx,vy] for i=1:Nq,j=1:Nθ]
#         Eq = [Eq[:,1:end-1] Eq[:,1:end-1]]
#     end
#     c = RGB{Float64}(1,1,1)
#     α = real.(abs.(S0[vx,vy]))/maximum(real.(abs.(S0)))
#     pyplot()
#     hm = heatmap(θs,rs,Eq(vx,vy)*(α<.2 ? 0 : 1),proj=:polar,aspect_ratio=:equal,
#     	legend=:none,grid=false,xticks=:none,yticks=:none,
#     	background_color_subplot=α*c)
# end
# plot_Pr(vx,vy,S,S0,SEQ) = begin
# 	Nq, Nθ = size(SEQ)
# 	rs = range(0,stop=1,length=10)*25e-6
# 	θs = range(0,stop=2π,length=10)
#     Pr(vx,vy) = begin
#         Eq = [S[i,j][vx,vy] for i=1:Nq,j=1:Nθ]
# 		b = [get_bvalue(SEQ[i,j])[1] for i=1:Nq,j=1:Nθ]
# 		n = [get_bvalue(SEQ[i,j])[2] for i=1:Nq,j=1:Nθ]
# 		B = [-b[i,j]*[n[i,j][1]^2
# 					2*n[i,j][1]*n[i,j][2]
# 					  n[i,j][2]^2] for i=1:Nq,j=1:Nθ]
# 		B = [B[i][j] for i=1:Nq*Nθ,j=1:3] #reshape
# 		y = log.(abs.(Eq))[:]
# 		# exp(-b*n'D n) <-> exp(-4*π^2*(Δ-δ/3)*q'*D*q)
# 		# min D || -b*n'D*n - log(S) ||_2 + λ ||D||_2
# 		# min D || B D - y ||_2 + λ ||D||_2 => D = (B'B+λI)^-1 B' y
# 		λ = 0
# 		Id = Matrix{Float64}(I,3,3); D = [0;0;0]
# 		for n = 0:2 #Tikhonov regularized and iterativly weighted
# 			W = n==0 ? diagm(0=>Eq[:].^2) : diagm(0=>exp.(2*B*D))
# 			D = (B'*W*B + λ*Id)^-1*B'*W*y
# 		end
# 		# Diffusion propagator
# 		D_inv = [D[1] D[2];
# 				 D[2] D[3]]^-1
# 		pr = exp.(-[([r*cos(θ) r*sin(θ)]*D_inv*[r*cos(θ);r*sin(θ)])[1] for r=rs, θ=θs])
#     end
#     c = RGB{Float64}(1,1,1)
#     α = real.(abs.(S0[vx,vy]))/maximum(real.(abs.(S0)))
#     pyplot()
#     hm = heatmap(θs,rs,Pr(vx,vy)*(α<.2 ? 0 : 1),proj=:polar,aspect_ratio=:equal,
#     	legend=:none,grid=false,xticks=:none,yticks=:none,
#     	background_color_subplot=α*c)
# end


text:
## Phase diagram (superconducting, armchair, optim)

function runfig3()
    p = Params(Ln = 1440, W = 1440, Ls = 0, scale = 40,
        λ = 5, α = 0, EZ = SA[0, 1e-3, 0], Δ = 0.3, d = 0, τ = 1);
    ha = nanoribbonSA(p);
    hz = nanoribbonSZ(p);
    pha = parametric(ha, @onsite!((o; dμ = 0, dEZy = 0) -> o - dμ * σ0τz +  dEZy*σyτ0))
    phz = parametric(hz, @onsite!((o; dμ = 0, dEZy = 0) -> o - dμ * σ0τz +  dEZy*σyτ0))
    mata = similarmatrix(pha, Quantica.flatten)
    matz = similarmatrix(phz, Quantica.flatten)
    dμs = range(-2.5, 2.5, 100)
    boundarySA = @showprogress [optimize(fA(dμ, mata, pha), 0.00001, 2).minimizer
        for dμ in dμs]
    boundarySZ = @showprogress [optimize(fZ(dμ, matz, phz), 0.00001, 2).minimizer
        for dμ in dμs]
    return savephasediagram("fig3", p, boundarySA, boundarySZ, dμs)
end

function gapA( mat, ph; kw...)
    bloch!(mat, ph(; kw...), 0)
    λ = maximum(real.(first(eigs(mat; sigma = 0.00001, nev = 4))))
    return λ
end

function gapZ(mat, ph; kw...)
    k0 = 2.5692321  # that's the Bloch phase where normal zigzag bands cross zero
    bloch!(mat, ph(; kw...), k0)
    λ = maximum(real.(first(eigs(mat; sigma = 0.00001, nev = 4))))
    return λ
end

fA(dμ, mat, ph) = x -> gapA(mat, ph; dEZy = x, dμ)
fZ(dμ, mat, ph) = x -> gapZ(mat, ph; dEZy = x, dμ)

max_stars_repo_name: fernandopenaranda/MBSinBLG

max_issues_repo_name: fernandopenaranda/MBSinBLG

max_forks_repo_name: fernandopenaranda/MBSinBLG


text:
module TestUnivariateFiniteArray

using Test
using MLJBase
using CategoricalArrays
import Distributions:pdf, support
import Distributions
using StableRNGs
import Random

rng = StableRNG(111)
n   = 10
c   = 3

@testset "constructing UnivariateFiniteArray objects" begin

    probs  = rand(rng, n)
    supp = ["class1", "class2"]

    @test_throws DomainError UnivariateFinite(supp, probs, pool=missing)
    u = UnivariateFinite(supp, probs, pool=missing, augment=true)
    @test length(u) == n
    @test size(u) == (n,)
    @test pdf.(u, "class2") ≈ probs

    # autosupport:
    @test_throws DomainError UnivariateFinite(probs, pool=missing)
    u = UnivariateFinite(probs, pool=missing, augment=true)
    @test length(u) == n
    @test size(u) == (n,)
    @test pdf.(u, :class_2) ≈ probs
    probs = probs ./ sum(probs)
    u = UnivariateFinite(probs, pool=missing)
    @test u isa UnivariateFinite
    @test pdf(u, :class_1) == probs[1]
    probs = rand(10, 2)
    probs = probs ./ sum(probs, dims=2)
    u = UnivariateFinite(probs, pool=missing)
    @test length(u) == 10
    u.scitype == Multiclass{2}
    pdf.(u, :class_1) == probs[:, 1]
    u = UnivariateFinite(probs, pool=missing, augment=true)
    @test length(u) == 10
    u.scitype == Multiclass{3}
    pdf.(u, :class_2) == probs[:, 1]

    probs = [-1,0,1]
    @test_throws(DomainError,
                 UnivariateFinite(probs, pool=missing, augment=true))

    v = categorical(1:3)
    @test_logs((:warn, r"Ignoring"),
               UnivariateFinite(v[1:2], rand(3), augment=true, pool=missing))
    @test_logs((:warn, r"Ignoring"),
               UnivariateFinite(v[1:2], rand(3), augment=true, ordered=true))

end

@testset "get and set" begin
    # binary
    s = rand(rng, n)
    u = UnivariateFinite([:yes, :no], s, augment=true, pool=missing)

    @test u[1] isa UnivariateFinite
    v = u[3:4]
    @test v isa UnivariateFiniteArray
    @test v[1] ≈ u[3]

    # set:
    u[1] = UnivariateFinite([:yes, :no], [0.1, 0.9], pool=u[1])
    @test pdf(u[1], :yes) == 0.1
    @test pdf(u[1], :no) == 0.9

    # multiclass
    P   = rand(rng, n, c)
    P ./= sum(P, dims=2)
    u   = UnivariateFinite(P, pool=missing)

    @test u[1] isa UnivariateFinite
    v = u[3:4]
    @test v isa UnivariateFiniteArray
    @test v[1] ≈ u[3]

    # set:
    s = Distributions.support(u)
    scalar = UnivariateFinite(s, P[end,:])
    u[1] = scalar
    @test u[1] ≈ u[end]
    @test pdf.(u, s[1])[2:end] == P[2:end,1]
end

n = 10
P = rand(n);
all_classes = categorical([:no, :yes], ordered=true)
u = UnivariateFinite(all_classes, P, augment=true)

# next is not specific to `UnivariateFiniteArray` but is for any
# abstract array with eltype `UnivariateFinite`:
@testset "piratical pdf" begin
    @test pdf(u, [:yes, :no]) == hcat(P, 1 .- P)
    @test pdf(u, reverse(all_classes)) == hcat(P, 1 .- P)
    @test pdf([u...], [:yes, :no]) == hcat(P, 1 .- P)
    @test pdf([u...], all_classes) == hcat(1 .- P, P)
end

@testset "broadcasting: pdf.(uni_fin_arr, scalar) " begin
    @test pdf.(u,:yes) == P
    @test pdf.(u,all_classes[2]) == P

    # check unseen probablities are a zero *array*:
    v = categorical(1:4)
    probs = rand(3)
    u2 = UnivariateFinite(v[1:2], probs, augment=true)
    @test pdf.(u2, v[3]) == zeros(3)
end

@testset "broadcasting: pdf.(uni_fin_arr, array_same_shape)" begin
    v = rand(classes(u), n)
    @test broadcast(pdf, u, v) == [pdf(u[i], v[i]) for i in 1:length(u)]
    @test broadcast(pdf, u, get.(v)) ==
        [pdf(u[i], v[i]) for i in 1:length(u)]
end

@testset "broadcasting: check indexing in `getter((cv, i), dtype)` see PR#375" begin
    c  = categorical([0,1,1])
    d = UnivariateFinite(c[1:1], [1 1 1]')
    v = categorical([0,1,1,1])
    @test broadcast(pdf, d, v[2:end]) == [0,0,0]
end

@testset "_getindex" begin
   @test MLJBase._getindex(collect(1:4), 2, Int64) == 2
   @test MLJBase._getindex(nothing, 2, Int64) == zero(Int64)
end

@testset "broadcasting mode" begin
    # binary
    rng = StableRNG(668)
    probs = rand(rng, n)
    u = UnivariateFinite(probs, augment = true, pool=missing)
    supp = Distributions.support(u)
    modes = mode.(u)
    @test modes isa CategoricalArray
    expected = [ifelse(p > 0.5, supp[2], supp[1]) for p in probs]
    @test all(modes .== expected)

    # multiclass
    rng = StableRNG(554)
    P   = rand(rng, n, c)
    P ./= sum(P, dims=2)
    u   = UnivariateFinite(P, pool=missing)
    expected = mode.([u...])
    @test all(mode.(u) .== expected)
end

@testset "cat for UnivariateFiniteArray" begin

    # ordered:
    v = categorical([:no, :yes, :maybe, :unseen])
    u1 = UnivariateFinite([v[1], v[2]], rand(5), augment=true)
    u2 = UnivariateFinite([v[3], v[2]], rand(6), augment=true)
    us = (u1, u2)
    u = cat(us..., dims=1)
    @test length(u) == length(u1) + length(u2)
    @test classes(u) == classes(u1)
    supp = Distributions.support(u)
    @test Set(supp) == Set([:no, :yes, :maybe])
    s1 = Distributions.support(u1)
    s2 = Distributions.support(u2)
    @test pdf(u, s1)[1:length(u1),:] == pdf(u1, s1)
    @test pdf(u, s2)[length(u1)+1:length(u1)+length(u2),:] ==
        pdf(u2, s2)
    @test pdf.(u, v[1])[length(u1)+1:length(u1)+length(u2)] ==
        zeros(length(u2))
    @test pdf.(u, v[3])[1:length(u1)] == zeros(length(u1))
    @test pdf.(u, v[4]) == zeros(length(u))

    # unordered:
    v = categorical([:no, :yes, :maybe, :unseen], ordered=true)
    u1 = UnivariateFinite([v[1], v[2]], rand(5), augment=true)
    u2 = UnivariateFinite([v[3], v[2]], rand(6), augment=true)
    us = (u1, u2)
    u = cat(us..., dims=1)
    @test length(u) == length(u1) + length(u2)
    @test classes(u) == classes(u1)
    supp = Distributions.support(u)
    @test Set(supp) == Set([:no, :yes, :maybe])
    s1 = Distributions.support(u1)
    s2 = Distributions.support(u2)
    @test pdf(u, s1)[1:length(u1),:] == pdf(u1, s1)
    @test pdf(u, s2)[length(u1)+1:length(u1)+length(u2),:] ==
        pdf(u2, s2)
    @test pdf.(u, v[1])[length(u1)+1:length(u1)+length(u2)] ==
        zeros(length(u2))
    @test pdf.(u, v[3])[1:length(u1)] == zeros(length(u1))
    @test pdf.(u, v[4]) == zeros(length(u))

    @test pdf([vcat(u1, u2)...], supp) ≈
        pdf(vcat([u1...], [u2...]), supp)
    h = hcat(u1, u1)
    h_nowrap = hcat([u1...], [u1...])
    @test size(h) == size(h_nowrap)
    # TODO: why does identity.(h) not work?
    @test h[3,2] ≈ h_nowrap[3,2]

    # errors
    v1 = categorical(1:2, ordered=true)
    v2 = categorical(v1, ordered=true)
    levels!(v2, levels(v2) |> reverse )
    probs = rand(3)
    u1 = UnivariateFinite(v1, probs, augment=true)
    u2 = UnivariateFinite(v2, probs, augment=true)
    @test_throws DomainError vcat(u1, u2)

    v1 = categorical(1:2)
    v2 = categorical(2:3)
    u1 = UnivariateFinite(v1, probs, augment=true)
    u2 = UnivariateFinite(v2, probs, augment=true)
    @test_throws DomainError vcat(u1, u2)

end

end

true


text:
γVals = [0.2,0.3,0.4,0.6,0.8]
nVals = [collect(3:20),collect(3:15),collect(3:10),collect(3:5),collect(3:5)]
cFracVals = collect(0.02:0.02:0.5)


text:
################################################################################
#
# Testcases for LinearTransformations.jl are added here.
#
################################################################################

using FluxHELayers.LinearTransformations

@testset "naive_matrix_multiply" begin
    SEALParams = setup_ckks(2^13, [40, 40, 40, 40, 40])
    initial_scale = 2.0^40

    vector = encrypt_array([1.0, 5.0, 9.0, 13.0], initial_scale, SEALParams)

    matrix = [1.0 2.0 3.0 4.0;
              5.0 6.0 7.0 8.0;
              9.0 10.0 11.0 12.0;
              13.0 14.0 15.0 16.0]

    output = naive_matrix_multiply(matrix,
                                   vector,
                                   SEALParams.evaluator,
                                   SEALParams.encoder,
                                   SEALParams.relinearization_keys)

    output = decrypt_array(output, SEALParams)

    @test output[1:4] ≈ [90, 202, 314, 426] atol=0.1
end

@testset "square_matrix_multiply" begin
    SEALParams = setup_ckks(2^13, [40, 40, 40, 40, 40])
    initial_scale = 2.0^40

    vector = encrypt_vector(repeat([1.0, 5.0, 9.0, 13.0], 1024),
                            initial_scale,
                            SEALParams)

    matrix = [1.0 2.0 3.0 4.0;
              5.0 6.0 7.0 8.0;
              9.0 10.0 11.0 12.0;
              13.0 14.0 15.0 16.0]

    output = square_matrix_multiply(matrix,
                                    vector,
                                    SEALParams.evaluator,
                                    SEALParams.encoder,
                                    SEALParams.galois_keys,
                                    SEALParams.relinearization_keys)

    output = decrypt_vector(output, SEALParams)

    @test output[1:4] ≈ [90, 202, 314, 426] atol=0.1
end

@testset "hybrid_matrix_multiply" begin
    SEALParams = setup_ckks(2^13, [40, 40, 40, 40, 40])
    initial_scale = 2.0^40

    vector = encrypt_vector(repeat([1.0, 5.0, 9.0, 13.0], 1024),
                            initial_scale,
                            SEALParams)

    matrix = [1.0 2.0 3.0 4.0;
              5.0 6.0 7.0 8.0]

    output = hybrid_matrix_multiply(matrix,
                                    vector,
                                    SEALParams.evaluator,
                                    SEALParams.encoder,
                                    SEALParams.galois_keys,
                                    SEALParams.relinearization_keys)

    output = decrypt_vector(output, SEALParams)

    @test output[1:2] ≈ [90, 202] atol=0.1
end


text:
using KitBase, FluxReconstruction, OrdinaryDiffEq, Langevin, LinearAlgebra, Plots
using ProgressMeter: @showprogress

function FR.positive_limiter(u::AbstractMatrix{T}, γ, weights) where {T<:AbstractFloat}
    # mean values
    u_mean = [
        sum(u[1, :] .* weights),
        sum(u[2, :] .* weights),
        sum(u[3, :] .* weights),
    ]
    t_mean = 1.0 / conserve_prim(u_mean, γ)[end]
    p_mean = 0.5 * u_mean[1] * t_mean

    # density corrector
    ϵ = min(1e-13, u_mean[1], p_mean)
    ρ_min = minimum(u[:, 1]) # density minumum can emerge at both solution and flux points
    t1 = min((u_mean[1] - ϵ) / (u_mean[1] - ρ_min + 1e-8), 1.0)
    @assert 0 < t1 <= 1 "incorrect range of limiter parameter t"

    for i in axes(u, 2)
        u[1, i] = t1 * (u[1, i] - u_mean[1]) + u_mean[1]
    end

    # energy corrector
    tj = Float64[]
    for i in axes(u, 2)
        prim = conserve_prim(u[:, i], γ)

        if prim[3] < ϵ
            prob = NonlinearProblem{false}(tj_equation, 1.0, (u[:, i], u_mean, γ, ϵ))
            sol = solve(prob, NewtonRaphson(), tol = 1e-6)
            push!(tj, sol.u)
        end
    end

    if length(tj) > 0
        t2 = minimum(tj)
        for j in axes(u, 2), i in axes(u, 1)
            u[i, j] = t2 * (u[i, j] - u_mean[i]) + u_mean[i]
        end
    end

    return nothing
end

function FR.positive_limiter(u::AbstractArray{T,3}, γ, wp, wq, ll, lr) where {T<:AbstractFloat}
    # tensorized quadrature weights
    weights = zeros(length(wp), length(wq))
    for i in axes(weights, 1), j in axes(weights, 2)
        weights[i, j] = wp[i] * wq[j]
    end

    # mean values
    u_mean = [
        sum(u[:, 1, :] .* weights),
        sum(u[:, 2, :] .* weights),
        sum(u[:, 3, :] .* weights),
    ]
    t_mean = 1.0 / conserve_prim(u_mean, γ)[end]

    # boundary variables
    ρb = zeros(2, length(wq))
    mb = zero(ρb)
    eb = zero(ρb)
    for j in axes(ρb, 2)
        ρb[:, j] .= [dot(u[:, 1, j], ll), dot(u[:, 1, j], lr)]
        mb[:, j] .= [dot(u[:, 2, j], ll), dot(u[:, 2, j], lr)]
        eb[:, j] .= [dot(u[:, 3, j], ll), dot(u[:, 3, j], lr)]
    end

    # density corrector
    ϵ = min(1e-13, u_mean[1], t_mean)
    ϵ < 0 && @show ϵ
    ρ_min = min(minimum(ρb), minimum(u[:, 1, :])) # density minumum can emerge at both solution and flux points
    t1 = min((u_mean[1] - ϵ) / (u_mean[1] - ρ_min + 1e-8), 1.0)
    @assert 0 <= t1 <= 1 "incorrect range of limiter parameter t"
    t1 < 1 && @show t1

    for i in axes(u, 1), j in axes(u, 3)
        u[i, 1, j] = t1 * (u[i, 1, j] - u_mean[1]) + u_mean[1]
    end

    # energy corrector
    tj = Float64[]
    for i = 1:2
        prim = conserve_prim([ρb[i], mb[i], eb[i]], γ)

        if 1 / prim[end] < ϵ
            prob = NonlinearProblem{false}(FR.tj_equation, 1.0, ([ρb[i], mb[i], eb[i]], u_mean, γ, ϵ))
            sol = solve(prob, NewtonRaphson(), tol = 1e-8)
            push!(tj, sol.u)
        end
    end
    for i in axes(u, 1), j in axes(u, 3)
        prim = conserve_prim(u[i, :, j], γ)

        if 1 / prim[end] < ϵ
            prob = NonlinearProblem{false}(FR.tj_equation, 1.0, (u[i, :, j], u_mean, γ, ϵ))
            sol = solve(prob, NewtonRaphson(), tol = 1e-8)
            push!(tj, sol.u)
        end
    end

    if length(tj) > 0
        t2 = minimum(tj)
        t2 < 1 && @show t2
        @assert 0 <= t2 <= 1 "incorrect range of limiter parameter t"
        for k in axes(u, 3), j in axes(u, 2), i in axes(u, 1)
            u[i, j, k] = t2 * (u[i, j, k] - u_mean[j]) + u_mean[j]
        end
    end

    return nothing
end

begin
    x0 = 0
    x1 = 1
    ncell = 100
    nface = ncell + 1
    dx = (x1 - x0) / ncell
    deg = 2 # polynomial degree
    nsp = deg + 1
    γ = 5 / 3
    cfl = 0.02
    dt = cfl * dx / (2.0)
    t = 0.0
    tspan = (0.0, 0.15)

    uqMethod = "galerkin"
    nr = 9
    nRec = 18
    opType = "uniform"
    parameter1 = 0.95
    parameter2 = 1.05
end
ps = FRPSpace1D(x0, x1, ncell, deg)
uq = UQ1D(nr, nRec, parameter1, parameter2, opType, uqMethod)
V = vandermonde_matrix(ps.deg,ps.xpl)
VInv = inv(Array(ps.V))
l2 = [uq.t2Product[j-1, j-1] for j = 1:uq.nm+1]

cd(@__DIR__)
include("rhs.jl")
include("../filter.jl")

begin
    isRandomLocation = true
    isPrefilter = true

    u = zeros(ncell, nsp, 3, uq.nm+1)
    if isRandomLocation
        # stochastic location
        for i = 1:ncell, j = 1:nsp
            prim = zeros(3, uq.nq)

            for k = 1:uq.nq
                if ps.x[i] <= 0.5 + 0.05 * uq.op.quad.nodes[k]
                    prim[:, k] .= [1.0, 0.0, 0.5]
                else
                    prim[:, k] .= [0.125, 0.0, 0.625]
                end
            end

            prim_chaos = zeros(3, uq.nm+1)
            for k = 1:3
                prim_chaos[k, :] .= ran_chaos(prim[k, :], uq)
            end

            u[i, j, :, :] .= uq_prim_conserve(prim_chaos, γ, uq)
        end
    else
        # stochastic density
        for i = 1:ncell, j = 1:nsp
            prim = zeros(3, uq.nm+1)
            if ps.x[i] <= 0.5
                prim[1, :] .= uq.pce
                #prim[1, 1] = 1.0
                prim[2, 1] = 0.0
                prim[3, 1] = 0.5
            else
                prim[:, 1] .= [0.125, 0.0, 0.625]
            end

            u[i, j, :, :] .= uq_prim_conserve(prim, γ, uq)
        end
    end

    if isPrefilter
        # pre-filtering
        for iter = 1:2
            for j = 1:size(u, 1)
                for s = 1:size(u, 3)
                    uModal = VInv * u[j, :, s, :]
                    #FR.modal_filter!(uModal, 15e-2, 10e-5; filter = :l2opt)
                    #FR.modal_filter!(uModal, 10, 10; filter = :exp)
                    FR.modal_filter!(uModal; filter = :lasso)
                    FR.modal_filter!(uModal, 5e-2, 5e-2; filter = :l2)
                    u[j, :, s, :] .= V * uModal
                end
            end
        end
    end
end

p = (ps.J, ps.ll, ps.lr, ps.dl, ps.dhl, ps.dhr, γ, uq)
prob = ODEProblem(dudt!, u, tspan, p)
nt = tspan[2] ÷ dt |> Int
itg = init(prob, Midpoint(), saveat = tspan[2], adaptive = false, dt = dt)

@showprogress for iter = 1:nt
    for i = 1:size(itg.u, 1)
        begin
            uNodal = zeros(nsp, 3, uq.nq)
            for idx = 1:nsp, jdx = 1:3
                uNodal[idx, jdx, :] .= chaos_ran(itg.u[i, idx, jdx, :], uq)
            end
            #positive_limiter(uNodal, γ, ps.wp/2, uq.op.quad.weights, ps.ll, ps.lr)
            for idx = 1:nsp
                tmp = @view uNodal[idx, :, :]
                positive_limiter(tmp, γ, uq.op.quad.weights)
            end
            for idx = 1:nsp, jdx = 1:3
                itg.u[i, idx, jdx, :] .= ran_chaos(uNodal[idx, jdx, :], uq)
            end
        end

        ũ = VInv * itg.u[i, :, 1, :]
        su = maximum([ũ[end, j]^2 / sum(ũ[:, j].^2) for j = 1:uq.nm+1])
        sv = maximum([ũ[j, end]^2 * uq.t2Product[uq.nm, uq.nm] / sum(ũ[j, :].^2 .* l2) for j = 1:nsp])
        isShock = max(shock_detector(log10(su), ps.deg), shock_detector(log10(sv), ps.deg))
        if isShock
            λ1 = dt * (su) * 5.0
            λ2 = dt * (sv) * 5.0

            for s = 1:size(itg.u, 3)
                û = VInv * itg.u[i, :, s, :]
                #FR.modal_filter!(û, λ1, λ2; filter = :l2)
                #FR.modal_filter!(û, λ1, λ2; filter = :l2opt)
                #FR.modal_filter!(û, 1e-2, 1e-2; filter = :l2)
                FR.modal_filter!(û; filter = :lasso)

                itg.u[i, :, s, :] .= ps.V * û
            end
        end
    end

    step!(itg)
end

begin
    x = zeros(ncell * nsp)
    w = zeros(ncell * nsp, 3, uq.nm+1)
    for i = 1:ncell
        idx0 = (i - 1) * nsp

        for j = 1:nsp
            idx = idx0 + j
            x[idx] = ps.xpg[i, j]

            w[idx, :, :] .= itg.u[i, j, :, :]
        end
    end

    sol = zeros(ncell*nsp, 3, 2)
    for i in axes(sol, 1)
        p1 = zeros(3, uq.nm+1)
        p1 = uq_conserve_prim(w[i, :, :], γ, uq)
        p1[end, :] .= lambda_tchaos(p1[end, :], 1.0, uq)

        for k = 1:3
            sol[i, k, 1] = mean(p1[k, :], uq.op)
            sol[i, k, 2] = std(p1[k, :], uq.op)
        end
    end

    pic1 = plot(x, sol[:, 1, 1], label="ρ", xlabel="x", ylabel="mean")
    plot!(pic1, x, sol[:, 2, 1], label="U")
    plot!(pic1, x, sol[:, 3, 1], label="T")
    pic2 = plot(x, sol[:, 1, 2], label="ρ", xlabel="x", ylabel="std")
    plot!(pic2, x, sol[:, 2, 2], label="U")
    plot!(pic2, x, sol[:, 3, 2], label="T")
    plot(pic1, pic2)
end

plot(x, sol[:, 1, 1], label="ρ", xlabel="x", ylabel="mean")
plot(x, sol[:, 1, 2], label="ρ", xlabel="x", ylabel="std")

max_stars_repo_name: CSMMLab/FRSG

max_issues_repo_name: CSMMLab/FRSG

max_forks_repo_name: CSMMLab/FRSG


text:
# Autogenerated wrapper script for SPIRV_Cross_jll for x86_64-linux-musl-cxx03
export libspirv_cross, spirv_cross

JLLWrappers.@generate_wrapper_header("SPIRV_Cross")
JLLWrappers.@declare_library_product(libspirv_cross, "libspirv-cross-c-shared.so.0")
JLLWrappers.@declare_executable_product(spirv_cross)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libspirv_cross,
        "lib/libspirv-cross-c-shared.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        spirv_cross,
        "bin/spirv-cross",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function kalman_smoother!(
    smoother_output::SmootherOutput,
    system::LinearUnivariateTimeInvariant,
    filter_output::FilterOutput,
)
    multivariate_system = to_multivariate_time_variant(system)
    return kalman_smoother!(smoother_output, multivariate_system, filter_output)
end

function kalman_smoother!(
    smoother_output::SmootherOutput,
    system::LinearMultivariateTimeVariant{Fl},
    filter_output::FilterOutput,
) where Fl
    # Load dimensions data
    n = size(system.y, 1)
    m = size(system.T[1], 1)

    # Load filter data
    v = filter_output.v
    F = filter_output.F
    a = filter_output.a
    P = filter_output.P

    # Smoothed state and its covariance
    L = Vector{Matrix{Fl}}(undef, n)
    r = Vector{Vector{Fl}}(undef, n)
    N = Vector{Matrix{Fl}}(undef, n)

    # Initialization
    N[end] = zeros(Fl, m, m)
    r[end] = zeros(Fl, m)

    @inbounds for t in n:-1:2
        if any(y_ -> isnan(y_), system.y[t, :])
            r[t - 1] = system.T'[t] * r[t]
            N[t - 1] = system.T'[t] * N[t] * system.T[t]
        else
            Z_transp_invF = system.Z[t]' * inv(F[t])
            K = system.T[t] * P[t] * Z_transp_invF
            L[t] = system.T[t] - K * system.Z[t]
            r[t - 1] = Z_transp_invF * v[t] + L[t]' * r[t]
            N[t - 1] = Z_transp_invF * system.Z[t] + L[t]' * N[t] * L[t]
        end
        smoother_output.alpha[t] = a[t] + P[t] * r[t - 1]
        smoother_output.V[t] = P[t] - P[t] * N[t - 1] * P[t]
    end

    # Last iteration
    if any(y_ -> isnan(y_), system.y[1, :])
        r0 = system.T'[1] * r[1]
        N0 = system.T'[1] * N[1] * system.T[1]
    else
        Z_transp_invF = system.Z[1]' * inv(F[1])
        K = system.T[1] * P[1] * Z_transp_invF
        L[1] = system.T[1] - K * system.Z[1]
        r0 = Z_transp_invF * v[1] + L[1]' * r[1]
        N0 = Z_transp_invF * system.Z[1] + L[1]' * N[1] * L[1]
    end
    smoother_output.alpha[1] = a[1] + P[1] * r0
    smoother_output.V[1] = P[1] - P[1] * N0 * P[1]
    return smoother_output
end


text:
export
    Env,
    step,
    reset,
    observation_space_spec,
    action_space_spec,
    render,
    obs_names,
    reward_names,
    vectorized,
    num_envs

abstract type Env end
Base.step(env::Env, action::Int) = error("Not implemented")
Base.step(env::Env, action::Float64) = error("Not implemented")
Base.step(env::Env, action::Array{Float64}) = error("Not implemented")
Base.reset(env::Env) = error("Not implemented")
Base.reset(env::Env; kwargs...) = error("Not implemented")
Base.reset(env::Env, dones::Union{Vector{Bool}, Void}; kwargs...) = reset(env; kwargs...)
observation_space_spec(env::Env) = error("Not implemented")
action_space_spec(env::Env) = error("Not implemented")
render(env::Env) = error("Not implemented")
obs_names(env::Env) = error("Not implemented")
vectorized(env::Env) = false
num_envs(env::Env) = 1


text:
using Base.Test
using ThreePhotons

#Postprocess runs

params = Dict("L"=>L, "K3_range" =>1:K, "N"=>N, "LMAX"=>LMAX, "K2_range"=>1:KMAX, "rmax"=>float(KMAX), "qmax"=>float(KMAX))
state = Dict{Any,Any}("intensity"=>deleteTerms(intensity, K, L))
postprocess_run(params, state, "../data/structures/crambin.pdb", false)


text:
using Documenter
using Peaks

DocMeta.setdocmeta!(Peaks, :DocTestSetup, :(using Peaks); recursive=true)
makedocs(
    sitename = "Peaks",
    modules = [Peaks],
    format = Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://halleysfifthinc.github.io/Peaks.jl",
    ),
    authors="Allen Hill <allenofthehills@gmail.com> and contributors",
    repo="https://github.com/halleysfifthinc/Peaks.jl/blob/{commit}{path}#L{line}",
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/halleysfifthinc/Peaks.jl",
    push_preview=true,
)

max_stars_repo_name: halleysfifthinc/Peaks

max_issues_repo_name: halleysfifthinc/Peaks

max_forks_repo_name: halleysfifthinc/Peaks


text:
const ColumnIndex = Union{Symbol,Integer}

"""
An RLEDataFrame extends DataFrame and contains a colection of like-length and like-type
    RLEVectors. In a way, this creates a type like an RLE matrix. But, we deliberately
    avoid the complexity of matrix operations, such as factorization. It is expected
    that most operations will be column-wise. Based on RleDataFrame from Bioconductor's
    `genoset` package (also by Peter Haverty).
### Examples
```julia
x = RLEDataFrame( [RLEVector([1, 1, 2]), RLEVector([2, 2, 2])], [:a, :b])
y = RLEDataFrame( [RLEVector([5])],[:a] )
z = RLEDataFrame( a=RLEVector([5,2,2]), b=RLEVector([4,4,4])
```
"""
mutable struct RLEDataFrame{T1,T2<:Integer}
    columns::Vector{RLEVector{T1,T2}}
    colindex::NamedTuple
    function RLEDataFrame{T1,T2}(columns::Vector{RLEVector{T1,T2}}, colnames::Vector{Symbol}) where {T1,T2}
        ncol = length(columns)
        nval = length(colnames)
        if ncol != nval
            throw(ArgumentError("Incoming columns collection and column names must have the same length."))
        end
        if ncol > 0
            lens = map(length, columns)
            (min, max) = Base.extrema(lens)
            if min != max
                throw(ArgumentError("All incoming columns must be of equal length."))
            end
        end
        c = Tuple( Symbol(x) for x in colnames )
        colindex = NamedTuple{c}(1:ncol)
        new(columns, colindex)
    end
end

==(x::RLEDataFrame, y::RLEDataFrame) = x.colindex == y.colindex && x.columns == y.columns

nrow(x::RLEDataFrame) = length(x.columns[1])
ncol(x::RLEDataFrame) = length(x.columns)
Base.length(x::RLEDataFrame) = length(x.columns)
index(x::RLEDataFrame) = x.colindex
columns(x::RLEDataFrame) = x.columns
Base.names(x::RLEDataFrame) = collect(keys(x.colindex))
Base.size(x::RLEDataFrame) = (nrow(x), ncol(x))

function Base.show(io::IO, x::RLEDataFrame)
    t = typeof(x)
    show(io, t)
    println()
    for (c,v) in zip(names(x),columns(x))
       println(io,"Column: $c")
        println(io,v)
    end
end

function RLEDataFrame(cols, colnames)
    f = cols[1]
    RLEDataFrame{eltype(f.runvalues),eltype(f.runends)}(cols, colnames)
end

function RLEDataFrame(; kwargs...)
    cnames  = [k for (k,v) in kwargs]
    cvalues = [v for (k,v) in kwargs]
    RLEDataFrame(cvalues,cnames)
end

Base.copy(x::RLEDataFrame) = RLEDataFrame(copy(x.columns), names(x))

### Get/set
## Just columns
Base.getindex(x::RLEDataFrame,j::Colon) = copy(x)
Base.getindex(x::RLEDataFrame,j::ColumnIndex) = columns(x)[index(x)[j]]
function Base.getindex(x::RLEDataFrame,j::AbstractArray)
    ind = index(x)
    inds = [ ind[x] for x in j ]
    RLEDataFrame( columns(x)[inds], names(x)[inds] )
end

function Base.setindex!(x::RLEDataFrame, value::AbstractVector, j::Integer)
    if length(value) != nrow(x)
        throw(ArgumentError("Length of incoming value must match existing columns."))
    end
    if j <= length(x)
        columns(x)[j] = value
    else
        throw(BoundsError())
    end
    x
end

function Base.setindex!(x::RLEDataFrame, value::AbstractVector, j::Symbol)
    if length(value) != nrow(x)
        throw(ArgumentError("Length of incoming value must match existing columns."))
    end
    if j in names(x)
        columns(x)[index(x)[j]] = value
    else
        x.colindex = merge(index(x), NamedTuple{(j,)}( Tuple(length(x) + 1) ))
        x.columns = push!(x.columns,value)
    end
    x
end

## with rows
function Base.getindex(x::RLEDataFrame, i, j)
    ind = index(x)
    j_inds = [ ind[x] for x in j ]
    cols = [ x.columns[j_ind][i] for j_ind in j_inds ]
    RLEDataFrame( cols, names(x)[j_inds] )
end
Base.getindex(x::RLEDataFrame, i::Integer, j::ColumnIndex) = x[j][i]
Base.getindex(x::RLEDataFrame, i::Integer, j) = x[ [i], j ]
Base.getindex(x::RLEDataFrame, i, j::ColumnIndex) = x[j][i]

function Base.setindex!(x::RLEDataFrame, value, i, j)
    ind = index(x)
    j_inds = [ ind[x] for x in j ]
    for j_ind in j_inds
        x.columns[j_ind][i] = value
    end
    x
end
function Base.setindex!(x::RLEDataFrame, value, i, j::ColumnIndex)
    x[j][i] = value
end

## Conversion
Base.convert(Matrix, x::RLEDataFrame) = hcat(map(collect,x.columns)...)

### Familiar operations over rows or columns from R

# Probably these are all a job for mapslice or slicedim. I need to RTM.
rowmap(x::Matrix,f::Function) = [ f( @view x[i,:] ) for i in 1:size(x)[1] ]
colmap(x::Matrix,f::Function) = [ f( @view x[:,j] ) for j in 1:size(x)[2] ]
rowMeans(x) = rowmap(x,mean)
rowMedians(x) = rowmap(x,median)
rowSums(x) = rowmap(x,sum)
colMeans(x) = colmap(x,mean)
colMedians(x) = colmap(x,median)
colSums(x) = colmap(x,sum)

#rowMeans(x::RLEDataFrame) = rowSum(x) ./ ncol(x)
colSums(x::RLEDataFrame) = map(sum, columns(x))
colMeans(x::RLEDataFrame) = colSums(x) ./ nrow(x)
colMedians(x::RLEDataFrame) = map(median, columns(x))


text:
@reexport module Threshold

using SpikingNNFunctions.Threshold: poisson
using Adapt
using Random

import ..SpikingNN: excite!, evaluate!, reset!, isactive

export AbstractThreshold, isactive

abstract type AbstractThreshold end

"""
    Ideal(vth::Real)

An ideal threshold spikes when `v > vth`.
"""
struct Ideal{T<:Real} <: AbstractThreshold
    vth::T
end

isactive(threshold::Ideal, t::Integer; dt::Real = 1.0) = false

"""
    evaluate!(threshold::Ideal, t::Real, v::Real; dt::Real = 1.0)
    (threshold::Ideal)(t::Real, v::Real; dt::Real = 1.0)
    evaluate!(thresholds::AbstractArray{<:Ideal}, t::Integer, v; dt::Real = 1.0)

Return `t` when `v > threshold.vth`.
"""
evaluate!(threshold::Ideal, t::Real, v::Real; dt::Real = 1.0) = (v >= threshold.vth) ? t : zero(t)
(threshold::Ideal)(t::Real, v::Real; dt::Real = 1.0) = evaluate!(threshold, t, v; dt = dt)
evaluate!(thresholds::T, t::Integer, v; dt::Real = 1.0) where T<:AbstractArray{<:Ideal} =
    Int.((v .>= thresholds.vth) .* adapt(typeof(v), fill(t, size(v))))

"""
    Poisson(ρ₀::Real, Θ::Real, Δᵤ::Real, rng::AbstractRNG)

Choose to output a spike based on a inhomogenous Poisson process given by

``X < \\mathrm{d}t \\: \\rho_0 \\exp\\left(\\frac{v - \\Theta}{\\Delta_u}\\right)``

where ``X \\sim \\mathrm{Unif}([0, 1])``.
Accordingly, `dt` must be set correctly so that the neuron does not always spike.

Fields:
- `ρ₀::Real`: baseline firing rate at threshold
- `Θ::Real`: firing threshold
- `Δᵤ::Real`: voltage resolution
- `rng`: random number generation
"""
struct Poisson{T<:Real,  RT<:AbstractRNG} <: AbstractThreshold
    ρ₀::T
    Θ::T
    Δᵤ::T
    rng::RT
end

function Poisson{T}(ρ₀::Real, Θ::Real, Δᵤ::Real; rng::RT = Random.GLOBAL_RNG) where {T<:Real, RT}
    Poisson{T,RT}(ρ₀, Θ, Δᵤ, rng)
end
Poisson(ρ₀::Real, Θ::Real, Δᵤ::Real; kwargs...) = Poisson{Real}(ρ₀, Θ, Δᵤ; kwargs...)

isactive(threshold::Poisson, t::Integer; dt::Real = 1.0) = true

"""
    evaluate!(threshold::Poisson, t::Integer, v::Real; dt::Real = 1.0)
    (::Poisson)(t::Integer, v::Real; dt::Real = 1.0)
    evaluate!(thresholds::AbstractArray{<:Poisson}, t::Integer, v; dt::Real = 1.0)

Evaluate Poisson threshold function. See [`Threshold.Poisson`](@ref).
"""
evaluate!(threshold::Poisson, t::Integer, v::Real; dt::Real = 1.0) =
    poisson(threshold.ρ₀, threshold.Θ, threshold.Δᵤ, v; dt = dt, rng = threshold.rng) ? t : zero(t)
(threshold::Poisson)(t::Integer, v::Real; dt::Real = 1.0) = evaluate!(threshold, t, v; dt = dt)
evaluate!(thresholds::T, t::Integer, v; dt::Real = 1.0) where T<:AbstractArray{<:Poisson} =
    Int.(poisson(thresholds.ρ₀, thresholds.Θ, thresholds.Δᵤ, v; dt = dt, rng = threshold.rng) .* adapt(typeof(v), fill(t, size(v))))

end


text:
"""
    get_preconditioner(model::AbstractModel{T,N},op::LinearMap{T}) where {T, N}

 Get precondtioner. Details are stored in a struct of type Preconditioner that can be passed to a solver.

"""
function get_preconditioner(model::AbstractModel{T,N},op::LinearMap{T}) where {T, N}

    @unpack gu,gv,wu,wv=model.fields
    @unpack params,solver_params=model

    m,n=size(op)
    @assert m == n == gu.n + gv.n

    n = gu.n + gv.n
    n_coarse = wu.n[] + wv.n[]

    restrict_fun(x) = restrictvec(model,x)
    prolong_fun(x) = prolongvec(model,x)

    restrict=LinearMap{T}(restrict_fun,n_coarse,n;issymmetric=false,ismutating=false,ishermitian=false,isposdef=false)
    prolong=LinearMap{T}(prolong_fun,n,n_coarse;issymmetric=false,ismutating=false,ishermitian=false,isposdef=false)

    op_diag=get_op_diag(model,op)

    #Four colour Jacobi preconditioner. Red-Black checkerboard Jacobi for each velocity component.
    sweep=[[1 .+ mod(i-j,2) for i=1:gu.nxu, j=1:gu.nyu][gu.mask];[3 .+ mod(i-j,2) for i=1:gv.nxv, j=1:gv.nyv][gv.mask]]
    sweep_order=[1,3,2,4]

    p=Preconditioner{T,N}(op=op, restrict=restrict, prolong=prolong,sweep=sweep, sweep_order=sweep_order,
            op_diag=op_diag, nsmooth=solver_params.nsmooth, tol_coarse = solver_params.tol_coarse,
            maxiter_coarse = solver_params.maxiter_coarse, smoother_omega=solver_params.smoother_omega)

    return p
end

"""
    precondition!(x, p, b)

Apply wavelet-based multigrid preconditioner using information stored in p.
"""
function precondition!(x, p, b)
    @unpack op,op_diag,nsmooth,sweep,sweep_order,smoother_omega,restrict,
            prolong,op_coarse,correction_coarse,tol_coarse,maxiter_coarse = p

    n=size(op,1)

    # Multigrid smooth
    x .= gauss_seidel_smoother!(x, op, b; iters = nsmooth, op_diag=op_diag,
                                sweep=sweep, sweep_order=sweep_order, smoother_omega = smoother_omega)

    resid=b-op*x;

    # Multigrid restriction
    b_coarse=restrict*resid

    # Multigrid solve for correction
    cg!(correction_coarse, op_coarse, b_coarse; reltol = tol_coarse, maxiter = maxiter_coarse)

    # Multigrid prolongation
    x .= x .+ prolong*correction_coarse

    # Multigrid smooth
    x .= gauss_seidel_smoother!(x, op, b; iters = nsmooth, op_diag=op_diag,
                                sweep=sweep, sweep_order=reverse(sweep_order), smoother_omega = smoother_omega)

    return x
end

"""
    get_op_diag(wavi::AbstractModel,op::LinearMap)

 Get diagonal of operator for use in preconditioner.

"""
function get_op_diag(model::AbstractModel,op::LinearMap)
    @unpack gu,gv=model.fields
    @unpack params,solver_params=model
    m,n=size(op)
    @assert m == n == gu.n + gv.n
    op_diag=zeros(eltype(op),n)
    sm=solver_params.stencil_margin
    sweep=[[1+mod((i-1),sm)+sm*mod((j-1),sm) for i=1:gu.nxu, j=1:gu.nyu][gu.mask];
           [1+sm^2+mod((i-1),sm)+sm*mod((j-1),sm) for i=1:gv.nxv, j=1:gv.nyv][gv.mask] ]
    e=zeros(Bool,n)
    for i = unique(sweep)
        e[sweep .== i] .= true
        op_diag[e] .= (op*e)[e]
        e[sweep .== i] .= false
    end
    return op_diag
end

"""
    LinearAlgebra.ldiv!(x::AbstractVecOrMat{T}, p::AbstractPreconditioner{T,N}, b::AbstractVecOrMat{T}) where {T,N}

Overload LinearAlgebra.ldiv! function so that the bespoke preconditioner is deployed in calls to the
conjugate gradient method if p has type  <: AbstractPreconditioner.

"""
function ldiv!(x::AbstractVecOrMat{T}, p::AbstractPreconditioner{T,N}, b::AbstractVecOrMat{T}) where {T,N}
    precondition!(x, p, b)
end

"""
    gauss_seidel_smoother!(x, op, b;
                                iters=5,
                                op_diag=diag(sparse(op)),
                                sweep=(1:size(op,1)),
                                sweep_order=unique(sweep),
                                smoother_omega=1.0)
Apply smoother used in multigrid preconditioner.
"""
function gauss_seidel_smoother!(x, op, b;
                                iters=5,
                                op_diag=diag(sparse(op)),
                                sweep=(1:size(op,1)),
                                sweep_order=unique(sweep),
                                smoother_omega=1.0)
    resid=b-op*x
    for i = 1:iters
        for j = sweep_order
              idx = sweep .== j
              x[idx] .= x[idx] .+ smoother_omega .* resid[idx]./op_diag[idx]
              resid .= b .- op*x
        end
    end
    return x
end


text:
Base.@propagate_inbounds function _swap_rows!(A::AbstractMatrix, i, j)
    @boundscheck @assert 1 ≤ i ≤ size(A, 1)
    @boundscheck @assert 1 ≤ j ≤ size(A, 1)

    i == j && return A

    @inbounds for col_idx in 1:size(A, 2)
        A[i, col_idx], A[j, col_idx] = A[j, col_idx], A[i, col_idx]
    end
    return A
end

Base.@propagate_inbounds function _mul_row!(
    A::AbstractMatrix,
    val,
    row_idx,
    starting_at = 1,
)
    @boundscheck @assert 1 ≤ starting_at ≤ size(A, 2)

    @inbounds for col_idx in starting_at:size(A, 2)
        A[row_idx, col_idx] *= val
    end
    return A
end

# Version over exact field:

Base.@propagate_inbounds function _find_pivot(
    A::AbstractMatrix,
    col_idx,
    starting_at = 1,
)
    k = findnext(!iszero, @view(A[:, col_idx]), starting_at)
    isnothing(k) && return false, col_idx
    return true, k
end

Base.@propagate_inbounds function _reduce_column_by_pivot!(
    A::AbstractMatrix,
    row_idx,
    col_idx,
    starting_at = 1,
)
    @boundscheck checkbounds(A, row_idx, col_idx)
    @boundscheck 1 ≤ starting_at ≤ size(A, 2)

    for ridx in 1:size(A, 1)
        ridx == row_idx && continue
        v = A[ridx, col_idx]
        iszero(v) && continue
        @inbounds for cidx in starting_at:size(A, 2)
            A[ridx, cidx] -= v * A[row_idx, cidx]
        end
    end
    return A
end

# version over AbstractFloat

Base.@propagate_inbounds function __findmax(
    f,
    A::AbstractArray{T};
    atol=eps(real(eltype(A)))*max(size(A)...)
) where {T<:Union{AbstractFloat,Complex{<:AbstractFloat}}}
    @assert !isempty(A)

    midx = first(eachindex(A))
    mval = f(A[midx])

    @inbounds for idx in eachindex(A)
        iszero(A[idx]) && continue
        v = f(A[idx])
        if v > mval
            midx, mval = idx, v
        end
        if v < atol
            A[idx] = zero(T)
        end
    end

    return midx, mval
end

Base.@propagate_inbounds function _find_pivot(
    A::AbstractMatrix{T},
    col_idx,
    starting_at = 1;
    atol = eps(real(eltype(A))) * size(A, 1),
) where {T<:Union{AbstractFloat,Complex{<:AbstractFloat}}}
    isempty(starting_at:size(A, 1)) && return false, starting_at
    @boundscheck @assert 1 ≤ starting_at
    # find the largest entry in the column below the last pivot
    @boundscheck @assert 1 ≤ col_idx ≤ size(A, 2)

    midx, mval = __findmax(abs, @view(A[starting_at:end, col_idx]), atol=atol)
    if mval < atol # the largest entry is below threshold so we zero everything in the column
        return false, starting_at
    end
    return true, oftype(starting_at, starting_at + midx - 1)
end

Base.@propagate_inbounds function _reduce_column_by_pivot!(
    A::AbstractMatrix{T},
    row_idx,
    col_idx,
    starting_at = 1,
    atol = eps(real(eltype(A))) * size(A, 1),
) where {T<:Union{AbstractFloat,Complex{<:AbstractFloat}}}
    @boundscheck checkbounds(A, row_idx, col_idx)
    @boundscheck 1 ≤ starting_at ≤ size(A, 2)

    for ridx in 1:size(A, 1)
        ridx == row_idx && continue
        @inbounds v = A[ridx, col_idx]
        if abs(v) < atol
            @inbounds A[ridx, col_idx] = zero(T)
            continue
        end
        for cidx in starting_at:size(A, 2)
            @inbounds A[ridx, cidx] -=
                cidx == col_idx ? A[ridx, cidx] : v * A[row_idx, cidx]
        end
    end
    return A
end

Base.@propagate_inbounds function row_echelon_form!(A::AbstractMatrix)
    pivots = Int[]
    row_idx = 0 # also: current_rank
    @inbounds for col_idx in 1:size(A, 2)
        found, j = _find_pivot(A, col_idx, row_idx + 1)
        found || continue
        row_idx += 1
        push!(pivots, col_idx)

        # swap the rows so that  A[row_idx, :] is the row with leading nz in
        # i-th column
        _swap_rows!(A, row_idx, j)

        w = inv(A[row_idx, col_idx])

        # multiply A[row_idx, :] by w
        _mul_row!(A, w, row_idx, col_idx)
        # A[current_rank, col_idx] is 1 now

        # zero the whole col_idx-th column above and below pivot:
        # to the left of col_idx everything is zero
        _reduce_column_by_pivot!(A, row_idx, col_idx, col_idx)
    end
    return A, pivots
end

row_echelon_form(A::AbstractMatrix) = row_echelon_form!(deepcopy(A))


text:
a_wgs84=6378137.
f_wgs84 = 1/298.257223563
b_wgs84 = a_wgs84*(1-f_wgs84)
e_wgs84 = sqrt((a_wgs84^2-b_wgs84^2)/a_wgs84^2)

#export ECEFtoGeodetic, GeodetictoECEF
export GeodetictoECEF

function GeodetictoECEF2(lat::Number, lon::Number, h::Number)
    # Auxiliary variables.
    sin_lat, cos_lat = sincos(lat)
    sin_lon, cos_lon = sincos(lon)

    # Radius of curvature [m].
    N = a_wgs84/sqrt(1 - e_wgs84^2*sin_lat^2 )

    # Compute the position in ECEF frame.
    SVector{3}((                     N + h)*cos_lat*cos_lon,
               (                     N + h)*cos_lat*sin_lon,
               ( (b_wgs84/a_wgs84)^2*N + h)*sin_lat)

    println(SVector)
end

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
using TensorOperations
using LoopVectorization

function RHF(molecule::Molecule, aoint::IntegralHelper, C::Array{Float64,2}, Λ, Alg::ConventionalRHF)
    @output "Computing integrals ..."
    t = @elapsed aoint["μ"]
    @output " done in {:>5.2f} s\n" t
    RHF(molecule,aoint,C,aoint["μ"],Λ)
end
function RHF(molecule::Molecule, aoint::IntegralHelper, C::Array{Float64,2}, Λ, Alg::DFRHF)
    @output "Computing integrals ..."
    t = @elapsed begin
        Amn = aoint["B"]
        J = aoint["J"]
        Qmn = Array{Float64}(undef,size(Amn))
        Fermi.contract!(Qmn,J,Amn,0.0,1.0,1.0,"Qmn","PQ","Pmn")
    end
    @output " done in {:>5.2f} s\n" t
    RHF(molecule,aoint,C,Qmn,Λ)
end
function RHF(molecule::Molecule, aoint::IntegralHelper, C::Array{Float64,2}, ERI::Array{Float64}, Λ::Array)
    Fermi.HartreeFock.print_header()

    #grab some options
    maxit = Fermi.CurrentOptions["scf_max_iter"]
    Etol  = 10.0^(-Fermi.CurrentOptions["e_conv"])
    Dtol  = Fermi.CurrentOptions["scf_max_rms"]
    do_diis = Fermi.CurrentOptions["diis"]
    oda = Fermi.CurrentOptions["oda"]
    oda_cutoff = Fermi.CurrentOptions["oda_cutoff"]
    oda_shutoff = Fermi.CurrentOptions["oda_shutoff"]

    #variables that will get updated iteration-to-iteration
    ite = 1
    E = 0.0
    ΔE = 1.0
    Drms = 1.0
    Drms = 1.0
    diis = false
    damp = 0.0
    converged = false

    #build a diis_manager, if needed
    do_diis ? DM = Fermi.DIIS.DIISManager{Float64,Float64}(size=Fermi.CurrentOptions["ndiis"]) : nothing
    do_diis ? diis_start = Fermi.CurrentOptions["diis_start"] : nothing

    #grab ndocc,nvir
    ndocc = try
        Int((molecule.Nα + molecule.Nβ)/2)
    catch InexactError
        throw(Fermi.InvalidFermiOption("Invalid number of electrons $(molecule.Nα + molecule.Nβ) for RHF method."))
    end
    nvir = size(C)[2] - ndocc
    nao = size(C)[1]

    @output " Number of AOs:                        {:5.0d}\n" nao
    @output " Number of Doubly Occupied Orbitals:   {:5.0d}\n" ndocc
    @output " Number of Virtual Spatial Orbitals:   {:5.0d}\n" nvir

    S = Array(Hermitian(aoint["S"]))
    J = aoint["J"]^-2
    T = aoint["T"]
    V = aoint["V"]
    dfsz = size(ERI,3)
    mn = zeros(size(S))
    Threads.@threads for m=1:nao
        for n=1:nao
            s = 0.0
            @avx for P=1:dfsz
                s += ERI[P,m,n]
            end
            mn[m,n] = s
        end
    end

    # Form the density matrix from occupied subset of guess coeffs
    Co = C[:, 1:ndocc]
    D = Fermi.contract(Co,Co,"um","vm")
    D_old = deepcopy(D)

    eps = zeros(Float64,ndocc+nvir)
    # Build the inital Fock Matrix and diagonalize
    F = zeros(Float64,nao,nao)
    build_fock!(F, T + V, D, ERI, Co, mn, J)
    F̃ = deepcopy(F)
    D̃ = deepcopy(D)
    @output " Guess Energy {:20.14f}\n" RHFEnergy(D,T+V,F)

    @output "\n Iter.   {:>15} {:>10} {:>10} {:>8} {:>8} {:>8}\n" "E[RHF]" "ΔE" "√|ΔD|²" "t" "DIIS" "damp"
    @output repeat("-",80)*"\n"
    t = @elapsed while ite ≤ maxit
        t_iter = @elapsed begin
            # Produce Ft
            if !oda || Drms < oda_cutoff
                Ft = Λ'*F*Λ
            else
                Ft = Λ'*F̃*Λ
            end

            # Get orbital energies and transformed coefficients
            eps,Ct = eigen(Hermitian(real.(Ft)))

            # Reverse transformation to get MO coefficients
            C = Λ*Ct
            #Ct = real.(Ct)

            # Produce new Density Matrix
            Co = C[:,1:ndocc]
            #D = Fermi.contract(Co,Co,"um","vm")
            Fermi.contract!(D,Co,Co,0.0,1.0,1.0,"uv","um","vm")

            # Build the Fock Matrix
            build_fock!(F, T + V, D, ERI, Co, mn, J)
            Eelec = RHFEnergy(D, T + V, F)

            # Compute Energy
            Enew = Eelec + molecule.Vnuc

            #branch for ODA vs DIIS convergence aids
            if oda && Drms > oda_cutoff && ite < oda_shutoff
                diis = false
                dD = D - D̃
                s = tr(F̃ * dD)
                c = tr((F - F̃) * (dD))
                if c <= -s/2
                    λ = 1.0
                else
                    λ = -s/(2*c)
                end
                F̃ .= (1-λ)*F̃ + λ*F
                D̃ .= (1-λ)*D̃ + λ*D
                damp = 1-λ
                do_diis ? err = transpose(Λ)*(F*D*aoint["S"] - aoint["S"]*D*F)*Λ : nothing
                do_diis ? push!(DM, F, err) : nothing
            elseif (!oda || ite > oda_shutoff || Drms < oda_cutoff) && do_diis
                damp = 0.0
                diis = true
                D̃ = D
                FP = transpose(Λ)*(F*D*aoint["S"] - aoint["S"]*D*F)*Λ
                do_diis ? err = transpose(Λ)*(F*D*aoint["S"] - aoint["S"]*D*F)*Λ : nothing
                do_diis ? push!(DM, F, err) : nothing
                #do_diis && ite > diis_start ? F,_ = Fermi.DIIS.extrapolate(DM) : nothing
                if do_diis && ite > diis_start
                    F = Fermi.DIIS.extrapolate(DM)
                end
            end
            FPrms = sqrt(sum(err.^2))

            # Compute the Density RMS
            ΔD = D - D_old
            Drms = sqrt(sum(ΔD.^2))

            # Compute Energy Change
            ΔE = Enew - E
            E = Enew
            D_old .= D
        end
        @output "    {:<3} {:>15.10f} {:>11.3e} {:>11.3e} {:>8.2f} {:>8}    {:5.4f}\n" ite E ΔE Drms t_iter diis damp
        ite += 1

        if (abs(ΔE) < Etol) & (Drms < Dtol) & (ite > 5)
            converged = true
            break
        end
    end

    @output repeat("-",80)*"\n"
    @output "    RHF done in {:>5.2f}s\n" t
    @output "    @E[RHF] = {:>20.16f}\n" E
    @output "\n   • Orbitals Summary\n"
    @output "\n {:>10}   {:>15}   {:>10}\n" "Orbital" "Energy" "Occupancy"
    for i in eachindex(eps)
        @output " {:>10}   {:> 15.10f}   {:>6}\n" i eps[i] (i ≤ ndocc ? "↿⇂" : "")
    end
    @output "\n"
    if converged
        @output "   ✔  SCF Equations converged 😄\n"
    else
        @output "❗ SCF Equations did not converge in {:>5} iterations ❗\n" maxit
    end
    @output repeat("-",80)*"\n"

    #occ = CanonicalOrbitals([CanonicalOrbital(Array{Float64,1}(C[:,i])) for i in 1:ndocc])
    #vir = CanonicalOrbitals([CanonicalOrbital(Array{Float64,1}(C[:,i])) for i in ndocc+1:ndocc+nvir])
    all = CanonicalOrbitals([CanonicalOrbital(Array{Float64,1}(C[:,i])) for i in 1:ndocc+nvir])
    aoint.orbs.ndocc = ndocc
    aoint.orbs.nvir = nvir
    aoint.orbs.frozencore = Fermi.CurrentOptions["drop_occ"]
    aoint.orbs.frozenvir = Fermi.CurrentOptions["drop_vir"]
    aoint.orbs["canonical"] = all
    aoint["F"]   = F

    return RHF(molecule, E, ndocc, nvir, eps, aoint)
end

function build_fock!(F::Array{Float64,2}, H::Array{Float64,2}, D::Array{Float64,2}, ERI::Array{Float64,4}, Co)
    F .= H
    Fermi.contract!(F,D,ERI,1.0,1.0,2.0,"mn","rs","mnrs")
    Fermi.contract!(F,D,ERI,1.0,1.0,-1.0,"mn","rs","mrns")
end

function build_fock!(F::Array{Float64,2}, H::Array{Float64,2}, D::Array{Float64,2}, ERI::Array{Float64,3}, Co, mn, J;build_K=true)
    F .= H
    sz = size(F,1)
    dfsz = size(ERI,1)
    no = size(Co,2)

    Fp = Array{Float64}(undef,dfsz)
    tj = @elapsed begin
        Fermi.contract!(Fp,D,ERI,0.0,1.0,2.0,"Q","rs","Qrs")
    end
    #Fermi.contract!(F,Fp,ERI,1.0,1.0,1.0,"mn","Q","Qmn")

    ρ = zeros(dfsz,sz,no)
    t1 = @elapsed begin
        Threads.@threads for m=1:sz
            for r=1:sz
                if abs(mn[m,r])^(1/2) > 1E-10
                s = 0.0
                @avx for Q=1:dfsz
                    eqmr = ERI[Q,m,r]
                    for c=1:no
                        ρ[Q,m,c] += Co[r,c]*eqmr
                    end
                    s += Fp[Q]*eqmr
                end
                F[m,r] += s
                end
            end
        end
        #Fermi.contract!(ρ,Co,ERI,0.0,1.0,1.0,"Qmc","rc","Qmr")
    end
    t2 = @elapsed Fermi.contract!(F,ρ,ρ,1.0,1.0,-1.0,"mn","Qmc","Qnc")
    @output "{} {}\n" tj t1+t2
end

function build_J!(J::Array{Float64,2},D::Array{Float64,2},ERI::Array{Float64,3})
    sz = size(ERI,2)
    dfsz = size(ERI,1)
    Fp = zeros(dfsz)
    Fermi.contract!(Fp,D,ERI,1.0,1.0,2.0,"Q","rs","Qrs")
    Fermi.contract!(J,Fp,ERI,1.0,1.0,1.0,"mn","Q","Qmn")
end

function build_K!(K::Array{Float64,2},D::Array{Float64,2},ERI::Array{Float64,3})
    sz = size(ERI,2)
    dfsz = size(ERI,1)
    Fp = zeros(dfsz,sz,sz)
    Fermi.contract!(Fp,D,ERI,0.0,1.0,-1.0,"Qrn","rs","Qns")
    Fermi.contract!(K,ERI,Fp,1.0,1.0,1.0,"mn","Qmr","Qrn")
end


text:
using Test
using SparseArrays
using DualDecomposition

const DD = DualDecomposition
const parallel = DD.parallel

parallel.init()

@testset "combine_dict with $(parallel.nprocs()) procs" begin
    x = Dict{Int,Float64}()
    xlen = 3
    for i = 1:xlen
        x[xlen*parallel.myid()+i] = xlen*parallel.myid()+i*0.2
    end

    x_combined = parallel.combine_dict(x)
    if parallel.is_root()
        for j = 1:parallel.nprocs(), i = 1:xlen
            @test x_combined[xlen*(j-1)+i] == xlen*(j-1)+i*0.2
        end
    end

    y = Dict{Int,SparseVector{Float64}}()
    ylen = 3
    for i = 1:ylen
        y[ylen*parallel.myid()+i] = [ylen*parallel.myid()+i*0.2+j for j=1:3]
    end
    y_combined = parallel.combine_dict(y)
    if parallel.is_root()
        for id = 1:parallel.nprocs(), i = 1:ylen
            @test y_combined[ylen*(id-1)+i] == sparsevec([ylen*(id-1)+i*0.2+j for j=1:3])
            @test typeof(y_combined[ylen*(id-1)+i]) == SparseVector{Float64,Int64}
        end
    end
end

parallel.finalize()


text:
#!/usr/bin/env julia

## =========================== ## ## =========================== ##
## 	   Euler's and Heun's methods for numerical integration 	 ##
##	   Ollin Demian Langle Chimal					             ##
##	                         						             ##
##     09 / 23 / 2018									         ##
## =========================== ## ## =========================== ##

"""
    euler(x,h;fun,kwargs...)

Uses Euler method to integrate function fun at x with step h
"""
function euler(x,h;fun=lv,kwargs...)
    try
        f = fun(x;kwargs...)
        y = x + f * h
        return y
    catch
        return @. x + fun(x;kwargs...) * h
    end
end

"""
    heun(x,h;fun,kwargs...)

Uses Heun method to integrate function fun at x with step h
"""
function heun(x,h;fun=lv,kwargs...)
    try
        y1 = euler(x,h;fun=fun,kwargs...)
        y2 = fun(x;kwargs...) .+ fun(y1;kwargs...)
        yf = x .+ h * 0.5 * y2
        return yf
    catch
        y1 = euler(x,h;fun=fun,kwargs...)
        y2 = @. fun(x;kwargs...) + fun(y1;kwargs...)
        yf = @. x + h * 0.5 * y2
        return yf
    end
end

"""
    lv(x;a,b,c,d)

Classical 2D Lotka-Volterra model
"""
function lv(x;a,b,c,d)
    y = [a*x[1]-b*x[1]*x[2],c*x[1]*x[2]-d*x[2]]
    y
end

"""
    lot3(x;α)

Generalized Lotka-Volterra model for 3 species as in
Flake's The Computational Beauty of Nature book.
"""
function lot3(x;α=0.75)
    A = [0.5 0.5 0.1; -0.5 -0.1 0.1; α 0.1 0.1]
    y = zeros(length(x))
    for i in 1:length(x)
        y[i] = x[i] * (A[i,:]' * (1 .- x))
    end
    y
end

"""
    my_integrator(ti,tf,h,y0;fun,algo,kwargs...)

Integrates a function fun from ti to tf with step h and initial condition
y0 with algorithm algo.
"""
function my_integrator(ti,tf,h,y0;fun=lv,algo=euler,kwargs...)
    t = ti:h:tf |> collect
    y = zeros(length(t),length(y0))
    y[1,:] .= y0
    for i ∈ 2:size(y)[1]
        value = algo(y[i-1,:],h;fun=fun,kwargs...)
        if any(@. !isless(value,1e4))
            y = y[1:i-1,:]
            t = t[1:size(y)[1]]
            break
        end
        y[i,:] .= value
    end
    t,y
end

# Test the integration in a Lotka-Volterra 2 species system and for a regular
# function like cosine.

#  @time t,y = my_integrator(0,500,1/2^4,[2/3,0.9];fun=lv,algo=heun,a=0.5,b=1.2,c=1,d=0.8)
#  @time t,y = my_integrator(0,500,0.125,[2/3,0.9];fun=cos,algo=heun)


text:
"""
Module to house an interface of vectors commonly used in physics.
"""
module PhysicalVectors

# Exports from PhysicalVector interface
export PhysicalVector, equal, add!, multiply!, ⋅, subtract!, divide!, dot, magnitude, magnitude2, normalize, normalize!, unit, unit!

"""
This is an interface for vectors used in physics.

# Additional Information

The only functions that the user must define for a
`VectorType <: PhysicalVector` are:

* Equality:                       `equal(vector1, vector2)`
* Addition:                       `Base.:+(vector1, vector2)`
* In-place Addition:              `add!(vector1, vector2)`
* In-place Subtraction:           `subract!(vector1, vector2)`
* Scalar Multiplication:          `Base.:*(scalar, vector)`
* In-place Scalar Multiplication: `multiply!(vector, scalar)`
* Scalar Product:                 `⋅(vector1, vector2) #\\cdot [tab]`

The other exported functions are standard vector
operations which can be written in terms of those
above:

* Scalar Multiplication:          `Base.:*(vector, scalar)`
* Subtraction:                    `Base.:-(vector1, vector2)`
* Divison by Scalar:              `Base.:/(vector, scalar)`
* In-place Divison by Scalar:     `divide!(vector, scalar)`
* Scalar Product:                 `dot(vector1, vector2)`
* Square Magnitude:               `magnitude2(vector)`
* Magnitude:                      `magnitude(vector)`
* Normalization:                  `normalize(vector)`
* In-place Normalization:         `normalize!(vector)`
* Unit:                           `unit(vector)`
* In-place Unit:                  `unit!(vector)`
"""
abstract type PhysicalVector end

"""
# REQUIRED FOR `PhysicalVector` INTERFACE

    equal(A::PhysicalVector, B::PhysicalVector)

Check for equality between two `PhysicalVector`s.

# Additional Information
* Note that Base.:(==) is not used for issues with mutability and equality versus egality.
* This function `error`s out in the default implementation to enforce its definition in subtypes of `PhysicalVector`.
"""
equal(A::PhysicalVector, B::PhysicalVector) = error("No implementation defined for vectors of type $(typeof(A)) and $(typeof(B)).")

"""
# REQUIRED FOR `PhysicalVector` INTERFACE

    +(A::PhysicalVector, B::PhysicalVector)

Addtion of two `PhysicalVector`s.

# Additional Information
* A call to a copy constructor is implied.
* This function `error`s out in the default implementation to enforce its definition in subtypes of `PhysicalVector`.
"""
Base.:+(A::PhysicalVector, B::PhysicalVector) = error("No implementation defined for vectors of type $(typeof(A)) and $(typeof(B)).")

"""
# REQUIRED FOR `PhysicalVector` INTERFACE

    add!(A::PhysicalVector, B::PhysicalVector)

In-place addition of two `PhysicalVector`s.

# Additional Information
* It is implied that the first argument, `A`, is modified in-place.
* This function `error`s out in the default implementation to enforce its definition in subtypes of `PhysicalVector`.
"""
add!(A::PhysicalVector, B::PhysicalVector) = error("No implementation defined for vectors of type $(typeof(A)) and $(typeof(B)).")

"""
# REQUIRED FOR `PhysicalVector` INTERFACE

    *(λ, A::PhysicalVector)

Multiplication of a `PhysicalVector` by the scalar `λ`.

# Additional Information
* A call to a copy constructor is implied.
* This function `error`s out in the default implementation to enforce its definition in subtypes of `PhysicalVector`.
"""
Base.:*(λ, A::PhysicalVector) = error("No implementation defined for vectors of type $(typeof(λ)) and $(typeof(A)).")

"""
# REQUIRED FOR `PhysicalVector` INTERFACE

    multiply!(A::PhysicalVector, λ)

In-place multiplication of a `PhysicalVector` by the scalar `λ`.

# Additional Information
* It is implied that the first argument, `A`, is modified in-place.
* This function `error`s out in the default implementation to enforce its definition in subtypes of `PhysicalVector`.
"""
multiply!(A::PhysicalVector, λ) = error("No implementation defined for vectors of type $(typeof(λ)) and $(typeof(A)).")

"""
# REQUIRED FOR `PhysicalVector` INTERFACE

    ⋅(A::PhysicalVector, B::PhysicalVector) #\\cdot [tab]

Scalar product of two `PhysicalVector`s. This is also referred to as the `dot` product.

# Additional Information
* This function `error`s out in the default implementation to enforce its definition in subtypes of `PhysicalVector`.
"""
⋅(A::PhysicalVector, B::PhysicalVector) = error("No implementation defined for vectors of type $(typeof(A)) and $(typeof(B)).")

#= ========================================================================================= =#
#  These are functions that follow the interface but are not
#  necessary to define for subtypes of PhysicalVector
#= ========================================================================================= =#
"""
    *(A::PhysicalVector, λ)

Multiplication of a `PhysicalVector` by the scalar `λ`.

# Additional Information
* This is an alias for `*(λ, A::PhysicalVector)`.
"""
Base.:*(A::PhysicalVector, λ) = λ * A

"""
    -(A::PhysicalVector, B::PhysicalVector)

Subtraction of two `PhysicalVector`s.

# Additional Information
* A copy constructor call is implied.
"""
Base.:-(A::PhysicalVector, B::PhysicalVector) = A + (-1 * B)

"""
    /(A::PhysicalVector, λ)

Division of a `PhysicalVector` by the scalar `λ`.

# Additional Information
* A copy constructor call is implied.
"""
Base.:/(A::PhysicalVector, λ) = A * (1 / λ)

"""
    divide!(A::PhysicalVector, λ)

Division of a `PhysicalVector` by the scalar `λ`.

# Additional Information
* It is implied that the first argument, `A`, is modified in-place.
"""
divide!(A::PhysicalVector, λ) = multiply!(A, 1 / λ)

"""
    dot(A::PhysicalVector, B::PhysicalVector)

Scalar product of two `PhysicalVector`s.

# Additional Information
* This is an alias call to `⋅(A::PhysicalVector, B::PhysicalVector) #\\cdot [tab]`
"""
dot(A::PhysicalVector, B::PhysicalVector) = A ⋅ B

"""
    magnitude2(A::PhysicalVector)

The square magnitude of a `PhysicalVector`.

# Additional Information
* This is an alias call to `⋅(A::PhysicalVector, A::PhysicalVector) #\\cdot [tab]`
"""
magnitude2(A::PhysicalVector) = A ⋅ A

"""
    magnitude(A::PhysicalVector)

The magnitude of a `PhysicalVector`.

# Additional Information
* This is an alias call to `sqrt( magnitude2(A::PhysicalVector) )`
"""
magnitude(A::PhysicalVector) = sqrt(magnitude2(A))

"""
    normalize(A::PhysicalVector)

Divide a `PhysicalVector` by its length.

# Additional Information
* A call to a copy constructor is implied.
"""
normalize(A::PhysicalVector) = A / magnitude(A)

"""
    normalize!(A::PhysicalVector)

Divide a `PhysicalVector` by its length.

# Additional Information
* It is implied that the argument, `A`, is modified in-place.
"""
normalize!(A::PhysicalVector) = divide!(A, magnitude(A))

"""
    unit(A::PhysicalVector)

Divide a `PhysicalVector` by its length.

# Additional Information
* This is an alias for `normalize(A::PhysicalVector)`.
* A call to a copy constructor is implied.
"""
unit(A::PhysicalVector) = normalize(A)

"""
    unit!(A::PhysicalVector)

Divide a `PhysicalVector` by its length.

# Additional Information
* This is an alias for `normalize!(A::PhysicalVector)`.
* It is implied that the argument, `A`, is modified in-place.
"""
unit!(A::PhysicalVector) = normalize!(A)

#= ========================================================================================= =#
#  Include implementations of the PhysicalVector interface.
#  Each implementation file is included and contains it's own export list.
#= ========================================================================================= =#
include("Vectors2D.jl")

end # module PhysicalVectors


text:
#include("../dataset_reader.jl")
const dataset1_labels = Float64[1, 2, 3, 4, 5]
const dataset1_nfeatures = 5
const dataset1_filename = "drug200.csv"

function get_numerical_dataset1(path::String="")
    data = get_data(joinpath(path, dataset1_filename))

    sex_dict = create_numerical_map(getindex.(data, 2))
    bp_dict = create_numerical_map(getindex.(data, 3))
    cholesterol_dict = create_numerical_map(getindex.(data, 4))
    drug_dict = create_numerical_map(getindex.(data, 6))

    conversion_pattern = [Float64, sex_dict, bp_dict, cholesterol_dict, Float64, drug_dict]
    convert_values(data, conversion_pattern)
end

function get_dataset1(path::String=""; shuffle::Bool=false, test_perc::Float64=0.2,
                      normalise::Bool=false, separate_test_data::Bool=true, balance_classes::Bool=false)
    data = get_numerical_dataset1(path)
    separate_dataset(data, test_percentage=test_perc, shuffle=shuffle, normalise=normalise,
                     separate_test_data=separate_test_data, balance_classes=balance_classes)
end

function get_dataset1_setup(path::String=""; normalise::Bool=false, separate::Bool=false,
                            balance::Bool=false)
    dataset = get_dataset1(path, normalise=normalise, separate_test_data=separate, balance_classes=balance)
    if separate
        SeparatedDataset(dataset..., dataset1_labels, dataset1_nfeatures, 2)
    else
        Dataset(dataset..., dataset1_labels, dataset1_nfeatures, 2)
    end
end

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
type ZeroElementException <: Exception end

#macro debug_stmt(x)
#  x
#end

macro debug_stmt(x)
end

# used to send some size information.
immutable TableSize{N}
  size::NTuple{N,Int}
end
TableSize(arr::AbstractArray) = TableSize(size(arr))

# some statements for compatibility.
const IS_JULIA_V05 = startswith(string(VERSION), "0.5") || startswith(string(VERSION), "0.6")

iosize_expr = if IS_JULIA_V05
  :(iosize_compat = Base.displaysize)
  #:(iosize_compat = Base.iosize)
else
  :(iosize_compat = Base.tty_size)
end

subarray_last_type = if IS_JULIA_V05
  :(const SUBARRAY_LAST_TYPE = true)
else
  :(const SUBARRAY_LAST_TYPE = 1)
end

eval(iosize_expr)
eval(subarray_last_type)


text:
using DiscreteFiltering
using LinearAlgebra
using SparseArrays
using Test

@testset "Discrete filtering" begin
    # Unknown domain type
    struct UnknownDomain <: DiscreteFiltering.Domain end

    @testset "Structures" begin
        include("structures/domain.jl")
        include("structures/filter.jl")
    end

    @testset "Matrix assembly" begin
        include("matrices/advection_matrix.jl")
        include("matrices/diffusion_matrix.jl")
        include("matrices/filter_matrix.jl")
        include("matrices/reconstruction_matrix.jl")
    end

    @testset "Equations" begin
        include("equations/equations.jl")
        include("equations/solve.jl")
        include("equations/solve_diffusion.jl")
        include("equations/solve_advection.jl")
        include("equations/solve_burgers.jl")
        include("equations/solve_adbc.jl")
    end

    @testset "Utils" begin
        include("utils/ridge.jl")
    end
end