text:
# Legit.jl -- Convert French law in LEGI format to Git & Markdown
# By: Emmanuel Raviart <emmanuel.raviart@data.gouv.fr>
#
# Copyright (C) 2015 Etalab
# https://github.com/etalab/Legit.jl
#
# The Legit.jl package is licensed under the MIT "Expat" License.

module Legit

using ArgParse
using Biryani
using Biryani.DatesConverters
using Compat
using Dates: Date, DateTime, datetime2unix, Day, year
using LibGit2
using LightXML
import Slugify

include("converters.jl")
include("articles.jl")

function main()
  args = parse_command_line()
  mode = args["mode"]
  readme_mode = args["readme"]

  @assert splitext(args["repository"])[end] == ".git"
  if ispath(args["repository"])
    @assert isdir(args["repository"])
    @assert iswritable(args["repository"])
    if args["erase"]
      rm(args["repository"]; recursive = true)
    end
  end
  if ispath(args["repository"])
    repository = repo_discover(args["repository"])
    head_reference = head(repository)
    latest_commit_oid = target(head_reference)
    latest_commit = lookup_commit(repository, latest_commit_oid)
    root_tree = GitTree(latest_commit)
  else
    mkpath(args["repository"])
    repository = init_repo(args["repository"]; bare = true)
    latest_commit = nothing
    root_tree = nothing
  end

  if mode == "all"
    documents_dir_walker = @task walk_documents_dir(joinpath(args["legi_dir"], "global"))
    node_by_nature = @compat Dict{String, SimpleNode}()
    node_by_year_by_container = @compat Dict{SimpleNode, Dict{String, SimpleNode}}()
    root_title = "Lois et règlements français"
  elseif mode == "codes"
    documents_dir_walker = @task walk_documents_dir(
      joinpath(args["legi_dir"], "global", "code_et_TNC_en_vigueur", "code_en_vigueur"),
      joinpath(args["legi_dir"], "global", "code_et_TNC_non_vigueur", "code_non_vigueur"),
    )
    root_title = "Codes juridiques français"
  else
    @assert mode == "non-codes"
    documents_dir_walker = @task walk_documents_dir(
      joinpath(args["legi_dir"], "global", "code_et_TNC_en_vigueur", "TNC_en_vigueur"),
      joinpath(args["legi_dir"], "global", "code_et_TNC_non_vigueur", "TNC_non_vigueur"),
    )
    node_by_nature = @compat Dict{String, SimpleNode}()
    node_by_year_by_container = @compat Dict{SimpleNode, Dict{String, SimpleNode}}()
    root_title = "Lois non codifiées et règlements français"
  end
  codes_en_vigueur_node = nothing
  codes_non_vigueur_node = nothing
  root_node = RootNode(root_title)
  skip_documents = args["start"] !== nothing
  for (document_index, document_dir) in enumerate(documents_dir_walker)
    if args["only"] !== nothing && basename(document_dir) != args["only"]
      continue
    end
    if skip_documents
      document_cid = basename(document_dir)
      if document_cid == args["start"]
        skip_documents = false
      else
        continue
      end
    end

println()
println()
println("=============================================================================================================")
println(document_index, " / ", document_dir)

    root_section = root_tree === nothing ?
      Section(root_node.title) :
      parse_section_commonmark(repository, entry_bypath(root_tree, "README.md"), root_node.title)

    version_dir = joinpath(document_dir, "texte", "version")
    version_filenames = sort(readdir(version_dir))
    struct_dir = joinpath(document_dir, "texte", "struct")
    struct_filenames = sort(readdir(struct_dir))
    @assert(length(struct_filenames) == length(version_filenames),
      "Directory $struct_dir doesn't contain the same number of files as directory: $struct_filenames.")

    articles_by_id = @compat Dict{String, Vector{Article}}()  # Articles are sorted by start date for each ID.
    changed_by_message_by_date = @compat Dict{Date, Dict{String, Changed}}()
    notes = nothing
    signataires = nothing
    visas = nothing
    for (version_filename, struct_filename) in zip(version_filenames, struct_filenames)
      version_xml_document = parse_file(joinpath(version_dir, version_filename))
      texte_version = Convertible(parse_xml_element(root(version_xml_document))) |> pipe(
        element_to_texte_version,
        require,
      ) |> to_value
      # free(version_xml_document)

      @assert(version_filename == struct_filename,
        "Filenames for struct and version differ: $struct_filename != $version_filename.")
      struct_xml_document = parse_file(joinpath(struct_dir, struct_filename))
      textelr = Convertible(parse_xml_element(root(struct_xml_document))) |> pipe(
        element_to_textelr,
        require,
      ) |> to_value
      # free(struct_xml_document)

      if mode == "all"
        nature = get(texte_version["META"]["META_COMMUN"], "NATURE", "nature_inconnue")
        document_container = get!(node_by_nature, nature) do
          return SimpleNode(root_node, nature)
        end
      elseif mode == "codes"
        document_container = root_node
      else
        @assert mode == "non-codes"
        nature = get(texte_version["META"]["META_COMMUN"], "NATURE", "nature_inconnue")
        document_container = get!(node_by_nature, nature) do
          return SimpleNode(root_node, nature)
        end
      end
      if contains(document_dir, "code_en_vigueur")
        if codes_en_vigueur_node === nothing
          codes_en_vigueur_node = SimpleNode(document_container, "Codes en vigueur")
        end
        document_container = codes_en_vigueur_node
      elseif contains(document_dir, "code_non_vigueur")
        if codes_non_vigueur_node === nothing
          codes_non_vigueur_node = SimpleNode(document_container, "Codes non en vigueur")
        end
        document_container = codes_non_vigueur_node
      else
        node_by_year = get!(node_by_year_by_container, document_container) do
          return @compat Dict{String, SimpleNode}()
        end
        start_date = get(texte_version["META"]["META_SPEC"]["META_TEXTE_VERSION"], "DATE_DEBUT", nothing)
        year_text = start_date === nothing ? "Année inconnue" : string(year(start_date))
        document_container = get!(node_by_year, year_text) do
          return SimpleNode(document_container, year_text)
        end
      end
      document = Document(document_container, texte_version, textelr)
@show node_title(document)
println("=============================================================================================================")
println()

      if notes === nothing
        notes_dict = get(document.texte_version, "NOTA", nothing)
        if notes_dict !== nothing
          notes = NonArticle(document, "notes", "Nota", notes_dict["CONTENU"])
        end
      end
      if signataires === nothing
        signataires_dict = get(document.texte_version, "SIGNATAIRES", nothing)
        if signataires_dict !== nothing
          signataires = NonArticle(document, "signataires", "", signataires_dict["CONTENU"])
        end
      end
      if visas === nothing
        visas_dict = get(document.texte_version, "VISAS", nothing)
        if visas_dict !== nothing
          visas = NonArticle(document, "visas", "", visas_dict["CONTENU"])
        end
      end

      parse_structure(document, articles_by_id, changed_by_message_by_date, document_dir)
    end
    link_articles(articles_by_id)

    for (date_index, date) in enumerate(sort(collect(keys(changed_by_message_by_date))))
      if args["last-date"] !== nothing && date > Date(args["last-date"])
        continue
      end
      changed_by_message = changed_by_message_by_date[date]
      for (message_index, message) in enumerate(sort(collect(keys(changed_by_message))))
        changed = changed_by_message[message]
println("-------------------------------------------------------------------------------------------------------------")
println("$date $message")
        epoch = int64(datetime2unix(DateTime(date)))
        if args["date"]
          # Modify dates before 1970-02-01 (epoch = 2678340) to ensure that they are >= 1970-01-01 (epoch = 0).
          if epoch < 2678340  # = 31 * 24 * 60 * 60 = 1970-02-01
            epoch /= 24 * 60 * 60  # Days become seconds
            epoch += 2678340 - 31  # 31 = days of 1970-01
          end
        end
        time_offset = 0
        author_signature = Signature("République française", "gitloi@data.gouv.fr", epoch, time_offset)
        committer_signature = author_signature
        commit_needed = false

        if date_index == 1 && message_index == 1
          if visas !== nothing
            unshift!(changed.articles, visas)
          end
          if signataires !== nothing
            push!(changed.articles, signataires)
          end
          if notes !== nothing
            push!(changed.articles, notes)
          end
        end

        # Update sections tree.
        for article in changed.articles
if isa(article, Article)
  println("Upserted: ", node_id(article), " ", node_git_file_path(article), ".")
else
  println("Upserted: non article ", node_git_file_path(article), ".")
end
          nodes = Node[]
          container = article.container
          while true
            unshift!(nodes, container)
            container = container.container
            if isa(container, RootNode)
              break
            end
          end
          section = root_section
          section_names = String[]
          for node in nodes
            dir_name = node_dir_name(node)
            push!(section_names, dir_name)
            child_section = get!(section.child_by_name, dir_name) do
              return Section()
            end
            if isa(child_section, UnparsedSection)
              child_section = parse_section_commonmark(repository,
                entry_bypath(root_tree, string(join(section_names, '/'), '/', "README.md")), child_section.short_title)
              section.child_by_name[dir_name] = child_section
            end
            section = child_section
            section.short_title = node_short_title(node)
            section.sortable_title = node_sortable_title(node)
            section.title = node_title(node)
          end
          section.child_by_name[node_name(article)] = article
        end
        deleted_articles = Article[]
        for article in changed.deleted_articles
          next_version = get(article.next_version, article)
          if next_version != article && node_git_file_path(article) == node_git_file_path(next_version)
            # Article is neither deleted nor moved => Keep it unchanged.
            continue
          end
println("Deleted: ", string(node_id(article), " ", node_git_file_path(article)))
          nodes = Node[]
          container = article
          while true
            unshift!(nodes, container)
            container = container.container
            if isa(container, RootNode)
              break
            end
          end
          section = root_section
          section_names = String[]
          sections = Node[]
          for node in nodes
            push!(sections, section)
            name = node_name(node)
            push!(section_names, name)
            child_section = get(section.child_by_name, name, nothing)
            if child_section === nothing
              warn("Unknown sub-section $(node_name(node))")
              sections = Node[]
              break
            elseif isa(child_section, UnparsedSection)
              child_section = parse_section_commonmark(repository,
                entry_bypath(root_tree, string(join(section_names, '/'), '/', "README.md")), child_section.short_title)
              section.child_by_name[name] = child_section
            end
            section = child_section
          end
          if !isempty(sections)
            push!(deleted_articles, article)
            for (node, section) in zip(reverse(nodes), reverse(sections))
              delete!(section.child_by_name, node_name(node))
              if !isempty(section.child_by_name)
                break
              end
            end
          end
        end

        if readme_mode == "single-page"
          root_tree_builder = GitTreeBuilder(repository, root_tree)
          single_section = first(first(root_section.child_by_name)[2].child_by_name)[2]
          single_section_filename = node_filename(single_section)
          blob = blob_from_buffer(repository, commonmark(single_section, readme_mode))
          blob_id = Oid(blob)
          entry = root_tree_builder[single_section_filename]
          if entry === nothing || Oid(entry) != blob_id
            insert!(root_tree_builder, single_section_filename, blob_id, int(0o100644))  # FILEMODE_BLOB
            commit_needed = true
            root_tree_oid = write!(root_tree_builder)
          end
        else
          tree_builder_by_git_dir_names = @compat Dict{Tuple, GitTreeBuilder}()

          for article in changed.articles
            git_dir_names = tuple(split(node_git_dir(article), '/')...)
            tree_builder = get!(tree_builder_by_git_dir_names, git_dir_names) do
              if root_tree === nothing
                latest_tree = nothing
              elseif isempty(git_dir_names)
                latest_tree = root_tree
              else
                entry = entry_bypath(root_tree, join(git_dir_names, '/'))
                latest_tree = entry === nothing ? nothing : lookup_tree(repository, Oid(entry))
              end
              return GitTreeBuilder(repository, latest_tree)
            end
            article_filename = node_filename(article)
            blob = blob_from_buffer(repository, commonmark(article, readme_mode))
            blob_id = Oid(blob)
            entry = tree_builder[article_filename]
            if entry === nothing || Oid(entry) != blob_id
              insert!(tree_builder, article_filename, blob_id, int(0o100644))  # FILEMODE_BLOB
              commit_needed = true
            end
          end

          for article in deleted_articles
            git_dir_names = tuple(split(node_git_dir(article), '/')...)
            tree_builder = get!(tree_builder_by_git_dir_names, git_dir_names) do
              if root_tree === nothing
                latest_tree = nothing
              elseif isempty(git_dir_names)
                latest_tree = root_tree
              else
                entry = entry_bypath(root_tree, join(git_dir_names, '/'))
                latest_tree = entry === nothing ? nothing : lookup_tree(repository, Oid(entry))
              end
              return GitTreeBuilder(repository, latest_tree)
            end
            article_filename = node_filename(article)
            entry = tree_builder[article_filename]
            if entry !== nothing
              delete!(tree_builder, article_filename)
              commit_needed = true
            end
          end

          if commit_needed
            root_tree_oid = nothing
            while !isempty(tree_builder_by_git_dir_names)
              git_dirs_names_to_build = Set(keys(tree_builder_by_git_dir_names))
              for (git_dir_names, tree_builder) in tree_builder_by_git_dir_names
                while !isempty(git_dir_names)
                  git_dir_names = git_dir_names[1 : end - 1]
                  pop!(git_dirs_names_to_build, git_dir_names, nothing)
                end
              end
              for git_dir_names in git_dirs_names_to_build
                tree_builder = pop!(tree_builder_by_git_dir_names, git_dir_names)
                if length(tree_builder) == 0 || readme_mode in ("deep", "flat") && length(tree_builder) == 1 &&
                    tree_builder["README.md"] !== nothing
                  tree_oid = nothing
                else
                  if readme_mode in ("deep", "flat")
                    section = root_section
                    for dir_name in git_dir_names
                      section = section.child_by_name[dir_name]
                    end
                    blob = blob_from_buffer(repository, commonmark(section, readme_mode))
                    insert!(tree_builder, "README.md", Oid(blob), int(0o100644))  # FILEMODE_BLOB
                  end
                  tree_oid = write!(tree_builder)
                end
                if isempty(git_dir_names)
                  root_tree_oid = tree_oid
                else
                  dir_name = git_dir_names[end]
                  git_dir_names = git_dir_names[1 : end - 1]
                  @assert !(git_dir_names in git_dirs_names_to_build)
                  tree_builder = get!(tree_builder_by_git_dir_names, git_dir_names) do
                    if root_tree === nothing
                      latest_tree = nothing
                    elseif isempty(git_dir_names)
                      latest_tree = root_tree
                    else
                      entry = entry_bypath(root_tree, join(git_dir_names, '/'))
                      latest_tree = entry === nothing ? nothing : lookup_tree(repository, Oid(entry))
                    end
                    return GitTreeBuilder(repository, latest_tree)
                  end
                  if tree_oid === nothing
                    delete!(tree_builder, dir_name)
                  else
                    insert!(tree_builder, dir_name, tree_oid, int(0o40000))  # FILEMODE_TREE
                  end
                end
              end
            end
          end
        end

        if commit_needed
          root_tree = lookup_tree(repository, root_tree_oid)
          latest_commit_oid = commit(repository, "HEAD", message, author_signature, committer_signature,
            root_tree, (latest_commit === nothing ? [] : [latest_commit])...)
          latest_commit = lookup_commit(repository, latest_commit_oid)
        end
      end
    end
  end
end

function parse_command_line()
  arg_parse_settings = ArgParseSettings()
  @add_arg_table arg_parse_settings begin
    "--date", "-d"
      action = :store_true
      help = "modify commit dates before 1970-02-01 to ensure that they are after 1970-01-01"
    "--erase", "-e"
      action = :store_true
      help = "erase existing Git repository"
    "--last-date", "-l"
      help = "remove commits that are too much in the future (more than 25 years) for compatibility with GitLab"
      range_tester = value -> Convertible(value) |> iso8601_input_to_date |> require |> is_valid
    "--mode", "-m"
      default = "all"
      help = "mode for generated tree of files in Git repository (all, codes)"
      range_tester = value -> value in ("all", "codes", "non-codes")
    "--only", "-o"
      help = "CID of single LEGI document to parse"
      range_tester = value -> Convertible(value) |> validate_cid |> require |> is_valid
    "--readme", "-r"
      default = "flat"
      help = "mode for README file (deep, flat, none, single-page)"
      range_tester = value -> value in ("deep", "flat", "none", "single-page")
    "--start", "-s"
      help = "CID of first LEGI document to parse"
      range_tester = value -> Convertible(value) |> validate_cid |> require |> is_valid
    "--verbose", "-v"
      action = :store_true
      help = "increase output verbosity"
    "legi_dir"
      help = "path of LEGI dir"
      required = true
    "repository"
      help = "path of Git repository to create"
      required = true
  end
  return parse_args(arg_parse_settings)
end

function walk_documents_dir(dir::String)
  filenames = readdir(dir)
  is_texte_version = endswith(dir, joinpath("texte", "version"))
  texte_version_files_found = false
  for filename in filenames
    file_path = joinpath(dir, filename)
    if isdir(file_path)
      walk_documents_dir(file_path)
    elseif is_texte_version && isfile(file_path)
      texte_version_files_found = true
    end
  end
  if texte_version_files_found
    document_dir = dirname(dirname(dir))
    produce(document_dir)
  end
end

function walk_documents_dir(dirs::String...)
  for dir in dirs
    walk_documents_dir(dir)
  end
end

main()

end # module


text:
# Question
# By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms
# 0, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

# Time
# O(n) any better solution??
# every third term is even (neglect 0 ) we can use it and remove if
# any better solution?

function main()
    limit = 4_000_000 # 4 million
    # series odd odd even types
    s = 2
    a = [1,2]
    b = 2
    i = 3
    while (b<limit)
        b = a[i-1]+a[i-2]
        push!(a, b)
        if(b %2 == 0)
            s+=b
        end
        i+=1
    end
    println(s)
end

@time main()

ext: jl


text:
using ElectronProgressBars

ElectronProgressBars._gset_progress(0.0, title="main")
ElectronProgressBars._gset_progress(0.1, :sub, title="main sub")
ElectronProgressBars._gset_progress(0.8, :subsub, title="main sub sub")

ElectronProgressBars._gset_progress(0.1, title="bg", taskid=UInt64(2))
ElectronProgressBars._gset_progress(0.9, :sub; title="bg sub", taskid=UInt64(2))


text:
using TraitSimulation
using GLM
using LinearAlgebra
using Random
using Test
using DataFrames
using VarianceComponentModels
using SnpArrays
using Suppressor
using DelimitedFiles

@testset "modelparameters" begin
	include("modelparameters.jl")
end

@testset "simulationmodels" begin
	include("simulationmodels.jl")
	include("matrixnormalsimulation.jl")
	include("writetopheno.jl")
end

@testset "simulatepower" begin
	include("vcm_power.jl")
	include("univariate_power.jl")
end


text:
function neighpairlist(x, ilist, neighlist, neighnum, rcutsq)

inum = length(ilist);

pairnum = zeros(eltype(ilist), inum);
pairlist = [];
for ii = 1:inum
    i = ilist[ii]; # atom i
    n1 = neighnum[i];
    n2 = neighnum[i+1];

    g = neighlist[(n1+1):n2]; # a list of neighbors around atom i
    d = x[:,g] .- x[:,i]; # distance between neighbors and atom i
    dsq = sum(d.*d, dims=1);
    dsq = dsq[:];
    ind = findall(dsq .<= rcutsq); # find neighbors within rcut from atom i
    # a list of neighbors within rcut from atom i
    pairnum[ii] = length(ind);
    pairlist = [pairlist; g[ind]];
end

pairnum = [0; cumsum(pairnum)];
pairlist = eltype(ilist).(pairlist)

return pairlist, pairnum

end

function neighpairlist(x, ilist, alist, atomtype, neighlist, neighnum, rcutsq, ntypes)

#typej = unique(typej);
inum = length(ilist);

pairnum = zeros(eltype(ilist), inum);
pairlist = [];

for ii = 1:inum
    i = ilist[ii]; # atom i
    n1 = neighnum[i];
    n2 = neighnum[i+1];

    # cut-off radius between atom of type i and atom of type j
    itype = atomtype[i];
    g = neighlist[(n1+1):n2]; # a list of neighbors around atom i
    j = alist[g];
    jtype = atomtype[j]; # types of neighbors around atom i
    rsq = rcutsq[jtype .+ (itype-1)*(ntypes)];

    d = x[:,g] .- x[:,i]; # distance between neighbors and atom i
    dsq = sum(d.*d, dims=1);
    dsq = dsq[:]
    ind = findall(dsq .<= rsq); # find neighbors within rcut from atom i

    # a list of neighbors with typej within rcut from atom i
    pairnum[ii] = length(ind);
    pairlist = [pairlist; g[ind]];
end

pairnum = [0; cumsum(pairnum)];
pairlist = eltype(ilist).(pairlist)

return pairlist, pairnum

end

function neighpairlisttypej(x, ilist, alist, atomtype, typej, neighlist, neighnum, rcutsq)

typej = unique(typej);
inum = length(ilist);

pairnum = zeros(eltype(ilist), inum);
pairlist = [];
for ii = 1:inum
    i = ilist[ii]; # atom i
    n1 = neighnum[i];
    n2 = neighnum[i+1];
    g = neighlist[(n1+1):n2]; # a list of neighbors around atom i
    j = alist[g];
    atype = atomtype[j]; # types of neighbors around atom i
    ind1 = findall(typej .== atype);  # choose neighbors of typej

    d = x[:,g] - x[:,i]; # distance between neighbors and atom i
    dsq = sum(d.*d, dims=1);
    dsq = dsq[:]
    ind2 = findall(dsq .<= rcutsq); # find neighbors within rcut from atom i

    # a list of neighbors with typej within rcut from atom i
    ind = intersect(ind1, ind2);
    pairnum[ii] = length(ind);
    pairlist = [pairlist; g[ind]];
end
pairnum = [0; cumsum(pairnum)];

pairlist = eltype(ilist).(pairlist)

return pairlist, pairnum

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
module test_julia_divisions

using Test

@test repr(5/3)  == "1.6666666666666667"
@test repr(5//3) == "5//3"
@test repr(5\3)  == "0.6"
@test repr(5÷3)  == "1"

#@info :divisions 5/3 5//3 5\3 5÷3
#=
┌ Info: divisions
│   5 / 3 = 1.6666666666666667
│   5 // 3 = 5//3
│   5 \ 3 = 0.6
└   5 ÷ 3 = 1
=#

end # module test_julia_divisions


text:
using Test
using PokerHandEvaluator
using PokerHandEvaluator: evaluate5
using PlayingCards
using BenchmarkTools
const PHE = PokerHandEvaluator

@testset "Straight flush (Ranks 1:10)" begin
    include("test_straight_flush.jl")
end
@testset "Test quads (Ranks 11:166)" begin
    include("test_quads.jl")
end
@testset "Test Full house (Ranks 167:322)" begin
    include("test_full_house.jl")
end
@testset "Test flush (Ranks 323:1599)" begin
    include("test_flush.jl")
end
@testset "Test off-suit straight (Ranks 1600:1609)" begin
    include("test_offsuit_straight.jl")
end
@testset "Test trips (Ranks 1610:2467)" begin
    include("test_trips.jl")
end
@testset "Test two-pair (Ranks 2468:3325)" begin
    include("test_two_pair.jl")
end
@testset "Test pair (Ranks 3326:6185)" begin
    include("test_pair.jl")
end
@testset "Test high card (Ranks 6186:7462)" begin
    include("test_high_card.jl")
end

@testset "N-methods" begin
    N_offsuit = length(PHE.hash_table_offsuit)
    N_flush = length(PHE.hash_table_suited)
    @test N_offsuit+N_flush == 7462
end

@testset "N-hands" begin
    L_1 = length(PHE.straight_ranks()) # flush
    L_2 = length(PHE.quad_ranks())
    L_3 = length(PHE.full_house_ranks())
    L_4 = length(PHE.flush_ranks())
    L_5 = length(PHE.straight_ranks()) # off-suit
    L_6 = length(PHE.trip_ranks())
    L_7 = length(PHE.two_pair_ranks())
    L_8 = length(PHE.pair_ranks())
    L_9 = length(PHE.high_card_ranks())

    @test L_1+L_2+L_3+L_4+L_5+L_6+L_7+L_8+L_9 == 7462
end

@testset "CompactHandEval & FullHandEval (6/7-card evaluate)" begin
    table_cards = (J♡,J♣,A♣,A♢)
    player_cards = (
      (A♠,2♠,table_cards...),
      (J♠,J♣,table_cards...),
    );
    che = CompactHandEval.(player_cards)
    winner = che[argmin(hand_rank.(che))]
    @test hand_rank(winner) == 47
    @test hand_type(winner) == :quads

    table_cards = (J♡,J♣,A♣,A♢,5♣)
    player_cards = (
      (A♠,2♠,table_cards...),
      (J♠,J♣,table_cards...),
    );
    fhe = FullHandEval.(player_cards)
    winner = fhe[argmin(hand_rank.(fhe))]
    @test hand_rank(winner) == 47
    @test hand_type(winner) == :quads
    # could have equally returned (J♠,J♣,J♡,J♣,A♢):
    @test best_cards(winner) == (J♠,J♣,J♡,J♣,A♣)
    @test all_cards(winner) == (J♠,J♣,table_cards...)

    # Additional interfaces
    fhe = FullHandEval((A♠,2♠), table_cards)
    che = CompactHandEval((A♠,2♠), table_cards)
    fhe = FullHandEval((A♠,2♠)..., table_cards...)
    che = CompactHandEval((A♠,2♠)..., table_cards...)

    fhe = FullHandEval((A♠,2♠), table_cards[1:end-1])
    che = CompactHandEval((A♠,2♠), table_cards[1:end-1])
    fhe = FullHandEval((A♠,2♠)..., table_cards[1:end-1]...)
    che = CompactHandEval((A♠,2♠)..., table_cards[1:end-1]...)

    fhe = FullHandEval((A♠,2♠), table_cards[1:end-2])
    che = CompactHandEval((A♠,2♠), table_cards[1:end-2])
    fhe = FullHandEval((A♠,2♠)..., table_cards[1:end-2]...)
    che = CompactHandEval((A♠,2♠)..., table_cards[1:end-2]...)
end

function compute_hand_type_simple(rank::Int)
         if 1 ≤ rank ≤ 10;     return :straight_flush
    elseif 11 ≤ rank ≤ 166;    return :quads
    elseif 167 ≤ rank ≤ 322;   return :full_house
    elseif 323 ≤ rank ≤ 1599;  return :flush
    elseif 1600 ≤ rank ≤ 1609; return :straight
    elseif 1610 ≤ rank ≤ 2467; return :trips
    elseif 2468 ≤ rank ≤ 3325; return :two_pair
    elseif 3326 ≤ rank ≤ 6185; return :one_pair
    elseif 6186 ≤ rank ≤ 7462; return :high_card
    else; error("Bad hand rank input")
    end
end

@testset "Hand types" begin
    ranks = 1:7462
    hrts = compute_hand_type_simple.(ranks)
    hrtbs = PHE.hand_type_binary_search.(ranks)
    @test all(hrts .== hrtbs)
    @test_throws AssertionError PHE.hand_type_binary_search(0)
    @test_throws AssertionError PHE.hand_type_binary_search(7463)
end

if VERSION >= v"1.4"
    @testset "Allocations" begin
        cards_5 = (J♡,J♣,A♣,A♢,5♣);
        cards_7 = (A♠,2♠,cards_5...);
        cards_6 = (J♣,cards_5...);

        # Compile first:
        alloc_1 = @allocated CompactHandEval(cards_5)
        alloc_2 = @allocated CompactHandEval(cards_6)
        alloc_3 = @allocated CompactHandEval(cards_7)
        alloc_4 = @allocated FullHandEval(cards_5)
        alloc_5 = @allocated FullHandEval(cards_6)
        alloc_6 = @allocated FullHandEval(cards_7)

        alloc_1 = @allocated CompactHandEval(cards_5)
        alloc_2 = @allocated CompactHandEval(cards_6)
        alloc_3 = @allocated CompactHandEval(cards_7)
        alloc_4 = @allocated FullHandEval(cards_5)
        alloc_5 = @allocated FullHandEval(cards_6)
        alloc_6 = @allocated FullHandEval(cards_7)

        @test alloc_1 ≤ 32
        @test alloc_2 ≤ 688
        @test alloc_3 ≤ 944
        @test alloc_4 ≤ 144
        @test alloc_5 ≤ 1440
        @test alloc_6 ≤ 1632
    end
end


text:
#!/usr/bin/env julia
array1 = [(2n + 1)^2 for n in 1:5]
array2 = [sqrt(i) for i in array1]
println(typeof(1:5), " ",
        typeof(array1),
        " ", typeof(array2))
1:5, array1, array2


text:
__precompile__()
module JFVMM

import Base: +, -, *, /, .*, ./
export MeshStructure, BoundaryCondition, CellValue, FaceValue, CellVector,
       arithmeticMean, geometricMean, harmonicMean, upwindMean, linearMean,
       tvdMean, createBC, boundaryConditionTerm, cellBoundary, solvePDE,
       divergenceTerm, gradientTerm, convectionUpwindTerm, createCellVector,
       convectionTerm, convectionTvdTerm, diffusionTerm, createCellVariable,
       createFaceVariable, copyCell, fluxLimiter, createMesh1D,
       createMesh2D, createMesh3D, createMeshRadial2D, createMeshCylindrical2D,
       createMeshCylindrical3D, createMeshCylindrical1D, solveLinearPDE,
       linearSourceTerm, constantSourceTerm, transientTerm,
       solveMUMPSLinearPDE, faceEval, cellEval, permfieldlogrndg, permfieldlogrnde,
       # plot, imshow, xlabel, ylabel, figure, legend, pcolor, contour, colorbar,
       JFVM_test, solveExplicitPDE, reshapeCell,
       cellVolume, reshapeInternalCell, internalCells, domainInt, convectionTvdRHS

# Ignore visualization for now
# visualizeCellVectors
# visualizeCells

include("jfvmm_1types.jl")
include("jfvmm_2meshstructure.jl")
include("jfvmm_3boundarycondition.jl")
include("jfvmm_4domain_variables.jl")
include("jfvmm_5diffusion_terms.jl")
include("transientTerms.jl")
include("domainOperators.jl")
include("convectionTerms.jl")
include("averagingTerms.jl")
include("calculusTerms.jl")
include("sourceTerms.jl")
include("solveVisualizePDE.jl")
include("JFVMtools.jl")
include("jfvm_test.jl")

end # module


text:
abstract type AbstractModel end

"""
```
struct Model <: AbstractModel
    objective::Objective
    eq_constraints::VectorOfFunctions
    ineq_constraints::VectorOfFunctions
    box_min::AbstractVector
    box_max::AbstractVector
end
```

The `Model` structs stores information about the nonlinear optimization problem.
- `objective`: the objective function of the problem of type [`Objective`](@ref).
- `eq_constraints`: the equality constraints of the problem of type [`VectorOfFunctions`](@ref). Each function in `ineq_constraints` can be an instance of [`IneqConstraint`](@ref) or [`AbstractFunction`](@ref). If the function is not an `IneqConstraint`, its right-hand-side bound is assumed to be 0.
- `ineq_constraints`: the inequality constraints of the problem of type [`VectorOfFunctions`](@ref). Each function in `ineq_constraints` can be an instance of [`IneqConstraint`](@ref) or [`AbstractFunction`](@ref). If the function is not an `IneqConstraint`, its right-hand-side bound is assumed to be 0.
"""
mutable struct Model{Tv <: AbstractVector} <: AbstractModel
    objective::Union{Nothing, Objective}
    eq_constraints::VectorOfFunctions
    ineq_constraints::VectorOfFunctions
    box_min::Tv
    box_max::Tv
end

"""
    Model()

Constructs an empty model. The decision variables are assumed to be of type `Vector{Float64}`.
"""
Model() = Model(Float64, nothing)

"""
    Model(f)

Constructs an empty model with objective function `f`. `f` can be an instance of `Base.Function` but must return a number,  or it can be an intance of [`Objective`](@ref). The decision variables are assumed to be of type `Vector{Float64}`.
"""
Model(f::Function) = Model(Float64, f)

"""
    Model(::Type{T}, f::Union{Nothing, Function}) where {T}

Constructs an empty model with objective function `f` and decision variable value type `Vector{T}`. `f` can be an instance of `Base.Function` but must return a number,  or it can be an intance of [`Objective`](@ref). `f` can also be `nothing` in which case, the objective function can be defined later using [`set_objective!`](@ref).
"""
Model(::Type{T}, f::Function) where {T} = Model(T, Objective(f))
function Model(::Type{T}, obj::Union{Nothing, Objective}) where {T}
    return Model(obj, VectorOfFunctions(EqConstraint[]), VectorOfFunctions(IneqConstraint[]), T[], T[])
end

getobjective(m::AbstractModel) = m.objective

getineqconstraints(m::AbstractModel) = m.ineq_constraints

geteqconstraints(m::AbstractModel) = m.eq_constraints

function getobjectiveconstraints(m::AbstractModel)
    return VectorOfFunctions((getobjective(m), getineqconstraints(m), geteqconstraints(m)))
end

getineqconstraint(m::AbstractModel, i::Integer) = getineqconstraints(m).fs[i]

geteqconstraint(m::AbstractModel, i::Integer) = geteqconstraints(m).fs[i]

getnineqconstraints(m::AbstractModel) = length(getineqconstraints(m))

getneqconstraints(m::AbstractModel) = length(geteqconstraints(m))

getnconstraints(m::AbstractModel) = getnineqconstraints(m) + getneqconstraints(m)

getnvars(m::AbstractModel) = length(getmin(m))

"""
    getdim(m::AbstractModel)

Returns a 2-tuple of the number of constraints and the number of variables in the model `m`.
"""
getdim(m::AbstractModel) = (getnconstraints(m), getnvars(m))

get_objective_multiple(model::Model) = getobjective(model).multiple[]

function set_objective_multiple!(model::Model, m::Real)
    getobjective(model).multiple[] = m
    return model
end

getmin(m::Model)= m.box_min
getmin(m::AbstractModel, i::Integer) = getmin(m)[i]

getmax(m::Model) = m.box_max
getmax(m::AbstractModel, i::Integer) = getmax(m)[i]

function getinit(m::AbstractModel)
    ma = getmax(m)
    mi = getmin(m)
    return map(1:length(mi)) do i
        _ma = ma[i]
        _mi = mi[i]
        _ma == Inf && _mi == -Inf && return 0.0
        _ma == Inf && return _mi + 1.0
        _mi == -Inf && return _ma - 1.0
        return (_ma + _mi) / 2
    end
end

function setmin!(m::AbstractModel, min)
    getmin(m) .= min
    return m
end
function setmin!(m::AbstractModel, i::Integer, min)
    getmin(m)[i] = min
    return m
end

function setmax!(m::AbstractModel, max)
    getmax(m) .= max
    return m
end
function setmax!(m::AbstractModel, i::Integer, max)
    getmax(m)[i] = max
    return m
end

# Box constraints
function setbox!(m::AbstractModel, minb::T, maxb::T) where {T}
    getmin(m) .= minb
    getmax(m) .= maxb
    return m
end
function setbox!(m::AbstractModel, i::Integer, minb, maxb)
    getmin(m)[i] = minb
    getmax(m)[i] = maxb
    return m
end

function addvar!(m::Model, lb, ub)
    append!(getmin(m), lb)
    append!(getmax(m), ub)
    return m
end

"""
    set_objective!(m::AbstractModel, f::Function)

Sets the objective of the moodel `m` to the function `f`. `f` must return a scalar.
"""
function set_objective!(m::AbstractModel, f::Function)
    m.objective = Objective(f)
    return m
end

function add_ineq_constraint!(m::AbstractModel, f::Function, s = 0.0; dim = length(f(getinit(m))))
    return add_ineq_constraint!(m, FunctionWrapper(f, dim), s)
end
function add_ineq_constraint!(m::AbstractModel, f::AbstractFunction, s = 0.0)
    return add_ineq_constraint!(m, IneqConstraint(f, s))
end
function add_ineq_constraint!(m::Model, f::IneqConstraint)
    push!(m.ineq_constraints.fs, f)
    return m
end
function add_ineq_constraint!(m::Model, fs::Vector{<:IneqConstraint})
    append!(m.ineq_constraints.fs, fs)
    return m
end

function add_eq_constraint!(m::AbstractModel, f::Function, s = 0.0; dim = length(f(getinit(m))))
    return add_eq_constraint!(m, FunctionWrapper(f, dim), s)
end
function add_eq_constraint!(m::AbstractModel, f::AbstractFunction, s = 0.0)
    return add_eq_constraint!(m, EqConstraint(f, s))
end
function add_eq_constraint!(m::Model, f::EqConstraint)
    push!(m.eq_constraints.fs, f)
    return m
end
function add_eq_constraint!(m::Model, fs::Vector{<:EqConstraint})
    append!(m.eq_constraints.fs, fs)
    return m
end


text:
using PrettyTables

#                                   Formats
# ==============================================================================

data = [true  100.0 0x8080 "String"
        false 200.0 0x0808 "String"
        true  300.0 0x1986 "String"
        false 400.0 0x1987 "String"]

# Default
# =======

run(`cp example.tex format_default.tex`)
include_pt_in_file("format_default.tex", "Table 1", data, backend = Val(:latex),
                   backup_file = false, tf = tf_latex_default)
run(`pdflatex format_default.tex`)
run(`pdflatex format_default.tex`)
run(`convert -density 150 -trim -strip format_default.pdf format_default.png`)
run(`rm format_default.tex`)
run(`rm format_default.pdf`)

# Simple
# ======

run(`cp example.tex format_simple.tex`)
include_pt_in_file("format_simple.tex", "Table 1", data, backend = Val(:latex),
                   backup_file = false, tf = tf_latex_simple)
run(`pdflatex format_simple.tex`)
run(`pdflatex format_simple.tex`)
run(`convert -density 150 -trim -strip format_simple.pdf format_simple.png`)
run(`rm format_simple.tex`)
run(`rm format_simple.pdf`)

# Modern
# ======

run(`cp example.tex format_modern.tex`)
include_pt_in_file("format_modern.tex", "Table 1", data, backend = Val(:latex),
                   backup_file = false, tf = tf_latex_modern)
run(`pdflatex format_modern.tex`)
run(`pdflatex format_modern.tex`)
run(`convert -density 150 -trim -strip format_modern.pdf format_modern.png`)
run(`rm format_modern.tex`)
run(`rm format_modern.pdf`)

# Booktabs
# ========

run(`cp example_booktabs.tex format_booktabs.tex`)
include_pt_in_file("format_booktabs.tex", "Table 1", data, backend = Val(:latex),
                   backup_file = false, tf = tf_latex_booktabs)
run(`pdflatex format_booktabs.tex`)
run(`pdflatex format_booktabs.tex`)
run(`convert -density 150 -trim -strip format_booktabs.pdf format_booktabs.png`)
run(`rm format_booktabs.tex`)
run(`rm format_booktabs.pdf`)

#                                 Highlighters
# ==============================================================================

t = 0:1:20;

data = hcat(t, ones(length(t))*1, 1*t, 0.5.*t.^2);

header = ["Time" "Acceleration" "Velocity" "Distance";
           "[s]"  "[m/s\$^2\$]"    "[m/s]"      "[m]"];

hl_v = LatexHighlighter( (data,i,j)->(j == 3) && data[i,3] > 9, ["color{blue}","textbf"]);

hl_p = LatexHighlighter( (data,i,j)->(j == 4) && data[i,4] > 10, ["color{red}", "textbf"])

hl_e = LatexHighlighter( (data,i,j)->(i == 10), ["cellcolor{black}", "color{white}", "textbf"])

run(`cp example.tex latex_highlighter.tex`)
include_pt_in_file("latex_highlighter.tex", "Table 1", data, header,
                   backend = Val(:latex), backup_file = false,
                   highlighters = (hl_e, hl_p, hl_v))
run(`pdflatex latex_highlighter.tex`)
run(`pdflatex latex_highlighter.tex`)
run(`convert -density 150 -trim -strip latex_highlighter.pdf latex_highlighter.png`)
run(`rm latex_highlighter.tex`)
run(`rm latex_highlighter.pdf`)

run(`sh cleanall.sh`)


text:
groups = [split(group, "\n"; keepempty=false) for group in split(
  open("a6.txt") do file
    read(file, String)
  end, "\n\n")]

part1 = sum([length(union([Set(collect(s)) for s in group]...)) for group in groups])
println("$part1")

part2 = sum([length(intersect([Set(collect(s)) for s in group]...)) for group in groups])
println("$part2")


text:
module NormalMeanPrecisionNodeTest

using Test
using ReactiveMP
using Random

import ReactiveMP: make_node

@testset "NormalMeanPrecisionNode" begin

    @testset "Creation" begin
        node = make_node(NormalMeanPrecision)

        @test functionalform(node)          === NormalMeanPrecision
        @test sdtype(node)                  === Stochastic()
        @test name.(interfaces(node))       === (:out, :μ, :τ)
        @test factorisation(node)           === ((1, 2, 3), )
        @test localmarginalnames(node)      === (:out_μ_τ, )
        @test metadata(node)                === nothing

        node = make_node(NormalMeanPrecision, meta = 1)

        @test functionalform(node)          === NormalMeanPrecision
        @test sdtype(node)                  === Stochastic()
        @test name.(interfaces(node))       === (:out, :μ, :τ)
        @test factorisation(node)           === ((1, 2, 3), )
        @test localmarginalnames(node)      === (:out_μ_τ, )
        @test metadata(node)                === 1

        node = make_node(NormalMeanPrecision, factorisation = ((1, ), (2, 3)))

        @test functionalform(node)          === NormalMeanPrecision
        @test sdtype(node)                  === Stochastic()
        @test name.(interfaces(node))       === (:out, :μ, :τ)
        @test factorisation(node)           === ((1, ), (2, 3, ))
        @test localmarginalnames(node)      === (:out, :μ_τ, )
        @test metadata(node)                === nothing
    end

    @testset "AverageEnergy" begin

        begin
            q_out = PointMass(1.0)
            q_μ   = NormalMeanPrecision(1.0, 1.0)
            q_τ   = GammaShapeRate(1.5, 1.5)

            for N in (NormalMeanPrecision, NormalMeanVariance, NormalWeightedMeanPrecision), G in (GammaShapeRate, GammaShapeScale)
                marginals = (Marginal(q_out, false, false), Marginal(convert(N, q_μ), false, false), Marginal(convert(G, q_τ), false, false))
                @test score(AverageEnergy(), NormalMeanPrecision, Val{(:out, :μ, :τ)}, marginals, nothing) ≈ 1.6034261002694663
            end
        end

        begin
            q_out = NormalMeanPrecision(1.0, 1.0)
            q_μ   = NormalMeanPrecision(1.0, 1.0)
            q_τ   = GammaShapeRate(1.5, 1.5)

            for N in (NormalMeanPrecision, NormalMeanVariance, NormalWeightedMeanPrecision), G in (GammaShapeRate, GammaShapeScale)
                marginals = (Marginal(q_out, false, false), Marginal(convert(N, q_μ), false, false), Marginal(convert(G, q_τ), false, false))
                @test score(AverageEnergy(), NormalMeanPrecision, Val{(:out, :μ, :τ)}, marginals, nothing) ≈ 2.1034261002694663
            end
        end

        begin
            q_out = PointMass(0.956629)
            q_μ   = NormalMeanPrecision(0.255332, 0.762870)
            q_τ   = GammaShapeRate(0.93037, 0.79312)

            for N in (NormalMeanPrecision, NormalMeanVariance, NormalWeightedMeanPrecision), G in (GammaShapeRate, GammaShapeScale)
                marginals = (Marginal(q_out, false, false), Marginal(convert(N, q_μ), false, false), Marginal(convert(G, q_τ), false, false))
                @test score(AverageEnergy(), NormalMeanPrecision, Val{(:out, :μ, :τ)}, marginals, nothing) ≈ 2.209338084063204
            end
        end

        begin
            q_out = NormalMeanPrecision(0.148725, 0.483501)
            q_μ   = NormalMeanPrecision(0.992776, 0.545851)
            q_τ   = GammaShapeRate(0.309396, 0.343814)

            for N in (NormalMeanPrecision, NormalMeanVariance, NormalWeightedMeanPrecision), G in (GammaShapeRate, GammaShapeScale)
                marginals = (Marginal(q_out, false, false), Marginal(convert(N, q_μ), false, false), Marginal(convert(G, q_τ), false, false))
                @test score(AverageEnergy(), NormalMeanPrecision, Val{(:out, :μ, :τ)}, marginals, nothing) ≈ 4.155913074139921
            end
        end

        begin
            q_out_μ = MvNormalMeanPrecision([0.2818402997601115, 0.0847764277628964], [0.7059042678955475 0.3595204552322394; 0.3595204552322394 0.22068491824258746])
            q_τ     = GammaShapeRate(0.49074414, 0.4071772)

            for N in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), G in (GammaShapeRate, GammaShapeScale)
                marginals = (Marginal(convert(N, q_out_μ), false, false), Marginal(convert(G, q_τ), false, false))
                @test score(AverageEnergy(), NormalMeanPrecision, Val{(:out_μ, :τ)}, marginals, nothing) ≈ 38.88138702883774
            end
        end

        begin
            q_out_μ = MvNormalMeanPrecision([0.8378350736808462, 0.41494396892699026], [1.0074451742986652 0.4369298270351709; 0.4369298270351709 0.19572138039218784])
            q_τ     = GammaShapeRate(0.435485, 0.5269575)

            for N in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), G in (GammaShapeRate, GammaShapeScale)
                marginals = (Marginal(convert(N, q_out_μ), false, false), Marginal(convert(G, q_τ), false, false))
                @test score(AverageEnergy(), NormalMeanPrecision, Val{(:out_μ, :τ)}, marginals, nothing) ≈ 138.6947657738283
            end
        end

    end

end
end


text:
"""
Generate figure 1 which shows the two parameter gibbs toy model
===================================================== """

using PyPlot
srcpath  = "/Users/ethananderes/Dropbox/BayesLense/src/"
savepath = "/Users/ethananderes/Dropbox/BayesLense/paper/"
seed = 1
srand(seed)

"""
Set up the Gibbs chain
------------------------------------------"""
rho = -0.99
sigma1 = [1 rho; rho 1]
sigmainv1 = inv(sigma1)
sigmainv2 = [1 0;-1 1]* sigmainv1 * [1 -1;0 1]
sigma2 = inv(sigmainv2)
pdfimag1(x) = (1/ (2 * pi)) * sqrt(det(sigmainv1)) * exp(-0.5 * (transpose(x) * sigmainv1 * x)[1])
pdfimag2(x) = (1/ (2 * pi)) * sqrt(det(sigmainv2)) * exp(-0.5 * (transpose(x) * sigmainv2 * x)[1])

xlim = 2.0
space = 0.01
gridpoints = Array{Float64,1}[[x,y] for x=(-xlim):space:xlim, y=(-xlim):space:xlim] #|> transpose |> flipud
imag1 = map(pdfimag1, gridpoints)
imag2 = map(pdfimag2, gridpoints)

sum(imag1 * ^(space,2))
function findContours(density, gridarea, probvec)
	hts = Float64[]
	prb = Float64[]
	for dh in density[1:50:end]
		push!(hts, dh)
		push!(prb, sum( gridarea*density[density .> dh] ))
	end
	answr = Float64[]
	for k in probvec
		toss, index = findmin(abs(prb .- k))
		push!(answr, hts[index])
	end
	answr
end

"""
Run the chain for 20 steps
------------------------------------------"""
thetachain1 = [0.0]
phichain1 = [0.0]
for k = 1:20
  push!(thetachain1, sigma1[1,2] * phichain1[end] * sqrt(sigma1[1,1] / sigma1[2,2]) + randn() * sqrt(1-sigma1[1,2]^2) * sqrt(sigma1[1,1]) )
  push!(phichain1, phichain1[end])
  push!(phichain1,  sigma1[1,2] * thetachain1[end] * sqrt(sigma1[2,2] / sigma1[1,1]) + randn() * sqrt(1-sigma1[1,2]^2) * sqrt(sigma1[2,2]) )
  push!(thetachain1, thetachain1[end])
end

tildethetachain2 = [0.0]
phichain2 = [0.0]
for k = 1:20
  push!(tildethetachain2, sigma2[1,2] * phichain2[end] * sqrt(sigma2[1,1] / sigma2[2,2]) + randn() * sqrt(1-sigma2[1,2]^2) * sqrt(sigma2[1,1]) )
  push!(phichain2, phichain2[end])
  push!(phichain2,  sigma2[1,2] * tildethetachain2[end] * sqrt(sigma2[2,2] / sigma2[1,1]) + randn() * sqrt(1-sigma2[1,2]^2) * sqrt(sigma2[2,2]) )
  push!(tildethetachain2, tildethetachain2[end])
end

probvec = [.5, .85]
cnt1 =  findContours(imag1, space^2, probvec)
cnt2 =  findContours(imag2, space^2, probvec)

"""
figure1a.pdf
------------------------------------------"""
contour(imag1,
	levels = cnt1,
	origin= "lower",
	extent=[-2.5,2.5,-2.5,2.5],
	colors="k",
	linewidth=1.5,
	label=L"$P(t, \varphi|data)$"
	)
plot(phichain1,thetachain1, linewidth=1.5)
ax = gca()
ax[:spines]["right"][:set_color]("none")
ax[:spines]["top"][:set_color]("none")
ax[:xaxis][:set_ticks_position]("bottom")
ax[:spines]["bottom"][:set_position](("data",0))
ax[:yaxis][:set_ticks_position]("left")
ax[:spines]["left"][:set_position](("data",0))
ax[:text](0.6, 0.6, L"$P$ $(t, \varphi|data)$", fontsize=27)
xticks([-2,-1,1,2], fontsize = 15)
yticks([-2,-1,1,2], fontsize = 15)
savefig(joinpath(savepath, "figure1a.pdf"), dpi=300, bbox_inches="tight", pad_inches=0)

"""
figure1b.pdf
------------------------------------------"""
contour(imag2, levels = cnt2, origin= "lower", extent=[-2.5,2.5,-2.5,2.5], colors="k", linewidth=1.5)
xlabel(L"$\varphi$", fontsize = 28),
ylabel(L"$\widetilde{t}$", fontsize = 28)
plot(phichain2, tildethetachain2, linewidth=1.5)
ax = gca()
ax[:spines]["right"][:set_color]("none")
ax[:spines]["top"][:set_color]("none")
ax[:xaxis][:set_ticks_position]("bottom")
ax[:spines]["bottom"][:set_position](("data",0))
ax[:yaxis][:set_ticks_position]("left")
ax[:spines]["left"][:set_position](("data",0))
ax[:text](0.7, 0.7, L"$P$ $(\widetilde{t}, \varphi|data)$", fontsize=27)
xticks([-2,-1,1,2], fontsize = 15)
yticks([-2,-1,1,2], fontsize = 15)
savefig(joinpath(savepath, "figure1b.pdf"), dpi=300, bbox_inches="tight", pad_inches=0)

max_stars_repo_name: EthanAnderes/BayesianCmbLensing

max_issues_repo_name: EthanAnderes/BayesianCmbLensing

max_forks_repo_name: EthanAnderes/BayesianCmbLensing


text:
# Factorize two different matrices (one real, one complex) and solve for multiple rhs and free memory
using Distributed
if nworkers()<3
    addprocs(3)
end
@everywhere using MUMPSjInv
using Test

include("getDivGrad.jl");
A = getDivGrad(24,20,22);
A2 = getDivGrad(30,31,33);
n = size(A,1);
n2 = size(A2,1);
A  = A + im*spdiagm(0 => rand(n));
nrhs = 10;
rhs = randn(n,nrhs) + im*randn(n,nrhs);
rhs2 = randn(n2,nrhs);

@test nworkers() >= 3
println("Factorize complex matrix on worker 1")
@time F1 = remotecall_fetch(factorMUMPS,workers()[1], A, 1)
dump(F1)

println("Factorize real matrix on worker 2")
@time F2 = remotecall_fetch(factorMUMPS, workers()[2], A2, 1)
dump(F2)

println("Solve complex system on worker 1")
@time x = remotecall_fetch(applyMUMPS, workers()[3], F1, rhs)
err = zeros(nrhs)
for i=1:nrhs
    err[i] =  norm(A * x[:, i]-rhs[:, i]) / norm(rhs[:, i])
end
@test maximum(err) < 1e-14

println("Solve real system on worker 2")
@time x2 = remotecall_fetch(applyMUMPS, workers()[2], F2, rhs2)
err = zeros(nrhs)
for i = 1:nrhs
    err[i] = norm(A2 * x2[:, i] - rhs2[:, i]) / norm(rhs2[:, i])
end
@test maximum(err) < 1e-14

println("Free memory on worker 1")
destroyMUMPS(F1)

println("Free memory on worker 2")
destroyMUMPS(F2)


text:
"""
I/O routines for Dark Energy Camera (DECam) images from the Dark Energy
Camera Legacy Survey (DECaLS).
"""
module DECALSIO

import FITSIO
import ..Celeste: SurveyDataSet, BoundingBox, load_images

struct DECALSDataSet <: SurveyDataSet
basedir::String
"FITS file with CCD metadata, relative to base dir. E.g., survey-ccds-dr5.kd.fits"
metadatafile::String
end

function _linear_pix_to_world(crpix1, crpix2, crval1, crval2,
                              cd1_1, cd1_2, cd2_1, cd2_2, x, y)
    dx = x .- crpix1
    dy = y .- crpix2
    ra = crval1 + cd1_1 .* dx + cd1_2 .* dy
    dec = crval2 + cd2_1 .* dx + cd2_2 .* dy
    return ra, dec
end

function _get_overlapping_ccds(dataset::DECALSDataSet, box::BoundingBox)
    f = FITSIO.FITS(joinpath(dataset.basedir, dataset.metadatafile))
    hdu = f[2]::FITSIO.TableHDU
    crpix1 = read(hdu, "crpix1")::Vector{Float32}
    crpix2 = read(hdu, "crpix2")::Vector{Float32}
    crval1 = read(hdu, "crval1")::Vector{Float64}
    crval2 = read(hdu, "crval2")::Vector{Float64}
    cd1_1 = read(hdu, "cd1_1")::Vector{Float32}
    cd1_2 = read(hdu, "cd1_2")::Vector{Float32}
    cd2_1 = read(hdu, "cd2_1")::Vector{Float32}
    cd2_2 = read(hdu, "cd2_2")::Vector{Float32}
    width = read(hdu, "width")::Vector{Int16}
    height = read(hdu, "height")::Vector{Int16}
    close(f)

    function _linear_pix_to_world(x, y)
        dx = x .- crpix1
        dy = y .- crpix2
        ra = crval1 .+ cd1_1 .* dx .+ cd1_2 .* dy
        dec = crval2 .+ cd2_1 .* dx .+ cd2_2 .* dy
        return ra, dec
    end

    # calculate RA, Dec of all CCD corners using a linear WCS approximation
    ra1, dec1 = _linear_pix_to_world(1.0, 1.0)
    ra2, dec2 = _linear_pix_to_world(width, 1.0)
    ra3, dec3 = _linear_pix_to_world(1.0, height)
    ra4, dec4 = _linear_pix_to_world(width, height)

    # Calculate potential overlap by checking the minimum and maximum
    # extent of each CCD in RA/Dec
    # overlap the target bounding box.
    #
    # It would be better to do this by checking whether each line
    # segment of the CCD boundary intersects any line segment of the
    # bounding box. But this is generally good enough.
    ramin = min.(ra1, ra2, ra3, ra4)
    ramax = max.(ra1, ra2, ra3, ra4)
    decmin = min.(dec1, dec2, dec3, dec4)
    decmax = max.(dec1, dec2, dec3, dec4)

    # We need to treat the discontinuity at RA = 0. We will use a
    # "good-enough", but not complete, solution. In particular, it doesn't
    # account for singularities at the poles and it doesn't treat
    # extremely large (~180 degrees) boxes.  A more complete treatment
    # would use spherical geometry. See
    # https://github.com/spacetelescope/spherical_geometry
    # for an example implementation in Python.
    #
    # Because the corner RA/Dec values (ra1, ra2,...) were derived
    # from crval, crpix, we might have situations where ramin=-0.1,
    # ramax=0.1, or ramin=359.9, ramax=360.1, but fortunately we
    # *don't* need to worry about situations where ramin=0.1,
    # ramax=359.9. (erroneously swapping the east and west limits of
    # the CCD).
    #
    # We also know that the target box doesn't overlap the
    # discontinuity because we required that ramax > ramin when
    # constructing the box.
    #
    # Our strategy is simply to rotate everything away from the
    # discontinuity as far as possible. We find the angle that will
    # center the box at ra = 180, and then rotate all our CCD
    # ramin/ramax by the same angle.
    box_mean_ra = (box.ramax + box.ramin) / 2.0
    offset = 180.0 - box_mean_ra
    box = BoundingBox(box.ramin + offset, box.ramax + offset,
                      box.decmin, box.decmax)
    ramin .= (ramin .+ offset) .% 360.0
    ramax .= (ramax .+ offset) .% 360.0

    # The comparisons are a little unintuitive, because we're looking
    # for *any* overlap.
    mask = ((ramax .> box.ramin) .&
            (ramin .< box.ramax) .&
            (decmax .> box.decmin) .&
            (decmin .< box.decmax))

    # we could return an array of file names and an array of hdu's
    # (from `image_filename` and `image_hdu` columns), but for now just return
    # indices)
    return find(mask)
end

function load_images(dataset::DECALSDataSet, box::BoundingBox)

    # This could return filenames and hdus in which to find images overlapping
    # the box:
    idx = _get_overlapping_ccds(dataset, box)

    # Rest of this function not implemented yet. In order to finish loading
    # Celeste.Image objects from the DECals data, we need to:
    #
    # - Create a function to read in both raw and calibrated images and
    #   determine calibration array (simply divide the calibrated image
    #   by the raw image)
    #
    # - Change the Image type to accept a 2-d rather than 1-d calibration
    #   array (in general, calibration will not be constant in one
    #   dimension like in SDSS.)
    #
    # - Add `has_background` to Image (make background optional in Image)
    #   Then after calling `load_images(dataset, box)`, check `has_background`
    #   and if false, calculate the background with SEP. This allows us to
    #   optionally read the background from the dataset (e.g., SDSS supplies
    #   one, but DECaLS does not).
    #
    # - Add option to SDSSDataSet to read background or not (e.g., use the
    #   SDSS-calculated background or let Celeste determine the background).
    #
    # - Add general config options for background subtraction,
    #   detection parameters (not actually necessary for this function
    #   but nice to have)
    #
    # - Do something for the PSF (needs investigation into DECaLS PSF model).
    #
    # Testing: make smaller versions of some calibrated and raw images
    # images in that ra, dec range (remove ccds from files) and put
    # them at NERSC for download with a Makefile (Are the raw and
    # calibrated images actually on disk at NERSC?)

    error("Not yet implemented")
end

end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_mono_m4_opt_rho0_69(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_mono_m4_opt_rho0_69!(A_copy)
end

@inline function exp_mono_m4_opt_rho0_69!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=7
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 T2k7
    # Computing Bb2 = x*I+x*A
    coeff1=0.29569733081129235
    coeff2=0.5579019761039045
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.29569733081129235
    coeff2=0.5579019761039045
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-0.22527364173809233
    coeff2=0.12283893416894799
    coeff3=0.7635839995489598
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.8313380028114623
    coeff2=0.1346607765677363
    coeff3=0.06933288385478588
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.10603592380569533
    coeff2=0.05824938295333092
    coeff3=0.5674218115639013
    coeff4=1.1157677193291784
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=1.1108810217330913
    coeff2=0.25686278900572646
    coeff3=-0.04714911885726358
    coeff4=0.10219516402901102
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.022799236610704456
    coeff2=-1.0359744996702613
    coeff3=0.6837326299054575
    coeff4=0.4707893269262639
    coeff5=0.8155394030489035
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.015675034646166878
    coeff2=1.260909852925862
    coeff3=0.3148486342463852
    coeff4=0.3364746558547452
    coeff5=0.034247913899630375
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing T2k7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=1.0216444204245163
    coeff2=0.5251211761135689
    coeff3=0.7800244667714745
    coeff4=-0.13237659081921324
    coeff5=0.48088503410167177
    coeff6=0.0674717167283821
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    return memslots1 # Returning T2k7
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
#From TextAnalysis.jl Package.
function PrepCorpus(dirname::String,DocType::Type)
    # Recursive descent of directory
    # Add all non-hidden files to Corpus

    docs = {}

    function add_files(dirname::String)
        if !isdir(dirname)
            error("DirectoryCorpus() can only be called on directories")
        end

        starting_dir = pwd()

        cd(dirname)
        for filename in readdir(".")
            if isfile(filename) && !ismatch(r"^\.", filename)
                push!(docs, FileDocument(abspath(filename)))
            end
            if isdir(filename) && !islink(filename)
                add_files(filename)
            end
        end
        cd(starting_dir)
    end

    add_files(dirname)

    crps=Corpus(docs)

    standardize!(crps,DocType)
    return crps
    println(DocType)
    #return crps
end


text:
using Test

include("./julie_test.jl")

max_stars_repo_name: ErikQQY/Julie

max_issues_repo_name: ErikQQY/Julie

max_forks_repo_name: ErikQQY/Julie


text:
module SummitCheckpoint

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Summit Checkpoint" => Ahorn.EntityPlacement(
        Maple.Checkpoint
    )
)

baseSprite = "scenery/summitcheckpoints/base02.png"

function Ahorn.selection(entity::Maple.Checkpoint)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(baseSprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Checkpoint, room::Maple.Room)
    checkpointIndex = get(entity.data, "number", 0)
    digit1 = floor(Int, checkpointIndex % 100 / 10)
    digit2 = checkpointIndex % 10

    Ahorn.drawSprite(ctx, baseSprite, 0, 0)
    Ahorn.drawSprite(ctx, "scenery/summitcheckpoints/numberbg0$digit1.png", -2, 4)
    Ahorn.drawSprite(ctx, "scenery/summitcheckpoints/number0$digit1.png", -2, 4)
    Ahorn.drawSprite(ctx, "scenery/summitcheckpoints/numberbg0$digit2.png", 2, 4)
    Ahorn.drawSprite(ctx, "scenery/summitcheckpoints/number0$digit2.png", 2, 4)
end

end


text:
module Clusters

export create_uniform_cluster_ncube, create_uniform_cluster_nsphere, compute_centroid

using LinearAlgebra, Plots

"""
    create_uniform_cluster_ncube(center, perturbation, num_points)

Creates a cluster of size `num_points`, where each point is uniformly randomly
selected from the volume of an n-cube centered at `center`, aligned with the
standard basis, and with edge length 2 * `perturbation`

Returns a `num_points` * length(`center`) matrix, where each row is a point in
the cluster.
"""
function create_uniform_cluster_ncube(center::AbstractVector{<:Number},
                                      perturbation::Number,
                                      num_points::Integer)::AbstractMatrix{<:AbstractFloat}
    # (num_points, 1)
    ones_vec = ones(Integer, (num_points, 1))
    # (1, length(center))
    center_T = transpose(center)
    # (num_points, length(center)), where each row is a copy of center
    center_vecs = ones_vec * center_T

    # Generate uniformly random values between [0, perturbation) and negate some
    # randomly in order to avoid floating point precision errors and preserve
    # symmetry of the distribution.
    random_perturbations =
        perturbation * rand(Float64, (num_points, length(center)))
    random_flips = rand([-1, 1], (num_points, length(center)))

    return center_vecs + (random_perturbations .* random_flips)
end

"""
    create_uniform_cluster_nsphere(center, perturbation, num_points)

Creates a cluster of size `num_points`, centered around `center`, where each
point is uniformly randomly selected from the volume of an n-sphere with radius
'perturbation'.

Returns a `num_points` * length(`center`) matrix, where each row is a point in
the cluster.
"""
function create_uniform_cluster_nsphere(center::AbstractVector{<:Number},
                                        perturbation::Number,
                                        num_points::Integer)::AbstractMatrix{<:AbstractFloat}
    # (num_points, 1)
    ones_vec = ones(Integer, (num_points, 1))
    # (1, length(center))
    center_T = transpose(center)
    # (num_points, length(center)), where each row is a copy of center
    center_vecs = ones_vec * center_T

    # num_points * length(center) standard normal random variables
    rand_vecs = randn(Float64, (num_points, length(center)))
    # Take the norm along each column.
    # num_points
    norms = map(LinearAlgebra.norm,
                [rand_vecs[n,:] for n in 1:size(rand_vecs, 1)])
    # num_points uniform in [0, 1)
    uniform_coeffs = rand(Float64, num_points)

    origin_cluster = perturbation *
        Diagonal(uniform_coeffs .^ (1/length(center))) *
        Diagonal(1 ./ norms) *
        rand_vecs

    return center_vecs + origin_cluster
end

"""
    compute_centroid(points)

Computes the centroid of a matrix of points.

`points` is an m * n matrix that describes m points (rows) of dimension n
(columns).

Returns a single vector of size `n` corresponding to the centroid.
"""
function compute_centroid(
        points::AbstractMatrix{<:Number})::AbstractVector{<:AbstractFloat}
    point_sum = sum(points, dims=1)[:]
    return point_sum ./ size(points, 1)
end

"""
    plot_cluster(points)

Plots a matrix of points, where each row represents the coordinates of a single
point.
"""
function plot_cluster(points::AbstractMatrix{<:Number})
    plt = plot([points[:,n] for n in 1:size(points, 2)]...,
               seriestype = :scatter,
               showaxis = :show,
               aspect_ratio = :equal)
    display(plt)
end

end # module

max_stars_repo_name: brianxie/automaton-edification

max_issues_repo_name: brianxie/automaton-edification

max_forks_repo_name: brianxie/automaton-edification


text:
# Helper functions

function list_boards()
    numboards = get_num_boards()
    numboards == 0 && error("No boards found!")
    for n = 0:numboards-1
        setup_board(n)
        println("=> Board $n")
        model = is_AD12(n) ? "AD12" : is_AD14(n) ? "AD14" : is_AD16(n) ? "AD16" : "unknown"
        println("  model:      $model")
        println("  serial:     $(get_serial(n))")
        println("  frequency:  $(get_frequency(n)) MHz")
        println("  resolution: $(get_adcresolution(n))-bit")
        println("  memory:     $(get_memsize(n)) MiB")
        println("  channels:   $(get_all_channels(n))")
        println("  microsynth: $(has_microsynth(n) ? "yes" : "no")")
    end
end

function configure_for_ttl_triggering(boardnum,clock=1,channels=3;num_seg=0,seg_len=8192)
    setup_board(boardnum)
	set_clock(boardnum,clock)
	set_channels(boardnum,3) # needed to make sure clock is updated properly when moving to single channel mode
	set_channels(boardnum,channels)
	set_trigger(boardnum,4,1)
	# if num_seg == 0, segmented capture is disabled
	set_segmented_capture(boardnum,num_seg,seg_len)
end

function configure_for_waveform_triggering(boardnum,clock=1,channels=3,trig_ch=1,thresh=2^11,hyst=128;num_avg=0,avg_len=8192)
    setup_board(boardnum)
	set_clock(boardnum,clock)
	set_channels(boardnum,3) # needed to make sure clock is updated properly when moving to single channel mode
	set_channels(boardnum,channels)
	set_trigger(boardnum,1,1,trig_ch)
	set_waveform_trigger_params(boardnum,thresh,hyst)
	# if num_avg == 0, averaging is disabled
	set_averager(boardnum,num_avg,avg_len)
end

function init_board(boardnum::Int)
    # Force initialization and setup of the board
    clear_setupdone_bit(boardnum)
    setup_board(boardnum)
end

function get_volts_12(boardnum::Int, numblocks::Int, v_offset::T=0f0, v_conv::T=0.350f0) where T<:AbstractFloat
    # optimized method to get 12-bit voltage samples without allocating memory for the UInt16 data
    block_buffer = BLOCK_BUFFER[]
    block_samples = reinterpret(UInt16, block_buffer)
    len_samples = length(block_samples)
    volts = Vector{T}(undef, numblocks*len_samples)
    @inbounds for b in 1:numblocks
        # Transfer 1 block from the board into the buffer
        mem_read(boardnum, block_buffer)
        # Save that block as voltage samples
        @simd for i in 1:len_samples
            volts[(b-1)*len_samples+i] = T(block_samples[i]&0x0fff) * (v_conv/2^12) - (v_conv/2 + v_offset)
        end
    end
    return volts
end

function get_volts_16(boardnum::Int, numblocks::Int, v_offset::T=0f0, v_conv::T=0.350f0) where T<:AbstractFloat
    # optimized method to get 16-bit voltage samples without allocating memory for the UInt16 data
    block_buffer = BLOCK_BUFFER[]
    block_samples = reinterpret(UInt16, block_buffer)
    len_samples = length(block_samples)
    volts = Vector{T}(undef, numblocks*len_samples)
    @inbounds for b in 1:numblocks
        # Transfer 1 block from the board into the buffer
        mem_read(boardnum, block_buffer)
        # Save that block as voltage samples
        @simd for i in 1:len_samples
            volts[(b-1)*len_samples+i] = T(block_samples[i]) * (v_conv/2^16) - (v_conv/2 + v_offset)
        end
    end
    return volts
end

function get_blocks(boardnum::Int, numblocks::Int)
    blocks = Vector{UInt8}(undef, numblocks*DIG_BLOCK_SIZE)
    return get_blocks!(blocks, boardnum, numblocks)
end

function get_blocks!(blocks::AbstractVector{UInt8}, boardnum::Int, numblocks::Int)
    checkbounds(blocks,numblocks*DIG_BLOCK_SIZE)
    block_buffer = BLOCK_BUFFER[]
    @inbounds for b = 1:numblocks
        # Transfer 1 block from the board into the buffer
        mem_read(boardnum,block_buffer)
        # Save that block
        blocks[(b-1)*DIG_BLOCK_SIZE+1:b*DIG_BLOCK_SIZE] = block_buffer
    end
    return blocks
end

function get_samples_12(boardnum::Int, numblocks::Int)
    # To get the 12-bit samples we use get_samples_16 and keep only the 12 least significant bits
    samples = get_samples_16(boardnum, numblocks)
    samples .= samples .& 0x0fff
    return samples
end

function get_samples_16(boardnum::Int, numblocks::Int)
    samples = Vector{UInt16}(undef, numblocks*DIG_BLOCK_SIZE÷2)
    get_blocks!(reinterpret(UInt8, samples), boardnum, numblocks)
    return samples
end

function get_samples_32(boardnum::Int, numblocks::Int)
    samples = Vector{UInt32}(undef, numblocks*DIG_BLOCK_SIZE÷4)
    get_blocks!(reinterpret(UInt8, samples), boardnum, numblocks)
    return samples
end


text:
#
# OIFITS.jl -
#
# Support for OI-FITS (Optical Interferometry data format) in Julia.
#
#------------------------------------------------------------------------------

module OIFITS

export
    OIDataSet,
    OIDataBlock,
    OITargetEntry,
    OI_ARRAY,
    OI_CORR,
    OI_FLUX,
    OI_INSPOL,
    OI_T3,
    OI_TARGET,
    OI_VIS,
    OI_VIS2,
    OI_WAVELENGTH

import Base:
    IndexStyle,
    axes,
    convert,
    copy,
    eachindex,
    eltype,
    empty!,
    firstindex,
    get,
    getindex,
    getproperty,
    haskey,
    isempty,
    iterate,
    keys,
    lastindex,
    length,
    merge,
    merge!,
    ndims,
    propertynames,
    push!,
    read,
    read!,
    setindex!,
    setproperty!,
    show,
    size,
    take!,
    values,
    write

using Base: @propagate_inbounds

import FITSIO
import CFITSIO
using FITSIO: FITS, HDU, TableHDU, FITSHeader

include("types.jl")
include("formats.jl")
import .Formats: FieldDefinition, get_format
include("methods.jl")
include("io.jl")

end # module


text:
using RobotDynamics: state, control, states, controls

@testset "Trajectories" begin
n,m,N = 5,3,11
nx = fill(n,N)
nu = fill(m,N)
t,dt = 0.0,0.1
x = @SVector rand(n)
u = @SVector rand(m)
z = RD.KnotPoint(x,u,t,dt)

#--- Empty contructor
Z = RD.SampledTrajectory(n,m,dt=dt,N=N)
@test all(isnan, state(Z[rand(1:N)]))
@test all(isnan, control(Z[rand(1:N)]))
@test length(controls(Z)) == N-1
@test length(states(Z)) == N
@test eltype(states(Z)) <: SubArray
@test eltype(controls(Z)) <: SubArray
@test RD.vectype(Z) <: Vector
@test Z[1].dt ≈ dt
@test Z[end].dt ≈ 0
@test Z[1].t ≈ 0
@test Z[N].t ≈ (N-1)*dt
@test RobotDynamics.is_terminal(Z[end]) == true
@test RobotDynamics.dims(Z) == (nx,nu,N)
@test RobotDynamics.num_vars(Z) == N*(n+m)-m

Z = RD.SampledTrajectory{n,m}(n,m, dt=dt, N=N, equal=true)
@test Z isa RD.AbstractTrajectory
@test all(isnan, state(Z[rand(1:N)]))
@test all(isnan, control(Z[rand(1:N)]))
@test length(controls(Z)) == N
@test length(states(Z)) == N
@test eltype(states(Z)) <: SubArray
@test eltype(controls(Z)) <: SubArray
@test Z[1].dt ≈ dt
@test Z[1].t ≈ 0
@test Z[N].t ≈ (N-1)*dt
@test RobotDynamics.is_terminal(Z[end]) == false
@test RobotDynamics.num_vars(Z) == N*(n+m)
@test RD.dims(Z) == (nx,nu,N)

#--- Try different time inputs
tf = (N - 1) * dt
Z = RD.SampledTrajectory(n, m, tf=tf, dt=dt)
@test Z[end].t ≈ tf
@test length(Z) == N
@test Z[1].t == 0.0
@test Z[1].dt ≈ dt

Z = RD.SampledTrajectory(n, m, tf=2tf, dt=dt)
@test Z[end].t ≈ 2tf
@test length(Z) == 2N - 1
@test Z[1].t == 0.0
@test Z[1].dt ≈ dt

Z = RD.SampledTrajectory(n, m, dt=dt, tf=tf)
@test length(Z) == N

@test_throws ErrorException RD.SampledTrajectory(n, m)
@test_throws ErrorException RD.SampledTrajectory(n, m, dt=dt)
@test_throws ErrorException RD.SampledTrajectory(n, m, tf=tf)
@test_throws ErrorException RD.SampledTrajectory(n, m, N=N)
@test_throws AssertionError RD.SampledTrajectory(n, m, dt=dt, N=N, tf=2tf)
@test_throws AssertionError RD.SampledTrajectory(n, m, dt=2dt, N=N, tf=tf)

#--- Vector constructor
X = [@SVector rand(n) for k = 1:N]
U = [@SVector rand(m) for k = 1:N-1]
Z = RD.SampledTrajectory(X,U, dt=fill(dt,N-1))
@test states(Z) ≈ X
@test controls(Z) ≈ U
@test RobotDynamics.gettimes(Z) ≈ range(0, length=N, step=dt)

dts = rand(N-1)
Z = RD.SampledTrajectory(X,U, dt=dts, t0=1.0)
@test RD.gettimes(Z) ≈ [1.0; cumsum(dts)]

@test_throws ErrorException RD.SampledTrajectory(X,U, N=2N, dt=0.1)
Z = RD.SampledTrajectory(X,U, dt=dt)
@test RobotDynamics.gettimes(Z) ≈ range(0, length=N, step=dt)

#--- Copying
Z2 = copy(Z)
RobotDynamics.setstate!(Z[1], x)
@test !(state(Z[1]) ≈ state(Z2[1]))
@test state(Z[2]) ≈ state(Z2[2])
X = 2 .* X
U = 3 .* U
X[1] = x
RobotDynamics.setstates!(Z2, X)
@test state(Z[1]) ≈ state(Z2[1])
@test !(control(Z2[2]) ≈ U[2])
RobotDynamics.setcontrols!(Z2, U)
@test control(Z2[2]) ≈ U[2]

X2 = hcat(3X...)
RobotDynamics.setstates!(Z2, X2)
@test state(Z2[4]) ≈ 3X[4]
U2 = hcat(2U...)
RobotDynamics.setcontrols!(Z2, U2)
@test control(Z2[2]) ≈ 2U[2]

RobotDynamics.setcontrols!(Z2, u)
@test control(Z2[3]) ≈ u

times = range(0,length=N,step=dt*2)
@test RobotDynamics.gettimes(Z) ≈ times ./2
RobotDynamics.settimes!(Z, times)
@test Z[end].t ≈ 2
@test Z[1].dt ≈ 2dt

# Test approx
RobotDynamics.settimes!(Z, times ./ 2)
@test !(Z ≈ Z2)
copyto!(Z, Z2)
@test Z[1] ≈ Z2[1]
@test Z ≈ Z2

@test RobotDynamics.is_terminal(Z[end]) == true
push!(U, 2*U[end])
Z = RD.SampledTrajectory(X,U, dt=fill(dt,N-1))
@test length(controls(Z)) == N
@test RobotDynamics.is_terminal(Z[end]) == false

#--- Test copyto!
Z0 = [RD.KnotPoint(3*X[k], 2*U[k], dt*(k-1), dt) for k = 1:N]
@test state.(Z0) ≈ 3 .* states(Z)
@test control.(Z0) ≈ 2 .* controls(Z)
copyto!(Z0, Z)
@test state.(Z0) ≈ states(Z)
@test control.(Z0) ≈ controls(Z)

Z2 = RD.SampledTrajectory(rand() .* X, rand() .* U, dt=dt)
@test !(states(Z) ≈ states(Z2))
@test !(controls(Z) ≈ controls(Z2))
copyto!(Z2, Z)
@test states(Z) ≈ states(Z2)
@test controls(Z) ≈ controls(Z2)

# Test iteration
Z = RD.SampledTrajectory(X,U, tf=tf)
@test Z[1] ≈ RD.KnotPoint(X[1],U[1],0.0,dt)
@test Z[end] ≈ RD.KnotPoint(X[end], U[end], dt, dt*(N-1))
@test Base.IndexStyle(Z) == IndexLinear()
Z_ = [z for z in Z]
@test Z_[1] === Z[1]
@test Base.IteratorSize(Z) == Base.HasLength()
@test Base.IteratorEltype(Z) == Base.HasEltype()

@test states(Z, 2) == [x[2] for x in X]
@test states(Z, 1:3) == [[x[i] for x in X] for i = 1:3]

#--- Test functions on trajectories
model = Cartpole()
dmodel = RD.DiscretizedDynamics{RD.RK4}(model)
n,m = RD.dims(model)
fVal = [@SVector zeros(n) for k = 1:N]
X = [@SVector rand(n) for k = 1:N]
U = [@SVector rand(m) for k = 1:N-1]
Z = RD.SampledTrajectory(X,U,dt=dt)

# Test shift fill
Z = RD.SampledTrajectory(X,U,dt=dt)
RobotDynamics.shift_fill!(Z)
@test state(Z[1]) ≈ X[2]
@test control(Z[2]) ≈ U[3]
@test Z[1].t ≈ dt
@test Z[end].t ≈ N*dt
@test Z[N].dt == 0
@test Z[N-1].dt ≈ dt

# # Test dynamics evaluation
# discrete_dynamics!(RK3, fVal, model, Z)
# @test fVal[1] ≈ discrete_dynamics(RK3, model, X[1], U[1], 0.0, dt)
# dyn = TO.DynamicsConstraint{RK3}(model, N)
# conval = TO.ConVal(n,m, dyn, 1:N-1)
# TO.evaluate!(conval, Z)
# @test conval.vals ≈ fVal[1:N-1] .- X[2:N]
# @test !(conval.vals ≈ [zeros(n) for k = 1:N-1])

# Test rollout
nx = fill(n,N)
nu = fill(m,N)
RD.rollout!(RD.StaticReturn(), dmodel, Z, X[1])
Zmut = RD.SampledTrajectory([RD.KnotPoint{n,m}(MVector{n+m}(z.z), z.t, z.dt) for z in Z])
RD.rollout!(RD.InPlace(), dmodel, Zmut, X[1])
@test Z ≈ Zmut

# Block constructors
Z = RD.SampledTrajectory{n,m}(randn(n,N), randn(m,N-1), tf=2.0)
@test RD.dims(Z) == (nx,nu,N)
@test Z[end].t == 2.0
@test all(z->z.dt ≈ 0.2, Z[1:end-1])

Z = RD.SampledTrajectory([randn(n) for k = 1:N], [randn(m) for k = 1:N], dt=0.1)
@test RD.dims(Z) == (nx,nu,N)
@test Z[end].t == 1.0
@test RD.gettimes(Z) ≈ range(0,tf,step=dt)

@test_throws ErrorException Z = RD.SampledTrajectory(randn(n,N), randn(m,N-1))
@test_throws AssertionError RD.SampledTrajectory(randn(n,N), randn(m,N-1), dt=0.1, tf=2.0)

Z2 = RD.SampledTrajectory(zeros(n,N), zeros(m,N), dt=0.01)
copyto!(Z2, Z)
@test RD.states(Z2) ≈ RD.states(Z)
@test RD.controls(Z2) ≈ RD.controls(Z)

Z = RD.SampledTrajectory([@SVector randn(n) for k = 1:N], [@SVector randn(m) for k = 1:N-1], dt=0.1)
@test RD.vectype(Z[1]) === SVector{n+m,Float64}
Z2 = RD.SampledTrajectory(zeros(n,N), zeros(m,N-1), dt=0.01)
copyto!(Z2, Z)
RD.gettimes(Z2) ≈ RD.gettimes(Z)

tf = RD.set_dt!(Z2, 0.01)
@test tf ≈ 0.1
@test RD.gettimes(Z2) ≈ range(0,step=0.01,length=N)

## Trajectory with different number of knot points
nx = [4,4,4,4,4, 4, 2,2,2,2,2]
nu = [2,2,2,2,4, 2, 1,1,1,1,1]
X = [randn(n) for n in nx]
U = [randn(m) for m in nu]
dt = 0.1
Z = SampledTrajectory(X, U, dt=dt)
@test Z isa SampledTrajectory{Any,Any}
@test RD.dims(Z) == (nx,nu,11)
@test RD.dims(Z[1]) == (4,2)
@test RD.dims(Z[end]) == (2,1)

end


text:
@doc raw"""
    AbstractAffineConnection

Abstract type for affine connections on a manifold.
"""
abstract type AbstractAffineConnection end

"""
    LeviCivitaConnection

The [Levi-Civita connection](https://en.wikipedia.org/wiki/Levi-Civita_connection) of a Riemannian manifold.
"""
struct LeviCivitaConnection <: AbstractAffineConnection end

struct MetricDecoratorType <: AbstractDecoratorType end

"""
    AbstractConnectionManifold{𝔽,M<:AbstractManifold{𝔽},G<:AbstractAffineConnection} <: AbstractDecoratorManifold{𝔽}

Equip an [`AbstractManifold`](@ref) explicitly with an [`AbstractAffineConnection`](@ref) `G`.

`AbstractConnectionManifold` is defined by values of [`christoffel_symbols_second`](@ref),
which is used for a default implementation of [`exp`](@ref), [`log`](@ref) and
[`vector_transport_to`](@ref). Closed-form formulae for particular connection manifolds may
be explicitly implemented when available.

An overview of basic properties of affine connection manifolds can be found in [^Pennec2020].

[^Pennec2020]:
    > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection
    > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image
    > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229.
    > doi: 10.1016/B978-0-12-814725-2.00012-1.
"""
abstract type AbstractConnectionManifold{𝔽} <:
              AbstractDecoratorManifold{𝔽,MetricDecoratorType} end

"""
    connection(M::AbstractManifold)

Get the connection (an object of a subtype of [`AbstractAffineConnection`](@ref))
of [`AbstractManifold`](@ref) `M`.
"""
connection(::AbstractManifold)

"""
    ConnectionManifold(M, C)

Decorate the [`AbstractManifold`](@ref) `M` with [`AbstractAffineConnection`](@ref) `C`.
"""
struct ConnectionManifold{𝔽,M<:AbstractManifold{𝔽},C<:AbstractAffineConnection} <:
       AbstractConnectionManifold{𝔽}
    manifold::M
    connection::C
end

@doc raw"""
    christoffel_symbols_second(
        M::AbstractManifold,
        p,
        B::AbstractBasis;
        backend::AbstractDiffBackend = diff_backend(),
    )

Compute the Christoffel symbols of the second kind in local coordinates of basis `B`.
For affine connection manifold the Christoffel symbols need to be explicitly implemented
while, for a [`MetricManifold`](@ref) they are computed as (in Einstein summation convention)

````math
Γ^{l}_{ij} = g^{kl} Γ_{ijk},
````

where ``Γ_{ijk}`` are the Christoffel symbols of the first kind
(see [`christoffel_symbols_first`](@ref)), and ``g^{kl}`` is the inverse of the local
representation of the metric tensor. The dimensions of the resulting multi-dimensional array
are ordered ``(l,i,j)``.
"""
christoffel_symbols_second(::AbstractManifold, ::Any, ::AbstractBasis)

@decorator_transparent_signature christoffel_symbols_second(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

@doc raw"""
    christoffel_symbols_second_jacobian(
        M::AbstractManifold,
        p,
        B::AbstractBasis;
        backend::AbstractDiffBackend = diff_backend(),
    )

Get partial derivatives of the Christoffel symbols of the second kind
for manifold `M` at `p` with respect to the coordinates of `B`, i.e.

```math
\frac{∂}{∂ p^l} Γ^{k}_{ij} = Γ^{k}_{ij,l}.
```

The dimensions of the resulting multi-dimensional array are ordered ``(i,j,k,l)``.
"""
christoffel_symbols_second_jacobian(::AbstractManifold, ::Any, B::AbstractBasis)
function christoffel_symbols_second_jacobian(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=diff_backend(),
)
    n = size(p, 1)
    ∂Γ = reshape(
        _jacobian(q -> christoffel_symbols_second(M, q, B; backend=backend), p, backend),
        n,
        n,
        n,
        n,
    )
    return ∂Γ
end
@decorator_transparent_signature christoffel_symbols_second_jacobian(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

"""
    connection(M::ConnectionManifold)

Return the connection associated with [`ConnectionManifold`](@ref) `M`.
"""
connection(M::ConnectionManifold) = M.connection

Base.copyto!(M::AbstractConnectionManifold, q, p) = copyto!(M.manifold, q, p)
Base.copyto!(M::AbstractConnectionManifold, Y, p, X) = copyto!(M.manifold, Y, p, X)

@doc raw"""
    exp(M::AbstractConnectionManifold, p, X)

Compute the exponential map on the [`AbstractConnectionManifold`](@ref) `M` equipped with
corresponding affine connection.

If `M` is a [`MetricManifold`](@ref) with a metric that was declared the default metric
using [`is_default_metric`](@ref), this method falls back to `exp(M, p, X)`.

Otherwise it numerically integrates the underlying ODE, see [`solve_exp_ode`](@ref).
Currently, the numerical integration is only accurate when using a single
coordinate chart that covers the entire manifold. This excludes coordinates
in an embedded space.
"""
exp(::AbstractConnectionManifold, ::Any...)

@decorator_transparent_fallback function exp!(M::AbstractConnectionManifold, q, p, X)
    tspan = (0.0, 1.0)
    A = get_default_atlas(M)
    i = get_chart_index(M, A, p)
    B = induced_basis(M, A, i, TangentSpace)
    sol = solve_exp_ode(M, p, X, tspan, B; dense=false, saveat=[1.0])
    n = length(p)
    return copyto!(q, sol.u[1][(n + 1):end])
end

"""
    gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend())

Compute the Gaussian curvature of the manifold `M` at the point `p` using basis `B`.
This is equal to half of the scalar Ricci curvature, see [`ricci_curvature`](@ref).
"""
gaussian_curvature(::AbstractManifold, ::Any, ::AbstractBasis)
function gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; kwargs...)
    return ricci_curvature(M, p, B; kwargs...) / 2
end
@decorator_transparent_signature gaussian_curvature(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

function injectivity_radius(M::AbstractConnectionManifold, p)
    return injectivity_radius(base_manifold(M), p)
end
function injectivity_radius(M::AbstractConnectionManifold, m::AbstractRetractionMethod)
    return injectivity_radius(base_manifold(M), m)
end
function injectivity_radius(M::AbstractConnectionManifold, m::ExponentialRetraction)
    return injectivity_radius(base_manifold(M), m)
end
function injectivity_radius(M::AbstractConnectionManifold, p, m::AbstractRetractionMethod)
    return injectivity_radius(base_manifold(M), p, m)
end
function injectivity_radius(M::AbstractConnectionManifold, p, m::ExponentialRetraction)
    return injectivity_radius(base_manifold(M), p, m)
end

"""
    ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend())

Compute the Ricci tensor, also known as the Ricci curvature tensor,
of the manifold `M` at the point `p` using basis `B`,
see [`https://en.wikipedia.org/wiki/Ricci_curvature#Introduction_and_local_definition`](https://en.wikipedia.org/wiki/Ricci_curvature#Introduction_and_local_definition).
"""
ricci_tensor(::AbstractManifold, ::Any, ::AbstractBasis)
function ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; kwargs...)
    R = riemann_tensor(M, p, B; kwargs...)
    n = size(R, 1)
    Ric = allocate(R, Size(n, n))
    @einsum Ric[i, j] = R[l, i, l, j]
    return Ric
end
@decorator_transparent_signature ricci_tensor(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

@doc raw"""
    riemann_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend=diff_backend())

Compute the Riemann tensor ``R^l_{ijk}``, also known as the Riemann curvature
tensor, at the point `p` in local coordinates defined by `B`. The dimensions of the
resulting multi-dimensional array are ordered ``(l,i,j,k)``.

The function uses the coordinate expression involving the second Christoffel symbol,
see [`https://en.wikipedia.org/wiki/Riemann_curvature_tensor#Coordinate_expression`](https://en.wikipedia.org/wiki/Riemann_curvature_tensor#Coordinate_expression)
for details.

# See also

[`christoffel_symbols_second`](@ref), [`christoffel_symbols_second_jacobian`](@ref)
"""
riemann_tensor(::AbstractManifold, ::Any, ::AbstractBasis)
function riemann_tensor(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=diff_backend(),
)
    n = size(p, 1)
    Γ = christoffel_symbols_second(M, p, B; backend=backend)
    ∂Γ = christoffel_symbols_second_jacobian(M, p, B; backend=backend) ./ n
    R = allocate(∂Γ, Size(n, n, n, n))
    @einsum R[l, i, j, k] =
        ∂Γ[l, i, k, j] - ∂Γ[l, i, j, k] + Γ[s, i, k] * Γ[l, s, j] - Γ[s, i, j] * Γ[l, s, k]
    return R
end
@decorator_transparent_signature riemann_tensor(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

@doc raw"""
    solve_exp_ode(
        M::AbstractConnectionManifold,
        p,
        X,
        tspan,
        B::AbstractBasis;
        backend::AbstractDiffBackend = diff_backend(),
        solver = AutoVern9(Rodas5()),
        kwargs...,
    )

Approximate the exponential map on the manifold over the provided timespan
assuming the default connection of the given manifold by solving the ordinary differential
equation

```math
\frac{d^2}{dt^2} p^k + Γ^k_{ij} \frac{d}{dt} p_i \frac{d}{dt} p_j = 0,
```

where ``Γ^k_{ij}`` are the Christoffel symbols of the second kind, and
the Einstein summation convention is assumed. The arguments `tspan` and
`solver` follow the `OrdinaryDiffEq` conventions. `kwargs...` specify keyword
arguments that will be passed to `OrdinaryDiffEq.solve`.

Currently, the numerical integration is only accurate when using a single
coordinate chart that covers the entire manifold. This excludes coordinates
in an embedded space.

!!! note
    This function only works when
    [OrdinaryDiffEq.jl](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl) is loaded with
    ```julia
    using OrdinaryDiffEq
    ```
"""
function solve_exp_ode(M, p, X, tspan, B::AbstractBasis; kwargs...)
    return error(
        "solve_exp_ode not implemented on $(typeof(M)) for point $(typeof(p)), vector $(typeof(X)), and timespan $(typeof(tspan)). For a suitable default, enter `using OrdinaryDiffEq` on Julia 1.1 or greater.",
    )
end

#
# Introduce transparency
# (a) new functions & other parents
for f in [
    christoffel_symbols_second_jacobian,
    exp,
    gaussian_curvature,
    flat,
    get_coordinates,
    get_vector,
    log,
    mean,
    median,
    project,
    ricci_tensor,
    riemann_tensor,
    sharp,
    vector_transport_along,
    vector_transport_direction,
    vector_transport_direction!, #since it has a default using _to!
    vector_transport_to,
]
    eval(
        quote
            function decorator_transparent_dispatch(
                ::typeof($f),
                ::AbstractConnectionManifold,
                args...,
            )
                return Val(:parent)
            end
        end,
    )
end

# (b) changes / intransparencies.
for f in [
    christoffel_symbols_second, # this is basic for connection manifolds but not for metric manifolds
    exp!,
    flat!,
    get_coordinates!,
    get_vector!,
    get_basis,
    inner,
    inverse_retract!,
    log!,
    mean!,
    median!,
    norm,
    project!,
    sharp!,
    retract!,
    vector_transport_along!,
    vector_transport_to!,
]
    eval(
        quote
            function decorator_transparent_dispatch(
                ::typeof($f),
                ::AbstractConnectionManifold,
                args...,
            )
                return Val(:intransparent)
            end
        end,
    )
end
# (c) special cases
function decorator_transparent_dispatch(
    ::typeof(exp!),
    M::AbstractConnectionManifold,
    q,
    p,
    X,
    t,
)
    return Val(:parent)
end
function decorator_transparent_dispatch(
    ::typeof(inverse_retract!),
    M::AbstractConnectionManifold,
    X,
    p,
    q,
    m::LogarithmicInverseRetraction,
)
    return Val(:parent)
end
function decorator_transparent_dispatch(
    ::typeof(retract!),
    M::AbstractConnectionManifold,
    q,
    p,
    X,
    m::ExponentialRetraction,
)
    return Val(:parent)
end


text:
"""
    ConfigTypes

Module containing ClimateMachine configuration types.
"""
module ConfigTypes

export ClimateMachineConfigType,
    AtmosLESConfigType, AtmosGCMConfigType, OceanBoxGCMConfigType

abstract type ClimateMachineConfigType end
struct AtmosLESConfigType <: ClimateMachineConfigType end
struct AtmosGCMConfigType <: ClimateMachineConfigType end
struct OceanBoxGCMConfigType <: ClimateMachineConfigType end

end


text:
#=
LineDiff(index = 0,
         text = "")
""
->
"hello"

+LineDiff(index = 0,
         text = "hello")

"hello"
->
"he"

-LineDiff(index = 2
         text = "llo")

=#

struct LineDiff
    index :: Int64
    text :: String
end

+(x::LineDiff, y::LineDiff) =

    function buffer(diff :: LineDiff)
end

function test_line_diff()
end


text:
import Pkg
Pkg.activate("./")

using PDMesh

x, v, y, vol = create(Cylinder(10.0, 3.0, 10.0), resolution=0.5, rand_=0.0)

write_data("./examples/ovito.data", x, vol, ones(Int64, size(x)[2]))

max_stars_repo_name: ravinderbhattoo/PDMesh

max_issues_repo_name: ravinderbhattoo/PDMesh

max_forks_repo_name: ravinderbhattoo/PDMesh


text:
using LeapSeconds
using Dates: DateTime, Month, year, month, day
using ERFA
using Test

@testset "Leap Seconds" begin
    msg = "UTC is not defined for dates before 1960-01-01."
    @test (@test_logs (:warn, msg) offset_tai_utc(DateTime(1959,1,1,))) == 0.0
    for dt = DateTime(1960,1,1):Month(1):DateTime(2018,12,1)
        y = year(dt)
        m = month(dt)
        d = day(dt)
        @test offset_tai_utc(dt) ≈ ERFA.dat(y, m, d, 0.0)
    end
end